Jean-Charles Quinton
email: jeancharles.quinton@loria.fr

Bernard Girau
email: bernard.girau@loria.fr

A sparse implementation of dynamic competition in continuous neural fields

This paper introduces a sparse implementation of the Continuum Neural Field Theory, promoting a trade-off in accuracy for higher computational efficiency and alleviated constraints on the underlying model. The sparse version reproduces the main properties of previous discrete 2D implementations, such as dynamic competition leading to localized focus activity or robustness to noise and distracters, with a much higher computational speed on standard computer architectures.

Introduction

Under adequate conditions, neural networks with lateral inhibition are able to maintain bubbles of activity in response to input excitation [START_REF] Wilson | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF][START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF][START_REF] Taylor | Neural bubble dynamics in two dimensions: Foundations[END_REF]. The lateral connectivity pattern classically takes the form of a difference of Gaussians function (DoG) with local excitation and large-scale inhibition. These developments fall under the Continuum Neural Field Theory (CNFT) and general convergence results have been reproduced extensively in experimental studies. Biologically inspired by the two dimensional topology of the cortical sheet, analyzed in terms of cortical maps and further decomposed in cortical columns [START_REF] Burnod | An adaptive neural network: The cerebral cortex[END_REF], simulations were limited to two dimensions and most of them use a discretized version of the continuous equations. Anyway, from a computational point of view, updating the activity of a regular mesh of units with massive lateral connectivity has a polynomial algorithmic complexity, making any attempt to simulate a high dimensional CNFT difficult.

The bubbles emerging from such networks, i.e. localized Gaussian-like patches of activity, are able to focus on spatiotemporally coherent stimuli, despite the presence of noise or distracters [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. Dynamic competition between units and their ex-citation by topologically organized inputs thus lead to the emergence of a robust attentional property. This property is quite generic and should be beneficial to any distributed system where decision or action is required.

A sparse implementation of the CNFT

For most experimental conditions, the focus map either rapidly converges to a global close-to-zero activity or to a set of distant bubbles. These bubbles must indeed be separated enough for the combined effect of their self-maintenance and input excitation to counterbalance the large-scale inhibition. In light of these observations, its seems reasonable to approximate the focus map by a sum of Gaussian fields, at least after the few steps needed for convergence. The key point of the model presented in this paper is to approximate any activity bubble by a Gaussian field, whose profile is determined by the Mexican hat function and CNFT parameters. The mesh of units used for each map in the standard discrete implementation is thus replaced by a set of Gaussian parameters (see Fig. 1).

As units generally rely on synaptic connections to modulate their activity and take fixed positions in topological neural networks, whether artificial or biological, the location parameter of the Gaussian distributions may similarly be constrained to take a finite set of values. This constraint is however not at all required by the socalled sparse implementation presented in the paper, which is spatially continuous. It could be anyway introduced to simulate the CNFT dynamics on Gas-nets models where Gaussian like diffusion patterns propagate activity between distant units [START_REF] Husbands | Better living through chemistry : Evolving gasnets for robot control[END_REF]. These remarks however suppose that the Gaussian fields actually map physical units, whereas they may simply model the activity over the field, independently of the substrate. where k is the number of fields.

Expected benefits

The goal of a sparse implementation of the CNFT is twofold. First, alleviating the constraints on the underlying network makes the CNFT compatible with any system where topologically organized distributed processes compete for decision or action. At the same time, a sparse implementation allows a dramatic speed up in the computations. Fast simulation of the network dynamics may be useful to easily tune the parameters of the model. Moreover, real-time processing of high dimensional networks is a requirement when manipulating sensorimotor systems. When controlling robots for instance, time constraints are inherent to the interaction loop and must be dealt with.

To further develop the potential benefits from manipulating compact representations, details of the perception-action loop must be stated. Both biological and artificial systems interact with their environment through a large number of sensory and motor dimensions. At least in artificial systems, many sensory inputs and motor commands take a unique and precise value (joint angles for example). When using cortical map representations, single values must be projected on discrete 2D maps, the activity at each unit of each map regularly updated, before synthesizing the unit activities on motor maps into a single command to be sent to effectors. Although this might be totally irrelevant for biological systems, using compact representations of the activity on the different maps not only avoids such conversions, but also allows the manipulation of high dimensional spaces, without using self-organizing maps to reduce the dimensionality [START_REF] Ménard | Model of multi-modal cortical processing: coherent learning in self-organizing modules[END_REF] and by limiting the combinatorial explosion of nD maps (with n > 2). At a theoretical level, turning to an implementation with continuous inputs and bubbles makes several aspects of the analysis closer to the case described in the original CNFT equations [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF]. A fixed spatial discretization is no longer required and the processing resolution can freely adapt to the task performed by the simulated system. Nevertheless, a fully continuous model apparently lacks the biological plausibility of discrete cortical maps with interconnected localized units. However, this might simply reflect a shift in focus from the underlying substrate to the map activity dynamics. In this perspective, it would remain useful as a practical tool, whatever its flaws or lack of accuracy for fine-grained models of cortical networks.

From discrete to sparse CNFT

In the following sections, we will adopt the notations introduced by Amari [1] and Rougier et al. [START_REF] Rougier | Emergence of attention within a neural population[END_REF], whose discrete implementation will serve as a reference.

Standard equations and discrete implementation

The focus neural field is represented by a manifold M in bijection with [-0.5, 0.5] 2 and the membrane potential at the position vector x and time t on this field by u(x,t). Similar notations are used for the input stimulation s(x,t). The dynamics of the membrane potential is described by the following single-layer field equation of lateral inhibition type:

τ ∂ u(x,t) ∂t = -u(x,t) + x ∈M w(x, x )u(x ,t)dx +s(x,t) + h ( 1
)
where h is the resting potential and w(x, x ) the lateral connection weight function satisfying the following equation:

w(x, x ) = Ae |x-x | 2 a 2 -Be |x-x | 2 b 2 (2)
In order to perform numerical simulations of the CNFT dynamics, the continuous equations have been discretized as follows:

τ ∂ u(x i j ,t) ∂t = -u(x i j ,t) + 1 n 2 ∑ (k,l)∈[0;n[ 2 w(x i j , x kl )u(x kl ,t) +s(x i j ,t) + h (3)
where x i j is the discrete position associated to the unit (i, j) and given by Eq. 4.

x i j = i n -0.5, j n -0.5 ∀(i, j) ∈ [0; n[ 2 (4)

Sparse implementation

To exclusively manipulate Gaussians in a sparse implementation, all terms in Eq. 1 must be transformed accordingly. At this point, it may be argued that stimulations s(x,t) will not necessarily take the form of normal distributions, especially when considering noisy inputs or overlapping stimuli. However, any discrete signal can be approximated by a sum of Gaussians, their number being positively correlated with the precision required [START_REF] Goshtasby | Curve fitting by a sum of gaussians[END_REF]. Although efficient algorithms to perform such a decomposition are available [START_REF] Childs | A fast, space-efficient algorithm for the approximation of images by an optimal sum of gaussians[END_REF], one goal of the current implementation is to make it compatible with sparse inputs.

In artificial systems, these inputs might be associated with bottom-up signals (coming from various measurement systems and modalities) as well as top-down signals competing for decision or action. In biological systems, sensory features are encoded by myriads of neurons, but population coding is often considered as it averages the variability of single neuron spike trains, even when each of them already takes inputs from large receptive fields [START_REF] Kandel | Principles of Neural Science[END_REF]. Receptive fields have been found to be

05 : =0 c ¢ 1 : =0 c ¢ 2 : =0 c ¢ 4 : =0 c ¢
Fig. 2 Comparison of the discrete and sparse algorithms on a few iterations with increasingly more distant Gaussian components (light gray). The discrete algorithm (dark gray) computes an output activity for each unit (1000 units uniformly sampled in [0; 1]). The sparse algorithm (black) computes the reciprocal inhibition between Gaussian distributions, merges them when close enough and outputs the sum of the resulting distributions. From left to right, the distance between the input centers ∆ c takes values in {0.05, 0.1, 0.2, 0.4}. The limit between excitation/merging and inhibition is 0.1 for both versions so that the maximal distortion appears for ∆ c 0.1.

approximately separable into a sum of amplitude modulated Gaussian components, for instance in the visual system [START_REF] Rodieck | Quantitative analysis of cat retinal ganglion cell response to visual stimuli[END_REF][START_REF] Wennekers | Separation of spatio-temporal receptive fields into sums of gaussian components[END_REF].

The effect of the CNFT computation on distributions whose centers are sufficiently separated relatively to the excitatory standard deviation a in Eq. 2 (∆ c >> a), is only inhibitory. Their interactions can then be well approximated by a point to point interaction, only considering activity distribution centers (as centers of mass are used in point mechanics). This can be proved using Taylor developments of the competition function w between positions x and x + dx with dx |x -x |. The same intuition is mathematically valid for ∆ c = 0, the sum of aligned normal distributions being a Gaussian distribution. Nevertheless, this approximation does not hold when ∆ c takes intermediate values. To evaluate the distortions resulting from such a rough approximation, numerical simulations have been performed. To get an idea of the potential quantitative differences between the discrete and sparse implementations for a range of ∆ c values, Fig. 2 compares the fields produced after a few iterations.

From now on, the following notations are adopted. A generic field G will be defined as the sum of k components g k . Determined by the set of parameters (x k , I k ), g k denotes a Gaussian field of amplitude I k centered on x k . Let g k (x) be the activity propagated by the Gaussian component g k at the point x satisfying Eq. 5.

g k (x) = I k * e |x k -x| 2 σ 2 (5)
with σ being equal to the excitatory standard deviation a of Eq. 2 for the focus field components. In most cases, this produces the stereotypical profile of the focus bubbles and is coherent with the choice of synthesizing each bubble by a single Gaussian field. Finally, the potential at any point of the field can then be computed as in Eq. 6. The choice of a field named G in this section is arbitrary and G could be

= u t 1 ∪ s t 1 ∪ c 1 )
, reflecting the convergence towards a single bubble of activity. substituted by any field introduced in the equations and paragraphs to come, such as the stimulation field S or focus field U. Let S t , U t and s t k , u t k be their respective state and components at time t.

g(x,t) = ∑ k g t k (x) (6)
To briefly introduce the algorithmic decomposition of Eq. 1 detailed in the following sections, we illustrate the expected dynamics of the model on Fig. 3. For each temporal step of the simulation, a competition field C is computed by propagating activity from the different components of the focus field U (see Sec. 2.3). The input, focus and competition fields are then integrated in a single field where spatiotemporally coherent stimuli form clusters of excitatory components (see Sec. 2.4). Activity resulting from the presence of noise or distracters (such as s t 2 on the figure) is not reinforced by focus components and counterbalanced by inhibitory competition components (c 3 ). Finally, close components are merged (u t+dt 1 ) and Gaussian with negative activity eliminated (see Sec. 2.5). In standard conditions, this ensures that a reduced set of distant bubbles will emerge, each represented by a single component.

Competition

The first step consists in computing an equivalent to the dynamic competition term in Eq. 1. For each focus and stimulation component, an activity I c k is computed by applying Eq. 7. For each u t k and s t k considered, an inhibitory component c k of parameters (x k , I c k ) is thus produced. All c k are combined in global inhibitory field C (for competition).

I c k = 1 n n ∑ i=1 w(x k , x i )I i ( 7
)
where n is the number of components on the focus field. The lateral competition weight function w has already been introduced in Eq. 2. Note that in the original equations, the dynamic competition was evaluated at each point of the field. In a sparse implementation, it is therefore required to account for the inhibition at any place where there is a Gaussian distribution. The stimulation S being independent of the internal computations performed, activity on the input field must also be inhibited (see Fig. 4). For example on Fig. 3,s t 2 would not be inhibited if c 3 was not introduced and only the interactions between u t 1 and u t 2 were taken into account.

Integration

The focus field U t , inhibitory field C resulting from the lateral competition and input field S must then be combined and integrated over time to reproduce the dynamics of Eq. 1. This integrative step is described by the following equation:

U t+dt = U t ∪ dt τ -U t ∪C ∪ S t + h (8)
where the ∪ operator applied to fields actually corresponds to the union of the component sets. When computing the potential at a given position by applying Eq. 6, this is actually equivalent to adding the contributions of the different fields. The scalar multiplication (by dt τ ) and addition (of h) are directly applied to the intensity of the Gaussian components. If we respectively denote by U and U the fields equivalent to -U t ∪C ∪ S t and dt τ (U + h), their component parameters (x k , I k ) and (x k , I k ) satisfy:

I k = dt τ (I k + h) (9)

Merging

Although the potential at any point of the field can be updated and computed at the end of the previous section, yet another step is required to account for the bubbles reinforcement and to avoid a combinatorial explosion. The unions in Eq. 8 indeed lead to an ever increasing number of components on the focus field. Even though in practice and for efficiency purpose, computations are factored so that the number of Gaussian fields centered on the same position remains minimal, stimulations are not constrained and constantly inject new components into the focus field. If the number of input Gaussians at t is n t , an underestimation of the number of components on the focus field at time T is ∑ T t=0 n t , which is quite problematic when considering the polynomial cost of computing the lateral competition. A merging algorithm is therefore introduced. Similar constraints appear in different research fields when manipulating large Gaussian mixture models or HMM-based models with continuous densities [START_REF] Xu | A new approach to merging gaussian densities in large vocabulary continuous speech recognition[END_REF]. The main difference being that the actual distribution of local stimulations is not here accessible to the model. Due to the robustness of the CNFT dynamics and efficiency constraint of the implementation, a simple Euclidian distance between the distribution locations is used as the merging criterion. This distance threshold is chosen to match the excitatory standard deviation a as to facilitate and reinforce the emergence of stereotyped bubbles of activity.

Algorithm 1 Merging algorithm for the Gaussian components on the focus field With the standard equations 1 and 3, bubbles track the moving stimuli on which they are focused because the renewed inputs introduce an asymmetry in the dynamic competition computations. The shifted bubble is in turn reinforced by the input, although always lagging behind the stimulation. The same tracking characteristic is here achieved by the merge function introduced in Algorithm 1 and defined by Eq. 10.

x new = I i I i + I j x i + I j I i + I j x j I new = I i + I j -I i × I j × |x i -x j | 2 α 2 (10)
where α is a constant higher than a. When components are perfectly aligned, as it happens when computing the inhibitory field I -from the focus field U, the Gaussian amplitudes are simply added during the merging operation, which is coherent with the discrete unit based computations of Eq. 3. When a focus component and a stimulation are close enough for the merging to occur, their intensities have respectively been scaled by 1 -dt τ and dt τ . This weighting allows the bubble to slightly shift towards the new stimulus while reinforcing or at least maintaining activity. Taking α = +∞ cancels the last term and introduces additional discontinuities but still allows the algorithm to perform correctly as long as a is chosen carefully. Components of negative amplitude at the end of the merging phase are simply removed from the field, thus limiting the focus potential to positive values.

Finally, when a single locus representing the entire focus field is needed, whether for undertaking actions or statistical analysis purpose, the barycenter c of the remaining components (i.e. those which have not been merged and represent distant bubbles of activity) can be rapidly computed by applying Eq. 11.

c = ∑ k I k × x k ∑ k I k (11)
3

Results and performance

The analysis of the model dynamics and its performance are relative to the competition and attentional properties of the CNFT. The discrete implementation described in [START_REF] Rougier | Emergence of attention within a neural population[END_REF] is used as a reference, and the same evaluation process is adopted. The authors commit to the assumption that a single bubble should emerge on the focus field and that it should robustly track the associated input stimulus despite noise and distracters. The scenarios posed by Rougier and Vitay guarantee that for adequate sets of parameters, there will be enough iterations for a single bubble to appear before noise and distracters are added. As a consequence, the performance is easily ascertained by measuring the distance between the input stimuli and focus bubble centers at any time.

Functional results

Although the discrete and sparse implementations share general principles and parameters, there are many qualitative differences in the algorithms, not to state the most obvious difference in the type of stimulation provided to the system. Since the two versions might not stand on equal footing, the goal here is mainly to illustrate how a Gaussian based implementation may reproduce the classical CNFT characteristics.

A fair comparison at least requires using optimal sets of parameters, which might greatly differ between implementations. To avoid complex hand-tuning, genetic algorithms are applied to evolve a population of parameter vectors independently for each version. The fitness function to minimize is on the cumulative error on a set of well chosen scenarios simulated in sequence. Indeed, depending on the task on which the optimization process is applied, optimal parameters may differ. For instance, the a value is highly correlated to the distance between the stimuli to focus on: a large value may improve the results only when unambiguous distant stimuli are presented. The parameters are thus optimized as to produce the best results for all three following scenarios: A) 2 bell-shaped distant stimuli s 1 and s 2 are introduced at time t = 0. Their intensity are governed by I 1 = 0.4 and I 2 (t) = 0.5 + 0.5 cos(π × (t/5)). B) 1 bell-shaped stimulus of standard deviation 0.1 and intensity 1.0 follow a circular trajectory of radius 0.2 around the point (0, 0) at 10 deg/s from t = 0. From t = 1, 5 distracters are added and take new random positions on the field every 1 s. C) 1 bell-shaped stimulus (same as above). At t = 1, Gaussian noise of amplitude 0.5 is added.

In the discrete version, noise is easily added by changing the stimulation intensity by a random amount for each x i j . This is not possible in the sparse version as each Gaussian component already corresponds to a correlated activity over the entire field. Adding a large number of random components would not do the trick as it would be quite equivalent to adding distracters. However, once a bubble is stabilized as it is guaranteed by the experiments reproduced, units outside the bubble are largely inhibited and the associated inputs have no effect on further computations. Noise introduces asymmetries under the bubble and destabilize it, which is roughly

Computational speed up

The following paragraphs do not take into account the cost decomposing a discrete input into a sum of Gaussian distributions and are thus mainly relevant for artificial systems where inputs are compatible, as argued in Sec. 1.2 and 2.2. A simple algorithmic complexity analysis helps understanding the computational benefits from the sparse implementation. Let n 2 be the size of the maps. Using infinite asymptotics, all terms are dominated by the convolution computed for the lateral competition, so that T (n) = O(n 4 ) using Landau notation1 . When using a singular value decomposition (SVD) on the lateral connection kernel, the 2D convolution can be decomposed in two orthogonal 1D convolution and the complexity drops to T (n) = O(n 3 ). The initial cost of computing the SVD can be neglected as long as the kernel remains constant throughout the entire simulation. When modeling higher dimensional maps of size n m with m > 2 and considering the non optimized version, the complexity becomes T (n) = O((n m ) 2 ) and projecting the data on the low dimensional space is required.

For the sparse implementation again, most computations can be neglected relatively to the inhibition and merging operations. The complexity of the lateral inhibition is in the worst case T c (n u ) = O(n u

2 ) if we let n u be the number of components on the focus field, supposed far greater than the number of input components. The complexity of the merging operation depends on the actual distribution of the components, and can therefore hardly be mathematically approximated. In practice, and as soon as a group of components get reinforced by the inputs, the lateral competition rapidly inhibits a large number of components and their number drastically decreases. In such a case, infinite asymptotics is no more valid but the computational cost of updating the fields anyhow drops far below the cost of the discrete version, that only depends on the size of the map. This is not surprising as making localized bubbles of activity emerge is the sole goal of the lateral competition, and was one reason for developing a sparse implementation. Experimental performance results are reproduced on Table 1 for different implementations and scenarios.

Conclusion

In this paper, we introduced a sparse version of the CNFT able to reproduce its main properties and speed up the computations with a relative independence to the number of dimensions manipulated. The strong assumptions and rough approximations introduced in the process must however be taken into account when applying it to in- puts not easily reduced to a sum of Gaussian components. Part of the problem might then be deviated to the input filtering and decomposition algorithms, instead of being treated directly by the CNFT. Perspectives include a fully distributed implementation for real-time high dimensional robotics applications, since the reduced number of components, their limited exchanges and low memory requirements make it compatible with parallel hardware that provides fast paced computations but limited resources.

Fig. 1

 1 Fig. 1 Focus map representations when two bubbles have emerged. (a) In the discrete version, as a regular grid of units with activity (a i j ) ∈ [0; 1] n×n where n is the size of the grid. (b) In the sparse version, as the sum of Gaussian fields g k defined by the parameters (x k , I k ) ∈ ([-0.5; 0.5] 2 × [0; 1]) kwhere k is the number of fields.

Fig. 3

 3 Fig. 3 Algorithmic decomposition of the sparse implementation.1 A competition field C is produced by propagating activity from {u t 1 , u t 2 } to {u t 1 , u t 2 , s t 1 , s t 2 }. 2 Components from the focus field {u t 1 , u t 2 }, input field {s t 1 , s t 2 } and competition field {c 1 , c 2 , c 3 } are integrated. 3 Close components are merged, and resulting components with negative intensity removed (s t 2 ∪ c 3 and u t 2 ∪ c 2 ). Only one component remains (u t+dt

Fig. 4

 4 Fig. 4 Competition between components. (a) The inhibitory activity I ck is computed at x k for the input field s k (supposing there is only one Gaussian field component u i on the focus). (b) The activity at x k is reduced. If I c k had not been taken into account (dashed line), s k would not have been negligible when computing u t+dt although far away from the only active bubble on the focus u t i .

Fig. 5

 5 Fig.5Functional results. (a) Normalized error dynamics for scenario A. Although the shifts in attention associated with a high error occur at slightly different times, the same kind of hysteresis appears for the sparse (dashed) and discrete (plain) versions relatively to s 2 intensity (cosine function). (b) Although the error is lower for the discrete version, both implementations remain below 2% for the scenarios A to C.

  Components from the focus field {u t 1 , u t 2 }, input field {s t 1 , s t 2 } and competition field {c 1 , c 2 , c 3 } are integrated. 3 Close components are merged, and resulting components with negative intensity removed (s t 2 ∪ c 3 and u t 2 ∪ c 2 ).

	Only one component remains
	(u t+dt 1

1 A competition field C is produced by propagating activity from {u t 1 , u t 2 } to {u t 1 , u t 2 , s t 1 , s t 2 }. 2

Table 1

 1 Comparison of mean computation time (in µs) to update the model for one timestep, exclusive of the input generation. Although there is a bigger variance with the sparse implementation since the computations depend on the number and distribution of components, the time remains over 50 times lower than for the SVD optimization of the discrete version of the model.

	Scenario	Discrete version	Discrete (SVD)	Sparse version
	A (alternation)	321131	37537	596
	B (distracters)	321084	37724	953
	C (noise)	321003	37491	632

Although n does not usually take very high values, a map of size 50 × 50 is common and already justifies the effect of n on performance.


Tom Leclerc
email: tom.leclerc@loria.fr

Laurent Ciarletta
email: laurent.ciarletta@loria.fr

André Schaff
email: andre.schaff@loria.fr

SLSF: Stable Linked Structure Flooding For Mobile Ad Hoc Networks

For some applications in ad hoc networks optimal dissemination is a key issue (e.g. service discovery, network management). In this paper, we are creating and exploiting stable (sub-)structures to achieve an efficient (as far as low network resource usage is concerned) dissemination by building a two-layer protocol. Firstly, single-hop clusters, among stableconnected devices, are created. Secondly, on top of those clusters, inter-cluster relays (ICR) are determined. This leads to an overall stable-connected structure. The results show that the proposed stable linked structure flooding (SLSF) protocol efficiently disseminates data among stable nodes. Interestingly with growing density both the number of forwarding nodes and the bandwidth used remain comparatively low. Therefore we plan to use SLSF as a basis for a stable service discovery.

I. INTRODUCTION

In this paper we consider large Mobile Ad hoc NETworks (MANET) where the wirelessly connected devices communicate spontaneously without any predefined infrastructure with each other [START_REF] Dousse | Connectivity in ad-hoc and hybrid networks[END_REF]. To reach a destination, nodes communicate using intermediate nodes as routers. Moreover, topology changes are common since joining and leaving of nodes occur dynamically. Our main research topic is service discovery in such networks where it is more important to provide stable services to contributing nodes than transient services to unstable nodes. Therefore optimal service discovery is more about a qualitative discovery (stable nodes/structures and their reachability) than a quantitative one (discovering all possible services). Firstly we build clusters (local groups) of one-hop stable-connected devices in a self-organizing manner using the NLWCA clustering protocol [START_REF] Andronache | Nlwca node and link weighted clustering algorithm for backbone-assisted mobile ad hoc networks[END_REF]. Secondly to create bigger stable-linked network structures, stable connections between nearby clusters are discovered. To do so we propose Inter-Cluster Relays (ICR) that are inspired by the Multi-Point Relays (MPR) of OLSR [START_REF]Optimized link state routing protocol (olsr), rfc3626[END_REF]. Moreover, extending WCPD [START_REF]Hytrace backbone-assisted path discovery in hybrid networks[END_REF], specific beacon formats are specified. We exploit the stablelinked structures within the network topology to streamline information exchange and to minimize the overhead.

A. NLWCA & WCPD

The Weighted Cluster-based Path Discovery protocol (WCPD) is designed to take advantage of the cluster topology built by the Node and Link Weighted Clustering Algorithm (NLWCA) in order to provide path discovery and broadcast mechanisms in mobile ad hoc networks.

NLWCA organizes ad hoc networks in one-hop clusters (Figure 1) by using only information available locally. Each device elects exactly one device as its clusterhead, i.e. the neighbor with the highest weight. So far, a topological chain can be formed by so called sub-head nodes. A sub-head is a node that elects a neighbor node as clusterhead but at the same time is elected as clusterhead by some other onehop neighbor nodes. However, sub-heads can lead to more than three hops between a source clusterhead and its nearby clusterheads which could lead to a complex communication protocol. To obtain strict one-hop clusters, thus simplifying the protocol, a rule was added to the original NLWCA algorithm: a node that already elected a foreign node as clusterhead is not eligible to be elected by another node as clusterhead. From a graph theory point of view, one-hop clusters form a dominating set.

The main goal of NLWCA is to avoid superfluous reorganization of the clusters, particularly when clusters cross each other. To achieve this, NLWCA assigns weights to the links between the own node and the network neighbor nodes. This weight is used to keep track of the connection stability to the one-hop network neighbors. When a link weight reaches a given stability threshold it is considered stable and the device is called stable neighbor device. The clusterhead is elected only from the set of stable neighbors which avoids the reorganization of the topology when two clusters are crossing for a short period of time.

WCPD, on top of NLWCA, discovers nearby stableconnected clusters in a pro-active fashion. For the nearby clusterheads discovery algorithm, WCPD uses the beacon to detect devices in communication range. NLWCA and WCPD combined provide to each node, through the beacon (Figure 2), following information about each stable one-hop neighbor: its weight, its clusterhead ID, the ID set of discovered The WCPD broadcasting algorithm is simple: the broadcast source node sends the message to the clusterhead, which stores the ID of the message and broadcasts it to the onehop neighborhood. After that, it sends it to all nearby clusterheads by multi-hop unicast. The inter-cluster destination nodes repeat the procedure except that the message source clusters are omitted from further forwarding. Additionally, the information about the ID of the broadcast messages and their sources is stored for a given period of time to avoid superfluous re-sending of the message.

B. OLSR

The Optimized Link State Routing Protocol (OLSR) is a well known routing protocol designed for ad hoc networks. It is a proactive protocol; hence it periodically exchanges topology information other nodes of the network. Onehop neighborhood and two-hop neighborhood are discovered using Hello Messages (similar to beacons). The multipoint relay (MPR) nodes are calculated by selecting the smallest one-hop neighborhood set needed to reach every two-hop neighbor node. The topology control information is only forwarded by the nodes which are selected as MPR. Every node possesses then a routing table containing the shortest path to every node of the network. OLSR enables optimized flooding of the network by building a tree-like topology for every node from a source (Figure 3). Therefore, MPR selection constructs an optimal connected dominating set [START_REF] Adjih | Computing connected dominated sets with multipoint relays[END_REF].

C. Proposed Approach: Inter-Cluster communication

NLWCA and WCPD provide a stable-connected cluster architecture, however the broadcasting algorithm of WCPD needs many improvement on reachability performances compared to OLSR which performs very well on reachability but lacks in scalability and uses a lot of bandwidth [START_REF] Leclerc | Olsr and wcpd as basis for service discovery in manets[END_REF]. Our SLSF (Stable Linked Structure Flooding) protocol replaces the inefficient broadcasting mechanism of WCPD with the ICR mechanism. SLSF combines the advantages of all the protocols NLWCA, WCPD and OLSR: scalability, stability, reachability, while keeping the drawbacks low (i.e. the bandwidth usage). SLSF forms a first level of hierarchy with a dominating set using NLWCA. Considering the dominant nodes of the underlying level (NLWCA), it forms an optimal connected dominating set with the ICR mechanism. The first level reduces the network to its dominant nodes and the second level insures shortest-path connectivity and minimal relay nodes among dominant nodes of the first level.

The remainder of this paper is organized as follows. Section II presents the related work. Sections III and IV describe the building blocks of our SLSF protocol. In Section V we evaluate our approach and present the simulation results. We conclude our paper and present the future work in section VI.

II. RELATED WORK

In ad hoc networks forwarding strategies should be employed to avoid broadcast storms (i.e. a message forwarded by all the nodes in the network). As depicted in [START_REF] Ni | The broadcast storm problem in a mobile ad hoc network[END_REF], broadcast storms can be counter-measured using several schemes i.e. probabilistic, counter-based, distance-based, location-based and cluster-based. We use the latter scheme, cluster-based, since it is the only one based on network topology information. The cluster architecture used in this paper solely relies on locally available information. [START_REF] Ni | The broadcast storm problem in a mobile ad hoc network[END_REF] proposes a clustering technique where the clusterhead is elected after a message exchange among the neighbors. The route is constructed on demand. The broadcasted request is sent among all the clusters being forwarded by all the nodes reaching nearby clusters (also called gateway nodes). Every crossed cluster adds its address to the message, so that after reaching the destination node, the response is sent back through the path collected during the request phase.

In [START_REF] Haas | The zone routing protocol (zrp) for ad hoc networks[END_REF] the authors present the Zone Routing Protocol (ZRP), which is a hybrid routing approach. It combines proactive routing inside a zone using bordercast and on-demand routing outside. A node, sending a message, first checks if the destination is inside its zone, if not it bordercasts the message to its gateway nodes. Those nodes repeat the same process until the message reaches its destination. As every zone is centered on the current node, in terms of dissemination, the ZRP results in plain bordercasting a message ahead its destination. To route inside a zone, ZRP needs k-hop information (k>1), which results in scalability issues similar to OLSR.

In [START_REF] Foroozan | A high performance cluster-based broadcasting algorithm for wireless ad hoc networks based on a novel gateway selection approach[END_REF], the authors construct elected clusters based on beacon information which provides the number of neighbors and their stability represented by the number of beacons received since the node became clusterhead. This clustering approach is similar to NLWCA but does not rely on both a link weight and a node weight. Thus NLWCA has more flexibility in terms of cluster selection. Another similarity to our approach is the forwarding node selection (named gateway selection). They present a protocol which enables nodes to be selected as gateways to be able to relay outgoing traffic from a cluster. The main differences is that our approach requires no message exchange except the payload broadcast itself. The comparison with the protocol presented in [START_REF] Foroozan | A high performance cluster-based broadcasting algorithm for wireless ad hoc networks based on a novel gateway selection approach[END_REF] will not be further analyzed since to our knowledge the gateway selection process is insufficiently described.

Many ad hoc protocols use the selection of forwarding nodes to reduce redundant messages. In [START_REF] Peng | Ahbp: An efficient broadcast protocol for mobile ad hoc networks[END_REF] and [START_REF] Lou | Double-covered broadcast (dcb): a simple reliable broadcast algorithm in manets[END_REF] broadcast relay gateways are selected with 2-hop knowledge. However we chose the already described OLSR protocol because of its popularity in ad hoc networks and Mesh networks, and also because it is the only one proved to optimize coverage of 2-hop nodes through MPRs [START_REF] Jacquet | Performance analysis of OLSR multipoint relay flooding in two ad hoc wireless network models[END_REF]. MPRs build optimal connected dominating sets [START_REF] Adjih | Computing connected dominated sets with multipoint relays[END_REF]. We take advantage of this property to build paths among the (not-connected) dominating sets built by the NLWCA clustering algorithm. Many other approaches that construct distributed connected dominating sets exist [START_REF] Wu | On calculating poweraware connected dominating sets for efficient routing in ad hoc wireless networks[END_REF], [START_REF] Blum | Connected dominating set in sensor networks and manets[END_REF]. However, our goal is not to create connected dominating sets, but to disseminate information over the stable structure built by NLWCA, using ICRs between the clusterheads to optimize nearby-cluster-paths. Thus, as a result of our structure, we have connected dominating sets, but only between a clusterhead and its nearby clusterheads.

III. SLSF -INTER-CLUSTER RELAY

Multiple paths to reach a given clusterhead requires choosing one path prior to another. We use a next-hop selection inspired by the MultiPoint Relay (MPR) mechanism of OLSR to select the forwarding neighbors. MPR mechanism ensures full optimized coverage of 2-hop nodes. In our cluster architecture the set of discovered nearby clusters represent the "2-hop nodes" of OLSR.

We name Inter-Cluster-Relays (ICR) the nodes selected as next-hop. The goal of ICR is to reach all nearby clusterheads with the minimal set of 1-hop neighbors while optimizing the hop-count. The ICR nodes are calculated by selecting the smallest one-hop neighborhood set (directly connected nodes) needed to reach every nearby clusterhead. ICR selection remains simple and straightforward because the possible inter-cluster configurations are restricted by the underlying one-hop cluster topology (see Figure 4 for examples). SLSF, on top of NLWCA, discovers the nearby clusters (similar to WCPD) by reading the neighbor beacons. The improvement and novelty relies on the ICR selection which avoids superfluous network communication overhead without any additional message exchange. SLSF keeps the last beacon of every one-hop neighbor in cache. Hence every node has the following information locally available about each stable 1-hop neighbor: its weigth, its clusterhead ID, the ID set of discovered clusterheads and their respective path length.

ICR selection occurs as follows: i. Select as ICR the neighbors that are the only ones reaching a particular nearby clusterhead (CH). ii. Remove the now covered clusters from the list. iii. Remove for every neighbor from the announced CHlist the entries with a worse (greater) hop count than the best one (i.e. keep only shortest path entries for ICR selection).

1. Calculate the cluster reachability for every onehop neighbor (i.e. number of foreign CHs the neighbor announces in its beacon). 2. Select the neighbor with the best reachability.

• Remove the now covered clusters from the list. 3. Else if equivalent: select the neighbor with the shortest path to the remaining to be covered CHs.

• Remove the now covered clusters from the list. 4. Else if equivalent: select the neighbor with the highest weight.

• Remove the now covered clusters from the list. 5. Else if equivalent: select the node with the biggest IP address.

• Remove the now covered clusters from the list. 6. While there is a not-covered CH, go back to 1.

A. 3-hop Inter-cluster case

NLWCA builds one-hop clusters, thus it permits up to three hops (two slave nodes) between clusterheads. ICR selection with two hops (one slave node) between clusterheads (Figure 4b) is straight selection by the clusterhead, however an additional hop (Figure 4d) requires additional attention.

A further hop involves an additional forward of the message to reach the nearby clusterhead. For example on Figure 4d, the source CH2 designates a node as ICR (here the blue slave neighbor of CH2). The designated ICR has to make a choice between one of the two (orange) slaves of CH1. This choice is simply computed by using ICR selection. Additionaly, the CH list aimed by the ICR selection is smaller since only 2 hops away CHs are to be considered, opposed to a list of all nearby CHs for the source CH's ICR selection (up to 3 hops away).

B. When to select ICR nodes?

ICR selection is done based on events. Every time a change in the stable neighborhood that influences the ICR calculation occurs, the ICR selection is re-calculated. Thus broadcasting or forwarding a message using ICRs is immediate: replace the ICR set in the message with the one locally pre-calculated.

Further detail on how broadcasts in network occur are given in section IV.

C. ICR: The big picture

To highlight the gain of ICR selection, Figure 6 shows an example with 5 clusters where the message source clusterhead S sends a broadcast. The broadcast of S will have the format shown on figure 5. On reception of this broadcast only nodes 1, 2 and 3 will forward the message, while the other neighboring nodes process the message silently. Note that node 1 selects 4 and 5 as ICR according to section III-A "3-hop Inter-cluster case".

We see that ICR selection reduces a lot the number of forwarding nodes. As an example, on Figure 6 there are 23 nodes in the network and only 10 nodes (including the clusterheads) are emitting to reach all the nodes in the network. Every clusterhead will emit the message once in order for their slave to receive it and if necessary include their local ICR selection for further forwarding in the network (see section IV). In comparison, there would be 15 nodes forwarding the message using OLSR. This is due to OLSR using only 2-hop information while SLSF uses 1-hop cluster information which represent information from up to 3-hops away. While 3-hop knowledge usually increases the amount of information to collect using clusters reduces drastically the nodes needed to keep track of for ICR selection.

IV. SLSF -BROADCAST

At this point, every communication occurs between one cluster and its nearby clusters. To enable communication with foreign clusters, we propose a simple broadcast mechanism. A more sophisticated foreign cluster-broadcast mechanism would be out of scope for this paper.

Our broadcast mechanism is simple now that we only need to deal at cluster level. A node willing to broadcast a message through the network will, unless it is its own CH, send it to its CH. The original message contains the source address, a corresponding sequence number and of course the payload data. The CH puts its own address as last crossed CH and adds a corresponding sequence number to the message. Finally it forwards it to all its nearby CHs using the ICR mechanism. On reception the nearby CHs replace the last crossed CH address with their own address and replace the corresponding sequence number. The ICR set is also updated, while excluding from the ICR selection the cluster the message came from.

To avoid superfluous re-sending of the message SLSF stores information about the ID of messages and their sources for a given period of time. As an example, on Figure 7 the node (N) 1 sends a message to its clusterhead (CH) 9. CH9 adds its own address, a corresponding sequence number (here 64) and the ICR set (here N3). While N2, N6 and N3 receive the message, only the latter will forward it. CH8 receives the message from N3 and replaces the last crossed address and sequence number with its own. The ICR set of CH8 is N7. N3 will not be selected again as the message came from CH9. N7, as an ICR, will forward the message but changes the ICR set to N2 since it is the "3-hop inter-cluster case" (section III-A). N2 then forwards the message to its CH11 which replaces the last crossed address and sequence number and leaves the ICR set empty. As a result, the message reaches all the nodes in the network.

V. SIMULATION & RESULTS

To evaluate the performances of our SLSF protocol, we implemented the three protocols (OLSR, NLWCA/WCPD and SLSF) on top of the JANE simulator [START_REF] Gorgen | Jane-the java ad hoc network development environment[END_REF] and performed several experiments.

For those experiments we used the Restricted Random Way Point mobility model [START_REF] Blažević | Self organized terminode routing[END_REF], whereby the devices move along defined streets on the map of Luxembourg City for 1000 seconds. For each device the speed was randomly varied between [0.5;1.5] units/s and the transmission range set to 25 units/s. At simulation startup, the devices are positioned at random selected crossroads and the movement to other crossroads is determined by the given random distribution seed. For each experiment 10 different random distribution seeds were used in order to feature results from different topologies and movement setups. For the used mobile environment where nodes move with low speeds between 1.8 and 5.4 km/h the NLWCA link-stability threshold is set on 2 [START_REF] Andronache | Nlwca node and link weighted clustering algorithm for backbone-assisted mobile ad hoc networks[END_REF]. Bandwidth used in order to build the topology for 100, 200 and 300 nodes Fig. 9. Overall number of sent messages and node receivers for 100 and 300 nodes. For visibility sake the results where smoothed with a polynomial equation of the 16th grade. Simulations were done to determine the bandwidth used by the protocols in order to build the topologies and the information dissemination performance on top of the different topologies. Then we compared the efficiency of the protocols and finally made a static evaluation to compare information dissemination solely on MPR and ICR performances.

OLSR exchanges the sets of one-hop neighbor nodes with every node in communication range. Similar to OLSR, SLSF uses beacons to exchange the list of the discovered nearbyclusterheads with the one-hop neighbor nodes. Since SLSF and WCPD have the same beacon structure, they use exactly the same bandwidth. To find out the network load produced during this phase, the size of the exchanged data sets were tracked every second of the simulation: for OLSR the size of the one-hop neighbor sets and for SLSF and WCPD the size of the discovered clusterheads.

In order to monitor the information dissemination performance (Reachability) a node was chosen to broadcast a message every 10 seconds during different simulation runs using different distribution seeds. The number of sent messages (i.e. broadcasts and unicasts) during the dissemination and the number of reached network nodes were tracked.

As shown in Figure 8 OLSR uses a higher bandwidth in both sparser (100 nodes) and denser networks (300 nodes). This was expected since OLSR is exchanging the set of onehop neighbors needed for the MPR nodes election, while SLSF only exchanges the set of locally discovered nearby clusterheads. The NLWCA protocol elects one clusterhead in each one-hop neighborhood, hence the number of clusterheads is a fractional amount of the total number of nodes.

SLSF uses approximatively the same amount of forwarders than WCPD but reaches from 5% (low density) to 20% (high density) more nodes (Figure 9). Note: the periodic behavior is due to the smoothing for visibility reasons. This is the pure gain of ICR selection which optimizes the number of forwarding nodes. The tracking results regarding the reachability show that for sparser networks OLSR and SLSF perform equivalently. When stabilized, for denser networks, OLSR reaches about 85% (40% for SLSF) of all nodes at the cost of 65% (5% for SLSF) of forwarders. SLSF reaches half the number of nodes compared to OLSR, but with 10 to 12 times less forwarders. Furthermore, while forwarders increase with OLSR proportionally to the number of receivers, their number is stable and low with SLSF. OLSR is optimal as far as reachability is concerned. When bandwidth and computation (active nodes) is concerned, OLSR may not be adequate.

Subsequently we calculated a "quality-cost" ratio, percentage of nodes reached divided by the bandwidth used, extracted from the results of Figures 8 and9. We see in Figure 10 that SLSF is in average three times more efficient than OLSR.

SLSF relies on stable structures built by NLWCA: only nodes considered as stable will receive the message. So finally, to compare the performances on equal levels, we experimented OLSR and SLSF in a static scenario where all nodes are considered stable connected. The experiments where done on a 300x300 units surface with 100 to 300 nodes randomly positioned using 100 different topology seeds. Again, the number of forwarding and receiving nodes using MPR and ICR selection where tracked. The results on Figure 11 show that SLSF outperforms OLSR in terms of ratio of receiving nodes over forwarding nodes. With increasing density on average with OLSR about 85% of the receivers are also forwarders, whereas in SLSF this amount decreases from 60% towards 30%.

As a further step, we started experimenting OLSR and SLSF, that we implemented in JAVA, on real world devices (Nokia N800). One of our research goal is to create stable context-aware service discovery. Therefore we used SLSF as basis for the service discovery protocol Zeroconf (mDNS/DNS-SD) [18]. Multicast DNS (mDNS) uses the mutlicast tree for DNS information dissemination. We replaced the mutlicast tree by our SLSF structure. The first results are positive and we plan as future work to further experiment and evaluate advanced usage scenarios with real devices.

VI. CONCLUSION & FUTURE WORK

In this paper we propose SLSF, a flooding protocol which selects ICRs (Inter-Cluster Relays) to optimize the communication among the stable-connected cluster architecture.

The goal of the ICR selection is to reach all nearby clusterheads with the minimal set of 1-hop neighbors while optimizing the hop-count. Further on, it reduces message collisions and hidden terminal problems since only a few nodes in the same neighborhood emit at the same time.

In our use case, ICR selection on top of the stable-cluster architecture reduces substantially the number of forwarding nodes compared to OLSR. Overall, SLSF is very efficient and performs well in high density networks while keeping the used bandwidth very low. We are concerned that this efficiency can have a drawback: one single lost message is much more penalizing than with OLSR. Therefore for the short term, we plan to add and study the impact of fault recovery mechanisms in order to improve reachability.

As future work we also plan to evaluate the performances using other mobility models and topology settings and also assess the results on real devices with and advanced usage scenario. We also are investigating the usage of SLSF as basis for the service discovery protocol Zeroconf with simulation and real world experiments.

Fig. 1 .

 1 Fig. 1. Example of two clusters built by NLWCA.

Fig. 2 .

 2 Fig. 2. NLWCA+WCPD Beacon.

Fig. 3 .

 3 Fig.3. OLSR topology for one source node in particular.

Fig. 4 .

 4 Fig. 4. Inter-cluster configuration examples where 1 and 2 are clusterheads.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Format of a broadcast message with payload.

Fig. 7 .

 7 Fig. 7. Foreign-cluster broadcast -Format and Path of a message sent from node 1

Fig. 8 .

 8 Fig. 8.Bandwidth used in order to build the topology for 100, 200 and 300 nodes

Fig. 10 .

 10 Fig. 10. Efficiency of Bandwidth usage for 100 and 300 nodes Fig. 11. Static scenario with 100 to 300 nodes

We plan to evaluate performances in the context of the French National Research Project SARAH, deploying a large scale ad hoc network inside a museum.

VII. ACKNOWLEDGEMENTS

This work was supported in part by the Conseil Régional de Lorraine and the ANR (SARAH Project).


Jean-Charles Quinton
email: quintonj@loria.fr

Exploring and Optimizing Dynamic Neural Fields Parameters Using Genetic Algorithms

The Continuous Neural Field Theory introduces biologically-inspired competition mechanisms in computational models of perception and action. This paper deals with the use of Genetic Algorithms to optimize its parameters, as to guarantee the emergence of robust cognitive properties. Such properties include the tracking of initially salient stimuli despite strong noise and distracters. Interactions between the parameter values, input dynamics and accuracy of model, as well as their implications for Genetic Algorithms are discussed. The fitness function and set of scenarios used to evaluate the parameters through simulation must be carefully chosen. Experimental results reflect an ineluctable tradeoff between generality and performance.

I. INTRODUCTION

A. Continuous Neural Field Theory (CNFT)

The cortex can be roughly approximated as a two dimensional structure made of myriads of elementary units called cortical columns [START_REF] Burnod | An adaptive neural network: The cerebral cortex[END_REF]. Despite this apparent homogeneity, the organization of the body and innate structure of the brain constrain the afferent connections in different parts of the cortex. Self-organization and competition within the cortex lead to differentiated areas dedicated to specific modalities, communicating through long-range connections and associative areas. But this process furthermore introduces a local topological organization within the cortical sheet, where close units tend to respond to similar input configurations [START_REF] Malsburg | Self-organization of orientation-sensitive cells in the striate cortex[END_REF]. This local correlation reinforces local excitatory connections, while balanced by global inhibition. Whether the connectivity patterns and input driven organization are partly innate or engaged in a form of circular causality, this phenomenon leads to the emergence of important perceptual and behavioral properties.

Extensive research has been done on how lateral interactions might evolve through the course of development [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF], [START_REF] Katz | Development of local circuits in mammalian visual cortex[END_REF], the Continuous Neural Field Theory (CNFT) also examines the dynamic properties of such neural maps [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF]. Properties include the emergence of bubbles of activity able to robustly track input stimuli in spite of strong noise or numerous distracters [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. This kind of attentional property combined with the dynamic competition resulting from the large scale inhibition is one of the characteristics needed to adopt complex goal-oriented behaviors. Bubbles can indeed shape the attractor landscape and constrain the fast dynamics of sensorimotor maps over an extended period of time. Radical changes constantly occur on low-level sensory maps during the course of many common actions and would lead to chaotic behavior without such modulations (consider for instance the optical flow when temporarily turning away from a visual target).

Although general convergence results have been demonstrated and parameters constrained to allow the emergence of such bubbles of activity [START_REF] Taylor | Neural bubble dynamics in two dimensions: Foundations[END_REF], many parameters need to be tuned through simulation. Mathematical analysis indeed provides simple rules to determine functional parameters for which numerical simulations can be run [START_REF] Alexandre | On practical neural field parameters adjustment[END_REF], and experimentation then allows to further test the robustness of the model, but the search for generic and optimal parameter values is limited for complex input dynamics. Parameters tightly interact in the equations and a slight change of a single parameter sometimes lead to dramatic changes in the behavior of the model on specific simulation scenarios. With the assumption that correct sets of parameters should lead to the emergence of bubbles able to robustly track salient input stimuli, many values are simply non functional: they either lead to the saturation of the field, the lack of any coherent activity, or the self-maintenance of static bubbles independently of the input dynamics. In this paper, we therefore focus on the selection of optimal parameters using Genetic Algorithms.

B. Genetic Algorithms (GA)

There are several reasons for choosing Genetic Algorithms (GA) as a metaheuristic for optimizing the attentional properties of the CNFT. Since we are interested in the robustness of the model, randomness is introduced in the input as to increase its variability and limit the assumptions to the kind of noise, shape and trajectory of distracters used. Whereas no random term is introduced in the equations governing the model and that its inner dynamics is entirely deterministic, several subsequent simulations may lead to totally different results. With inadequate parameter values and noisy inputs, the initial phase leading to the potential emergence of a bubble is highly chaotic, the reinforcement of the bubble requiring spatiotemporally coherent stimulations. In this sense, GA allow to test parameter values again and again for each generation, not simply relying on a single successful trial for a single individual. Purely coincidental survival of a model configuration will lead to its elimination and that of its progeny with a high probability. From experimental results and hand-tuning of parameters, there may be several distant and local optima, wide plateaus and even ridges for some subregions of the parameter space. This makes simple methods like the steepest descent inadequate for the problem, since a convex function was not found to evaluate the performance of the model. Even if it must underlined that many local minima lead to acceptable results, the search for optimal parameters also make simulated annealing techniques not very efficient though perfectly adapted to large search spaces, as restarts are often needed. Once groups within the global population have converged to local minima, GA keep the adaptive pressure increasing while still allowing short and long jumps within the search space respectively through mutations and cross-overs. Finally, it is not obvious that fixed parameter values can genetically provide a good functional basis on which ontogenetically reinforced synaptic connections might selforganize. Using GA to evolve the parameters is therefore a way to test the biological plausibility of having fixed global inhibition and local excitation, instead of having them co-emerge with or simply adapt to other cortical and thalamic connections through learning and development.

II. METHOD

A. CNFT equations and parameters to optimize

Adopting the notations introduced by Amari [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF], we will present a version of the CNFT using a single-layer field equation of lateral inhibition type. The focus neural field on which bubbles may form is represented by a manifold M in bijection with [-0.5, 0.5] 2 and the membrane potential at the position vector x and time t on this field by u(x, t). Similar notations are used for the stimulation s(x, t), that is the only external input to the system. The dynamics of the membrane potential is described by Equation [START_REF] Alecu | An empirical evaluation framework for qualifying dynamic neural fields[END_REF].

τ ∂u(x, t) ∂t = -u(x, t) + x ′ ∈M w(x, x')u(x', t)dx' + s(x, t) + h ( 1
)
where h is the resting potential (taken equal to 0) and w(x, x') the lateral connection weight function satisfying the following equation:

w(x, x') = Ae -|x-x'| 2 a 2 -Be -|x-x'| 2 b 2 (2)
At first glance, Equation (1) may seem to lack degrees of freedom, as the relaxation term -u(x, t) and stimulation s(x, t) share the same coefficient. However, whereas introducing non-linearities in the lateral competition term (such as a sigmoidal function) or integration of receptive fields in the stimulation term might totally change the dynamics of the CNFT, introducing a simple multiplicative coefficient for the input will only scale the focus (as long as h is null). Indeed, if we replace s(x, t) by α * s(x, t) in Equation ( 1) and let u ′ = u/α, u ′ then satisfies the original equation (i.e. without α).

When discretizing this equation on a regular mesh of units for simulation purpose, a spatial resolution parameter n is additionally introduced. However, the dynamics of the model is theoretically independent of both the spatial and temporal discretization parameters, except when chosen too coarse (taking a low n or high dt). This was confirmed in experiments with n ranging from 10 to 100, with n = 30 being sufficient for the expected properties to emerge. In the following sections, let n = 50 units and dt = 0.1 seconds. These parameters will therefore not be optimized and should not be referred to anymore.

Remaining parameters of Equation ( 1) and ( 2) that need to be tuned are synthesized in Table I. Although most parameters are limited to positive values, mathematical analysis of the equations further constrains the parameters to allow the emergence and disappearance of bubbles of activity, avoiding at the same time the saturation of the field [START_REF] Taylor | Neural bubble dynamics in two dimensions: Foundations[END_REF]. A Mexican hat shape for the lateral weight function is thus adopted here, implying the conditions A > B and a < b to be satisfied. As to introduce such constraints in the evolutionary process, we let K = B/A and k = a/b with (K, k) ∈ [0, 1]×]0, 1]. This does not reduce the number of parameters to be optimized but limits the dramatic effect that slight changes in parameters might have on the dynamics (for instance if a mutation leads from A > B to A + dA < B).

B. GA implementation

GA are applied to evolve a population of models, each individual genotype being encoded by a parameter vector

v i = (p k i ) k∈[1.
.s] , s being the number of parameters. In our case, this vector must account for all CNFT parameters and is thus defined by

v i = (A i , k i , K i , b i , τ i ).
The phenotype translation corresponds to the straightforward generation of the lateral connectivity pattern and the use of differential equation parameter τ i in Equation [START_REF] Alecu | An empirical evaluation framework for qualifying dynamic neural fields[END_REF].

A new population g + 1 is then composed of the individuals that were the most successful at generation g (see the definition of the fitness below) as well as new individuals obtained through the cross-over and mutation of two individuals taken randomly from population g. This guarantees that the best results ever obtained during the artificial evolutionary process will never be discarded from the population. The size of the population remains the same over generations as to maintain a constant evaluation cost. The cross-over operation produces a new vector where each parameter is assigned to a random weighted sum of the parent values. Mathematically, a genotype

v c = (γ k p k 1 + (1 -γ k )p k 2 )
is created where v 1 and v 2 would be the genotypes of the parents and (γ k ) random coefficients taken in (0; 1) s . Mutations are random variations ∆p k c applied to random components of the child vector (probability

p mut ). Let ∆p k c = π k |p k max -p k min | with π k randomly taken in [-m, m]
and m a factor introduced to constrain the mutation maximum amplitude. However, the algorithm also ensures that the associated parameter value remains within its imposed limits [p k min , p k max ], as K ∈ [0, 1] for instance. Although dependent on the chosen mutation range, these operations allow fast jumps over the entire parameter space, expansion through mutations and contraction through cross-overs.

Whereas GA are quite generic, the difficulty generally lies in the choice of the right genes and correct fitness function. In our case, the genes selection and encoding is intuitive and straightforward, but the fitness function goes through the evaluation of complex dynamical simulations. The overall process and its exact formulation must therefore be specified with great care. The following two sections are dedicated to the underlying assumptions and resulting consequences of the choice of fitness and associated scenarios used for its evaluation.

C. Fitness thoughts

Since the focus is on the distributed and emergent properties of the CNFT, a static evaluation of the parameters is not sufficient and simulations of the model must be run. The problem is now the high dependence of performance on the chosen scenario. If the goal is to optimize the entire spectrum of emergent properties displayed by the CNFT models, all these properties must be reflected in the selected scenarios.

As the model might easily bifurcate in the first few seconds of simulation based on initial conditions and amount of noise introduced, the fitness function is computed on a sequence of n stereotyped scenarios. This allows to easily account for the various properties of the model, but also for the intersimulation variability. The global fitness f (v) for a given parameter vector v is defined by:

f (v) = 1 n n k=1 f (s k , v)
where f (s k , v) corresponds to the evaluation of v for a single scenario s k . This subfunction must be defined as to synthesize the accuracy and efficiency of the model over the entire simulation in a unique scalar value, defined as follows.

f (s k , v) = error(s k , v) × conv(s k , v) × shape(s k , v) (3)
In Equation ( 3), error corresponds to the mean normalized distance between the center of the focus bubble and the center of the input stimulus the model is supposed to track and focus on. An instantaneous error (e t ) is thus computed for each iteration in the simulation, and then averaged on the last seconds of simulation to rule out the initial chaotic dynamics. This computation already assumes that the quality of the model is correlated to its ability to make a single stereotyped bubble emerge. For the details of the approximation of the center of the focus bubble, please refer to [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. This error measure also supposes that we a priori know or can easily determine which input stimulus is attended, and that so-called input stimuli merely exist. To make all these points clear, we determine the input flow independently of the model dynamics, as designers and evaluators. The model has indeed no mean of changing the input dynamics by acting on it (as would a typical sensorimotor system engaged in closed loop interactions with its environment). We therefore define what is considered to be coherent stimuli for the model and optimize the parameters in relation to this definition. As observers of the model dynamics, we can also compute 1. Simulation steps for one iteration. 1 The input map is generated from the current scenario parameters. 2 The model is updated, only using the input map and CNFT parameters. 3 Statistics are computed based on the known input positions and focus map produced by the model. statistics on the model and determine the closest stimulus to the focus bubble, information the model has of course no access to. Figure 1 proposes an illustration of these considerations.

Again in Equation ( 3), conv is defined as the convergence time, i.e. the time needed for e t to drop below a given relative threshold. In other words, let conv be the minimal time such that whatever t c ∈ [conv, t end ], e tc < α × min t (e t ) + (1α) × max t (e t ), where α is the threshold parameter set to 0.2 (see Figure 2). If the condition cannot be satisfied, conv is set to the ending time of the simulation t end . An additional variable shape is also taken into consideration to constrain the profile of the focus bubble. A lack of inhibition on the focus field or a limited coherence of the activity under the bubble may lead to a focus field with a wide support (in the sense of measure theory) and non stereotyped bubbles. Although the computation of the bubble center is not altered by such distortions, spread out activations may be inadequate for decision or motor control. At a given instant t of simulation, the instantaneous shape computation satisfies:

shape t = x∈M |u * (x, t) -u(x, t)| dx
To put it briefly, this equation computes a similarity between the actual focus u and an hypothetical ideal focus field u * where a single stereotyped bubble would emerge in the presence of the tracked stimulus, defined by:

u * (x, t) = I t A -B w + (c t , x) 2
where w + is the positive part of the lateral weight function w defined in Equation ( 2), c t the bubble center at time t and I t the tracked stimulus amplitude (used as a normalization factor for the CNFT profile whose maximum value is A-B).

D. Scenarios

As stated in the previous section, several complementary scenarios are needed to account for the various emergent properties of the CNFT (already described in [START_REF] Amari | Dynamics of pattern formation in lateral-inhibition type neural fields[END_REF], [START_REF] Rougier | Emergence of attention within a neural population[END_REF] and [START_REF] Chevallier | Visual focus with spiking neurons[END_REF]). For instance, the a value is highly correlated to the distance between the stimuli to focus on: a large value may improve the results only when unambiguous distant stimuli are presented. If they are too close, the CNFT will merge their activity under a single bubble that will not track either of the stimuli. But the model should not be dependent on the input dynamics within a given functional range, even if a multi-scale, multi-resolution approach should be adopted for a wider range of conditions. The parameters are thus optimized as to produce the best results for all three following scenarios, although non optimal for each scenario taken separately: C 2 bell-shaped distant stimuli s 1 and s 2 are introduced at time t = 0. Their intensity are refreshed every 5 seconds and governed by I 1 = 0.9 and I 2 (t) = 0.5+0.5 cos(π × (t/5)). D 1 bell-shaped stimulus of standard deviation 0.1 and intensity 1.0 follow a circular trajectory of radius 0.2 around the point (0, 0) at 10 deg/s from t = 0. From t = 1, 5 distracters are added and take new random positions on the field every 1 s. E 1 bell-shaped stimulus (same as above). At t = 0, Gaussian noise with a standard deviation of 0.5 is added.

Scenario C instantiates the competition property of the CNFT. The model should initially focus on the most salient target s 2 before switching when its intensity will drop to 0 for t = 5. The shifts in attention introduce an hysteresis phenomenon in the dynamics, as the resilience to noise and distracters partially comes from the combined effect of the feedback and input. For t > 5, the focus should remain on s 1 as s 2 will never reach a significantly higher intensity (at most 0.1 higher). This scenario also tests the ability of the model to relax when the stimulation disappears, i.e. not to sustain saturated bubbles independently of the input. Scenario D tests the capability of the model to track a single stimulus despite distracters of the same intensity. As to guarantee their spatiotemporal coherence, distracters are only randomly moved every second, allowing them to destabilize the focus for several iterations. Since the model cannot distinguish between the stimulus and the distracters, these are introduced with a delay of 1 second after the beginning of the experiment, as to allow the convergence of the focus field on the desired input. Notice that movement is introduced to make the experiment valid. Indeed, staticity is the default in Equation ( 1), as the feedback at x and time t + dt comes from the same field position x at time t. The speed of the stimulus is thus relatively high in order to force the model dynamics to be reactive.

Finally, scenario E checks if the model is robust to random noise on the entire field. The bubble is therefore highly distorted as shown on Figure 3. The noise is added from the start, as to ensure the selection of parameters allowing the convergence of the CNFT model under such conditions, and not simply its resilience to noise. When interacting with real environments, noise is constantly present and must be taken into account at any time. It would therefore be highly artificial to take its absence for granted during the first iterations. In order to show the qualitative differences in results obtained when removing this constraint, let E' be a new scenario identical to E except for the introduction of noise at t = 1 instead of t = 0. Results using different combinations of scenarios are detailed and discussed in the following section.

III. RESULTS

The GA described in the previous section has been applied to a population of 20 to 100 individuals. Their number did not quantitatively change the results as long as enough generations were produced to allow the convergence of the genome (>20 generations) and that mutation had a sufficient effect on the parameters (m = 0.1, p mut = 0.1). The parameters were randomly chosen with a uniform distribution in the following ranges:

A ∈ [0.1, 2], K ∈ [0.1, 1], b ∈ [0.01, 2], k ∈ [0.1, 1], τ ∈ [0.1, 2]
. Some of these boundaries were selected to avoid divisions by zero, the others to get rid of trivial or absurd dynamics. For each new generation, the 40% top parents were kept unchanged (the ones with the highest fitness), while 60% of the individuals were randomly selected as parents and randomly combined to produce the children.

The discrete implementation and equations used for the computer simulations were inspired by Rougier and Vitay [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. A singular value decomposition (SVD) of the lateral connection kernel leads to the decomposition of the 2D convolution into two orthogonal 1D convolutions. The resulting drop in complexity allows a ten fold speed-up on 50 × 50 maps. On an old Intel Pentium M 735 based laptop configuration, each update of the CNFT model takes approximately 37ms, each simulation therefore lasting for 3.7 seconds on average with dt = 0.1 and a time limit of 10 seconds. The optimization of the parameters in 20 generations using the GA on 20 individuals then takes about 74 minutes for one set of 3 scenarios and one fitness definition.

As stated previously, GA run on different sets of scenarios lead to different optimal parameters. To also show the importance of the expression used for the fitness function, results are given whether the shape variable is used. Optimal sets of parameters after convergence of the GA are then used to test the associated model on all four scenarios C , D , E and E' . The results are synthesized in Table II and discussed below.

A. Discussion

When looking at Table II, one can observe the large variations in optimal parameters and performance on the different scenarios (several orders of magnitude). The parameters obtained when using neither E nor the shape variable are the most efficient in terms of reactivity, convergence speed and mean error, except for E where the model never converges (mean error = 0.27). These measures are however correlated as a low value for τ and high values for A and B lead to large variations of the potential when no bubble is present to inhibit the rest of the field. This initial chaotic behavior, not suitable when noise is introduced at the beginning of the simulation, is illustrated on Figure 4. Thus, although the GA produces a lower error on several scenarios when compared to the hand tuned parameters proposed in [START_REF] Rougier | Emergence of attention within a neural population[END_REF], the optimal solution does not reflect all the constraints an expert might take into account.

The shape variable introduces such an additional constraint (bubble stereotyping), whose effect can be seen on Figure 5. The increased spatial coherence of the bubble requires an excitatory core on the lateral connection kernel, therefore resulting in changes in all its parameters. This leads to a more "Mexican hat like" lateral weight function compared to the mainly inhibitory function of the previous set of GA scenarios. The dynamics is then smoother and the system is still quite reactive to input changes, but still not able to cope with noise introduced from start (mean error = 0.26).

When using scenario E in the GA, τ jumps from 0.29/0.27 to 0.45 and is accompanied by a proportional change of A and B to compensate the slower temporal integration. The results well display the tradeoff between robustness and reactivity, as to follow rapidly moving inputs, the model needs to take into account any stimulations, including noise and distracters. Setting τ thus constrains the ability to switch between targets and the amount of spatiotemporally coherent stimulations required to destabilize the bubble.

At this point, a last parameters optimization is performed to check if the shape variable and scenario E would not be sufficient to account for all the CNFT emergent properties. Although still functional in most cases, the GA converges to a different optimal solution, where τ takes a high value of 0.64. The model is of course optimal on scenarios E and E' , but is slow to switch between targets or simply to converge on one. Whatever the set of scenarios chosen, the results are however quite robust, as correlated variations up to 20% in the parameters values relatively to the optimal result only lead to mean error fluctuations lower than 0.03 (the results are obtained using the best parameter sets found for the last generations). This has to be contrasted with random changes of the same amplitude which can produce non converging simulations.

As a concluding remark for this section, whatever values the parameters may take in Equation (1), the error cannot converge to zero with moving stimuli. This is true except if there is only one symmetric stimulus, no noise and no distracter, in which case the degenerated model with τ =dt, A=B=0 will no longer hold any CNFT property as u(x, t) = s(x, t). This explains the efficiency of the model on scenario C , where the staticity of the competitive term in the equation (the lateral connectivity kernel for position x being centered on x) is reflected by the staticity of the input. As a consequence, the model is unable to track a stimulus passing over a static distracter, as staticity is the default. Whether we include such a scenario in the evaluation of the fitness function or simply test the optimal solution with it, no set of parameter values can succeed as the lack of anticipation is inherent to the equation.

IV. CONCLUSIONS

The optimization of the Continuous Neural Field Theory parameters using Genetic Algorithms produces improved results compared to hand tuned values. Even with expertise and intuition, the effect of parameter changes are hard to predict as the system involves a large number of interconnected units, feedback and non linear dynamics. All these ingredients of deterministic chaos indeed generate complex dynamics prior to the emergence of a bubble. In addition to the robustness of the model emergent properties, this optimization process improves the understanding of the interactions between the parameters. The analysis of the results and inclusion of new scenarios also points to the limits of the considered equation and tradeoffs between the reactivity, robustness and generality of the model.

This not only opens the possibility to change the equation when limits are encountered, but the automatic exploration of the parameter space also allows a fair comparison of different implementations, as the number of parameters and their optimal values might greatly differ. Such extensions to classical Dynamic Neural Field equations include the introduction of delayed inhibition [START_REF]A dynamic neural field mechanism for self-organization[END_REF], the use of spiking units rather than rate coding [START_REF] Chevallier | Visual focus with spiking neurons[END_REF] or the introduction of a predictive term (to address the limitation presented at the end of the discussion). Other differential equations governing the neural fields have been proposed [START_REF] Folias | Breathers in two-dimensional neural media[END_REF] and the same metaheuristic applies. The method may also be useful to compare different implementations of the same equations, for instance approximating the maps with sparse representations for increased performance.

Experimental simulations are required as no mathematical framework is yet able to formally characterize the dynamics of an arbitrary DNF equation on arbitrary scenarios. Such an exploration is thus a way to determine if the behavioral differences observed during simulations come from side effects, badly tuned parameters or fundamental functional differences. Others have already acknowledged the need to optimize the parameters [START_REF] Alecu | An empirical evaluation framework for qualifying dynamic neural fields[END_REF] and in this perspective, we presented a low cost alternative to brute force approaches. For these later ones, coarsely sampling the 5-dimensional space explored in this paper ensures completeness, but already requires many times more simulations than Genetic Algorithms.

τ

  Inertial parameter. How fast can the activity change, whatever the dynamics of the stimulation. A Excitatory amplitude. How much close units should form coherent patches of activity. a Excitatory range. This constrains the size of the emerging bubbles. B Inhibitory amplitude. How much distant units should inhibit each others (relatively to A). b Inhibitory range. If b ≥ 1 (size of the neural field), we can consider the inhibition to be global.

  Fig.1. Simulation steps for one iteration.1 The input map is generated from the current scenario parameters.2 The model is updated, only using the input map and CNFT parameters.3 Statistics are computed based on the known input positions and focus map produced by the model.

10 LFig. 2 .Fig. 3 .

 1023 Fig.2. Typical profile of the error when stimuli are competing for the focus. After a few iterations where noise dominates and the focus center is close to the center of the map, the non-linearity introduced by the lateral competition induces a rapid convergence to one of the input stimulus. The instantaneous error et drastically decreases and definitely fall below the threshold used for the convergence time conv. Once the initial chaotic behavior has settled, the mean error is computed on the last seconds of simulation.

Fig. 4 .Fig. 5 .

 45 Fig.4. Comparison of the lateral connection kernels and error dynamics on scenario C . This scenario is the only one without a random component in it, allowing reproducible deterministic dynamics without using statistics over a large number of trials. Each row represents the results on a different set of scenarios and variables used to evolve the parameter vector population. When using the shape variable and scenario E , the model looses in reactivity but gains in smoothness and stability.

TABLE I DESCRIPTION

 I OF THE CNFT PARAMETERS USED TO OPTIMIZE THE MODEL BEHAVIOR.

TABLE II OPTIMAL

 II PARAMETERS AND SIMULATION RESULTS OBTAINED WITH DIFFERENT SETS OF FITNESS VARIABLES AND SCENARIOS (CELLS IN GRAY CORRESPOND TO THE SIMULATIONS WITH QUALITATIVELY LOWER PERFORMANCE THAN THE OTHERS)

	GA scenarios			Parameters				error for test scenarios
		A	a	B	b	τ	C	D	E	E'
	C, D, E ′ , ¬shape 0.135 0.34 0.128 1.00 0.29 0.05×10 -4 0.019 0.270 0.009 C, D, E ′ , shape 0.069 0.14 0.041 1.25 0.27 0.12×10 -4 0.019 0.256 0.015
	C, D, E, shape	0.074 0.28 0.062 0.88 0.45 0.75×10 -4 0.040 0.016 0.017
	E, shape	0.065 0.38 0.061 0.93 0.64 2.08×10 -4 0.071 0.016 0.017