Adam Gaier
email: adam.gaier@h-brs.de

Alexander Asteroth
email: alexander.asteroth@h-brs.de

Jean-Baptiste Mouret
email: jean-baptiste.mouret@inria.fr

Data-Efficient

Keywords: Exploration, Optimization, and Modeling of Diverse Designs through Surrogate-Assisted Illumination MAP-Elites, Surrogate Modeling, Quality Diversity

published or not. The documents may come

INTRODUCTION

Computational techniques for design optimization are often thought of by their creators as a nal step in the design process. Imagining their techniques will be used to push the limits of performance, algorithm designers judge success by the ability to re ne a design to its most optimal form [START_REF] Hornby | Automated Antenna Design with Evolutionary Algorithms[END_REF].

If, however, the goal is truly to support designers, this sole emphasis on optimality may be misplaced. Autodesk [START_REF] Bradner | Parameters tell the design story: Ideation and abstraction in design optimization[END_REF] recently conducted an interview to better understand how professional designers, engineers, and architects use design optimization tools.

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro t or commercial advantage and that copies bear this notice and the full citation on the rst page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). GECCO '17, Berlin, Germany © 2017 Copyright held by the owner/author(s). 978-1-4503-4920-8/17/07. . . $15.00 DOI: http://dx.doi.org/10.1145/3071178.3071282 2) Maximize the acquisition function, which balances exploitation and exploration, in every region of the feature space, producing an acquisition map.

0) Sample design space 1) Construct model 2) Maximize acquisition function 3) Sample acquisition map 4) Maximize performance estimation

3) Draw next samples to test on the objective function from the acquisition map. Repeat steps 1-3. 4) Maximize tness as predicted by the resulting model, producing a prediction map.

They found that optimization was most commonly used not at the end of the design process, but the beginning. Rather than using optimization to solve design problems, they were more commonly used to explore them.

Generating a range of candidate solutions that represent di erent design alternatives allows designers to explore various design concepts, and examine the trade o s they represent. These generated designs provide insight into the assumptions and consequences inherent to the problem de nition and constraints. Once constraints and objectives are reconsidered and adjusted, new designs are then generated and the process repeated.

This generative cycle allows designers to explore and describe complex design spaces, with high-performing solutions acting as concrete way points. They can then iterate on the designs found through this collaborative human-computer exploration and, after consideration of intangibles such as aesthetics, a design is nalized.

Multi-objective optimization is perhaps the most commonly used tool to produce a variety of designs. When objectives are in con ict, each design in the Pareto front represents a trade-o between them [START_REF] Deb | Unveiling innovative design principles by means of multiple con icting objectives[END_REF]. However, during the explorative process interest for designers often lies not only in the maximization of objectives, but the e ect of di erent design features on performance.

To probe the search space for interesting designs and design principles, new algorithms created speci cally for design space exploration should be applied. One such algorithm, MAP-Elites [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] explicitly explores the relationship between user-de ned features and performance. Designers select a few features deemed interesting or important, such as weight or structural strength, and MAP-Elites produces high-performing solutions which span the possible variations of those features. This illumination process reveals the performance potential of the feature space.

While e ective at nding a variety of high-performing solutions, the number of evaluations required by MAP-Elites is immense. The illumination process which produced the repertoire of hexapod controllers in [START_REF] Cully | Robots that can adapt like animals[END_REF], for example, required twenty million evaluations. In applications such as structural optimization or uid dynamics, where a single evaluation can take hours, this is unrealistic.

In computationally expensive problems it is common to make use of surrogate models, approximate models of the objective function, that are based on previously evaluated solutions [START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF][START_REF]A comprehensive survey of tness approximation in evolutionary computation[END_REF][START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF]. These models are constructed through the sampling of solutions based on an acquisition function, which balances exploitation and exploration to improve accuracy in high tness regions. These computationally e cient models can be used in place of the objective function during optimization, greatly accelerating the process. Incorporating surrogate-assistance into the evaluation-heavy illumination process has the potential to make MAP-Elites e cient enough to be used even in computationally expensive problems.

We present the Surrogate-Assisted Illumination (SAIL) algorithm to improve the e ciency, and so expand the applicability, of MAP-Elites. The value of integrating surrogate models into illumination relies on reducing computational cost while maintaining MAP-Elites' original capabilities, resulting in an algorithm that:

• Divergent -Produces a diversity of solutions which vary across a user-de ned continuum; • Accurate -Predicts behavior of the objective function in high-performing regions; • Optimal -Produces high-performing solutions; • E cient -Performs under computational constraints.

In broad terms SAIL works as follows (Figure 1, previous page): a surrogate model is constructed based on a set of initial solutions and their measured performance. MAP-Elites is used to produce solutions that maximize the acquisition function in every region of feature space, producing an acquisition map. New samples are then drawn from the acquisition map and evaluated, and these additional observations are used to improve the model. This process is repeated to produce increasingly accurate models of the high tness regions of the feature space. Performance predictions of the model can then be used by MAP-Elites in place of the objective function to produce a prediction map of estimated optimal designs.

RELATED WORK 2.1 Quality Diversity and MAP-Elites

Quality diversity (QD) algorithms [START_REF] Pugh | Quality Diversity: A New Frontier for Evolutionary Computation[END_REF] use evolutionary methods to produce an archive of diverse, high quality solutions within a single run. Rather than seeking a single global optimum, QD algorithms discover as many types of solutions to a problem as possible, and produce a best possible example of each type. For this reason they are also referred to as illumination algorithms, as they illuminate the performance potential of di erent regions of the solution space.

Among the few illumination algorithms, novelty search with local competition (NSLC) [START_REF] Lehman | Evolving a diversity of virtual creatures through novelty search and local competition[END_REF] uses a multiobjective approach to combine rewards for performance and novelty. The population is divided into niches based on similarity and their performance judged in relation to other members of their niche. Novelty is judged globally, with individuals rewarded based on their dissimilarity to their neighbors. In this way both exploration of the search space, as well as exploitation of existing niches is accomplished.

The MAP-Elites algorithm [START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] is designed to produce highperforming solutions across a continuum of n user-de ned feature dimensions. It rst divides the feature space into a grid, or map, of n-dimensional bins. The map houses the population of solutions, with each bin holding a single solution. When the map is visualized, with each bin colored according to the performance of the solution it contains, it provides an intuitive overview of the performance potential of each region of the feature space.

To initialize MAP-Elites a set of random solutions are rst evaluated and assigned to bins. The bin location of a solution is based on its features. If, for example, the feature space is 2D with one dimension for weight and another for cost, a low cost and low weight solution would be placed in the low cost, low weight bin location of the map. If the bin is empty, the solution is placed inside. If another solution is already occupying the bin, the new solution replaces it if it has a higher tness, otherwise it is discarded. As a result, each bin contains the best solution found so far for each combination of features. These solutions are known as elites.

To produce new solutions, parents are chosen randomly from the elites, mutated, and then evaluated and assigned a bin based on their features. Child solutions have two ways of joining the breeding pool: discovering an unoccupied bin, or out-competing an existing solution for its bin. Repetition of this process produces an increasingly explored feature space and an increasingly optimal collection of solutions, illuminating the performance potential of the entire feature space. MAP-Elites is summarized in Alg. 1.

MAP-Elites has been shown to be e ective in exploration and optimization in a variety of domains including: the design of walking soft robot morphologies [START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF], the generation of images that fool deep neural networks [START_REF] Nguyen | Deep neural networks are easily fooled: High con dence predictions for unrecognizable images[END_REF], and the evolution of robot controllers capable of adapting to damage [START_REF] Cully | Robots that can adapt like animals[END_REF].

SAIL uses MAP-Elites rather than NSLC for illumination. While the niche de nitions of NSLC are emergent, and neither even or consistent across runs, MAP-Elites de nes a xed structure of feature space boundaries, which greatly simpli es the process of sampling new solutions for inclusion in the surrogate model. Additionally, for design space exploration, this consistency allows designers to easily visualize and compare the e ect of altered constraints and conditions on the feature space.

Data-E icient Exploration, Optimization, and Modeling of Diverse Designs through SAIL GECCO '17, July 15-19, 2017, Berlin, Germany

Algorithm 1 MAP-Elites 1: function MAP E (objecti e_f unction(), X init ial ) 2:
X ← ∅ empty map for genome 3:

P ← ∅ empty map for performance 4:

X ← X init ial , P ← objecti e_f unction(X init ial )

5:

for iter = 1 → I do 6:

x ← random_selection(X)

7:

x ← random_ ariation(x)

8: b ← f eature_descriptor (x ) 9:
p ← objecti e_f unction(x )

10: if P (b ) = ∅ or P (b ) < p then 11: P (b ) ← p , X(b ) ← x 12:
end if

13:
end for 14:

return (X, P) Return illuminated map 15: end function

Surrogate-Assisted Optimization

Evolutionary approaches typically require a large number of evaluations before acceptable solutions are found. In many applications these performance calculations are far from trivial, and the computational cost becomes prohibitive. In these cases approximate models of the tness function, or surrogates, are used in their place.

Surrogate-assisted optimization has been a particularly useful approach in the computationally demanding context of computational uid dynamics [START_REF] Dumas | CFD-based optimization for automotive aerodynamics[END_REF][START_REF] Kc Giannakoglou | Aerodynamic shape design using evolutionary algorithms and new gradient-assisted metamodels[END_REF]. In the context of MAP-Elites, even when evaluations are inexpensive, due to their sheer number surrogateassistance has the potential to accelerate the illumination of the search space dramatically.

Modern surrogate-assisted optimization often takes place within the framework of Bayesian optimization (BO) [START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF][START_REF] Cully | Robots that can adapt like animals[END_REF][START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF]. BO approaches the problem of optimization not just as nding the most optimal solution, but of modeling the underlying objective function in high-performing regions. To estimate the objective function probabilistic models are used, giving each sample a predicted objective value and a certainty in that prediction. New samples are chosen where the model predicts a high objective value (exploitation) and where prediction uncertainty is high (exploration). The relative emphasis on exploitation and exploration is determined by the acquisition function. The sample which maximizes the acquisition function is chosen as the next observation.

A variety of data-driven machine learning techniques such as polynomial regression, support vector machines, and arti cial neural networks, can be used to construct surrogate models [START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF][START_REF]A comprehensive survey of tness approximation in evolutionary computation[END_REF], however as BO requires a probabilistic model, Gaussian processes (GP) [START_REF] Rasmussen | Gaussian Process for Machine Learning[END_REF] are typically used.

Gaussian Process Models. In the presented implementation of SAIL, Gaussian process (GP) models [START_REF] Rasmussen | Gaussian Process for Machine Learning[END_REF] are chosen for tness approximation. GP models are e ective even with a small number of samples and their predictions include a measure of certainty. In the active learning context of surrogate-assisted optimization a measure of model uncertainty is particularly useful, as this allows for the balancing of exploration and exploitation.

Gaussian process models are a generalization of the Gaussian distribution: where a Gaussian distribution describes random variables, de ned by mean and variance, a Gaussian process describes a random distribution of functions, de ned by a mean function m, and covariance function k.

f (x ) ∼ GP (m(x ), k (x, x )) (1)
In much the same way as an arti cial neural network can be thought of as a function that returns a scalar given an arbitrary input vector x, a GP model can be thought of as a function that, given x returns the mean and variance of a normal distribution, with the variance indicating the certainty of our prediction.

Gaussian process models make their predictions based on locality in the input space, a relationship de ned by a covariance function. A common choice is the squared exponential function: the closer the points are in input space the more closely correlated they are in the output space:

k (x i , x j ) = exp - 1 2 x i -x j 2 (2)
Given observations D = (x 1:t , f 1:t ) where f 1:t = f (x 1:t ), we can build a matrix of covariances. In the simple noise-free case we can then construct the kernel matrix:

K =          k (x 1 , x 1 ) • • • k (x 1 , x t ) . . . . . . . . . k (x t , x 1 ) • • • k (x t , x t )          (3)
Considering a new point (x t +1 ) we can derive the value (f t +1 = f (x t +1 )) from the normal distribution (for simplicity we assume a zero mean function m(x ) = 0):

f 1:t f t +1 ∼ N 0, K k k T k (x t +1 , x t +1 ) (4)
where k = [k (x t +1 , x 1 ), k (x t +1 , x 2 ), . . . , k (x t +1 , x t )] T allowing us to compute the GP as:

P ( f t +1 |D 1:t , x t +1 ) = N µ t (x t +1 ), σ 2 t (x t +1 ) (5)
where:

µ t (x t +1 ) = k T K -1 f 1:t (6) σ 2 t (x t +1 ) = k (x t +1 , x t +1 ) -k T K -1 k (7)
gives us the predicted mean and variance for a normal distribution at the new point x t +1 . If we were then to evaluate the objective function at this point, we would add it to our set of observations D, reducing the variance at x t +1 and at other points near to x t +1 . In this pure generalized form, our GP model weighs variations in every dimension equally, applying the same squared exponential relationship regardless of input dimension. For higher dimensional problems each dimension's e ect on the output is also weighted via a technique known as automatic relevance detection (ARD). The hyperparameters which weigh each dimension are set by maximizing the likelihood of the model given the data [START_REF] Rasmussen | Gaussian Process for Machine Learning[END_REF]. This increases model accuracy, and the weighting provides an understandable estimation of the relative importance of each input dimension.

SURROGATE-ASSISTED ILLUMINATION

To understand the relationship between features and performance, SAIL models the underlying objective function in di erent regions of the feature space. Sampling of the objective function in order to model its behavior in the best performing regions is also the goal of Bayesian optimization [START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF][START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF], and we adopt similar methods.

Bayesian optimization has two components. The rst is a probabilistic surrogate model of the objective function, which in SAIL takes the form of a Gaussian process (GP) model (see Section 2.2). The second is an acquisition function, which describes the utility of sampling a given point. The point with maximal utility is evaluated and its performance added to an observation set. The updated set of observations is then used to produce a more informed GP model. As we are not looking to model the objective function only at the global optimum, but at optima in all regions of the feature space, we must produce points which maximize utility in every region of the feature space.

Evaluating new solutions is expensive, making the de nition of "utility" critical to performance. Balance must be maintained between exploration, sampling points with high uncertainty, and exploitation, sampling of points which are likely to perform better than our current solutions.

The acquisition function de nes how the balance between exploration and exploitation is determined. In SAIL, the upper con dence bound (UCB) [START_REF] Srinivas | Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design[END_REF] is used. Proposed as part of the GP-UCB algorithm, use of UCB has been shown to minimize regret and maximize information gain in multi-armed bandit problems [START_REF] Srinivas | Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design[END_REF]. UCB judges potential observations optimistically, favoring uncertainty under the assumption that higher uncertainty hides a potentially higher reward. A high mean (µ (x )) and large uncertainty (σ (x )) are both favored, with relative emphasis tuned by the parameter κ.

UCB(x) = µ (x) + κσ (x) (8)
UCB performs competitively with more complex acquisition functions such as Expected Improvement (EI) and Probability of Improvement (PI) [START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF][START_REF] Calandra | An Experimental Evaluation of Bayesian Optimization on Bipedal Locomotion[END_REF]. These acquisition functions rely on comparisons to the current optimum, while UCB is based only on the underlying model. As SAIL is used to solve numerous localized problems in parallel, it requires an acquisition function independent of the global optimum. If compared globally, solutions in less optimal regions of the design space would have a vanishingly small probability of improving on the global optimum, and as bins are likely not to contain any precisely evaluated solutions, it will not always be possible to perform local comparisons against optima within a bin.

To estimate the relationship between features and performance, SAIL models the objective function not only around a global optimum, but around high-performing solutions over the entire feature space. To accurately predict performance in this slice of the search space, we produce potential observations with every combination of features. By dividing the feature space into bins and using MAP-Elites to produce a solution which maximizes the acquisition function in each, we produce an acquisition map.

It is from the acquisition map that we draw new observations. To reduce uncertainty over the entire feature space we use a Sobol sequence [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] to select which bins to draw the next samples from. Sobol sequences iteratively divide the range into ner uniform partitions, allowing for even sampling across the feature space. In the case that a sampled point results in an invalid solution, the next in the sequence can be drawn. Once evaluated the performance of these samples can be added to our set of observations and a new GP model constructed. A new acquisition map can then be created using this updated model, and the process repeated. GP ← Gaussian_process_model (D)

8:
acquisition() ← UCB(GP (x ))

9:

(X acq , P acq ) = MAP E (acquisition(), X) Select solutions from acquisition map for PE 10:

x ← X acq (Sobol it er ) 11:

X ← X ∪ x, P ← P ∪ PE(x) 12: end for 13: 14: 2) Produce Prediction Map 15: D ← (X, P)

Observation Set: Genome, Performance 16: GP ← Gaussian_process_model (D) 17: prediction() ← mean(GP (x )) 18: (X pr ed , P pr ed ) = MAP E (prediction(), X)

The SAIL algorithm is more precisely de ned in (Alg.2). An initial set of individuals is created using a Sobol sequence [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] to ensure initially even coverage of the parameter space. These individuals and their performance form a set of observations D, which is used to construct a GP model. An empty acquisition map is then created and lled with the individuals from D, along with their utility as judged by the acquisition function. These individuals are taken as the starting population for MAP-Elites (Alg.1) which then illuminates the map as described in Section 2.1: an elite is selected and mutated to produce a child, it is assigned a bin based on its features, and it nally competes for the bin if it is not occupied. This illumination process repeats for a number of iterations, and results in an acquisition map of elite individuals who maximize the acquisition function in their bin.

From the acquisition map we select the next samples for evaluation. To ensure even coverage of the feature space, we again employ a Sobol sequence to direct the sampling, this time producing coordinates in feature space rather than parameter values. These coordinates indicate the bin to be sampled, and the individual stored is precisely evaluated. Once evaluated these new individuals and tness pairs are added to our observation set D and the process can be repeated.

The mean prediction of the resulting GP model can then be taken as the tness function of MAP-Elites, and a prediction map produced. This map is an estimate of the relationship between features and performance, including an optimal design for each bin. As only the surrogate model of the objective is used for evaluation, this prediction map can be produced with minimal computation.

EXPERIMENTAL SETUP 4.1 Objectives and Constraints

We evaluate the performance of SAIL on a classic design problem, 2D airfoil optimization. Fitness is de ned as minimal drag while maintaining the same area and not decreasing lift compared to a base airfoil. Quadratically increasing penalties are introduced into the tness function to ensure that these constraints are followed with little deviation. The high-performing RAE2822 airfoil was chosen as our base, with foils evaluated at an angle of attack of 2.7 • , at Mach 0.5 and Reynolds number of 10 6 . Evaluation criteria are formally de ned for a solution x as:

tness(x ) = drag(x ) × penalty li (x ) × penalty area (x ) (9)
where drag(x ) = -lo (C D (x ))

penalty li (x ) =        C L (x ) li base 2 , if C L (x ) < li f t base 1, otherwise (10)
penalty area (x ) = 1 - |area -area base | area base 7 ( 11
)
While the area of the foil can be directly measured without aerodynamic tests, the drag 1 (C D ) and lift (C L ) must both be approximated. The UCB of the drag prediction is taken as the drag component of our tness function:

dra (x ) = µ drag (x ) + κσ drag (x ) (12)
As individuals are not rewarded for having high lift, but are only expected to maintain performance, we treat the prediction problem as one of classi cation rather than regression. Individuals are penalized based on the probability that they will have a lower lift than our base foil, based on the mean and variance supplied by our GP model:

penalty li (x ) = 1 -P (C L (x ) < li f t base ) (13)

Representation

We encode the airfoil using a variation of the the airfoil-speci c PARSEC parameterization [START_REF] Sobieczky | Parametric airfoils and wings. Recent Development of Aerodynamic Design[END_REF]. PARSEC uses polynomial expressions to encode design features, such as the radius of the leading edge or the curvature of the upper surface, requiring a small number of design parameters to express a large variety of designs.

We restrict the design space to foils with trailing edges which have the same end point and sharpness as our base foil. We also add an additional degree of freedom by splitting the leading edge radius into an upper (r LE up ) and lower leading edge radius (r LE l o ). The ten parameters used to de ne an airfoil are shown in Figure 2.

Dimensions of Variation

Illumination algorithms allow us to de ne dimensions of variation in which we would like to explore. We choose two of our PARSEC descriptors: the height of the highest point on the top side of the foil (Z up ), and the location along the length of the wing of this highest point (X up ). In early tests these parameters were found 1 As C D values are very small, they are converted to log scale in our tness calculation to be highly predictive of the drag. The range of Z up and X up are discretized into 25 partitions, giving us a 25×25 grid, or 625 bins.

In practice the dimensions of variation do not have to be parameter values, and in fact it is desirable that they not be. De ning dimensions of variation which do not align with the representation, but rather correspond to more abstract feature measures, allows for search in a low-dimensional feature space even with a high dimensional representation. Low level features should be chosen based on characteristics that the user would like to explore or, through their own experience, know are important or interesting. In this case parameter values were chosen as dimensions of variation for ease of analysis and comparison with other algorithms.

Zxx lo rLE up rLE lo Z up Z lo α TE β TE X lo X Zxx up

Baseline and Hyperparameters

To evaluate the optimality of the prediction maps produced by SAIL and how e ciently they are produced we compare to 1) standard MAP-Elites without surrogate assistance, and two variants of a traditional convergent search algorithm: 2) the covariance matrix adaptation evolution strategy (CMA-ES), and 3) surrogate-assisted CMA-ES (SA-CMA-ES). The unit of comparison used is the number of precise evaluations (PE), i.e. actual calls to the simulator.

We provide the SAIL algorithm a computational budget of 1000PE. 50PE is used to evaluate the initial pool of individuals which form the basis of the GP model. The remaining 950PE are spent in the course of the algorithm, with 10 new individuals added to the observation set at every iteration (Alg. 2 lines 4-12). This was compared to the standard MAP-Elites algorithm with a budget of 10 5 PE.

We are unaware of any other similar design space exploration techniques and so for a better understanding of the di culty of the task and the optimality of the solutions produced by SAIL we compare to the results of traditional convergent search algorithms, algorithms which are designed to nd a single optimum solution. As we have chosen parameter values as our dimensions of variation, it is possible to con ne a search within one bin of the map by restricting the valid parameter ranges of X up and Z up . Each bin can then be thought of as a single search problem. We solve each of these subproblems with the well-established covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. A budget of 1000PE per bin is given to nd optimal solutions.

A surrogate-assisted variant of CMA-ES (SA-CMA-ES) is also applied to solve the subproblem in each bin. A GP model is produced with 25 initial individuals drawn from a Sobol sequence, sampling in the same way as SAIL. CMA-ES is then used to maximize the acquisition function, computed with the same UCB-based tness criteria as SAIL, described in Section 4.1. The found optimum is added to the set of observations and the optimization process repeats with an updated model. This process is repeated 75 times, for a total of 100PE. Each bin is considered a distinct subproblem, and models and samples are not shared across bins.

Runs of CMA-ES, SA-CMA-ES, SAIL, and MAP-Elites were each replicated 20 times. 2 As optimal performance varies depending on the bin, in some comparisons tness will be reported as a percentage of the optimum value found in all experiments, i.e. 0% -100% of the optimum. Unless otherwise mentioned all values are medians across all experiments. Valid initial designs with a highest point at the leading edge of the wing (high Z up and low X up ) could not be found due to geometric constraints inherent in the PARSEC representation [START_REF] Padulo | Airfoil Design Under Uncertainty with Robust Geometric Parameterization[END_REF]. Only the remaining 577 bins were considered. Beyond our own implementation3 standard implementations were used for CMA-ES [START_REF] Hansen | CMAES Version 3.61[END_REF], Gaussian Processes [START_REF] Rasmussen | Gaussian Process Regression and Classi cation Toolbox[END_REF], and airfoil simulation [START_REF] Drela | XFOIL airfoil simulator[END_REF]. The prediction map of the feature space produced by SAIL in Figure 3 visualizes the e ect of the explored features (X up and Z up ) on performance. The height of the airfoil (Z up ) has the strongest e ect on tness, with taller airfoils performing worse than atter airfoils. The location on the wing of the highest point (X up ) has a more nuanced e ect, increasing or decreasing tness depending on the height of the airfoil. The best performing foils are not at the extremes of the feature space, but at a peak within the mid ranges. Similar designs and trends were also found by CMA-ES.

RESULTS

X up Z up

Accuracy

To evaluate the accuracy of the produced models, after the nal sample was collected a prediction map was produced. Each design in the prediction map was then precisely evaluated and the true C D and C L compared to the prediction of the model. The median results are shown in Figure 4. On the majority of samples the surrogate is reliably accurate, with more than 90% of drag (lo (C D )) predictions and more than 80% of lift (C L ) predictions within 5% of their true value. Drag errors are clustered in the same region of design space, a region where the ow simulator was less likely to converge and produce valid results. The purpose of our models is to estimate performance in the optimal regions of the search space. To test their accuracy in this high tness slice, we measure their ability to predict the performance of the best designs found by CMA-ES in each bin. We compare models built using a naive sampling of the parameter space with a Sobol sequence [START_REF] Niederreiter | Low-discrepancy and low-dispersion sequences[END_REF] to sampling done using acquisition maps produced by SAIL. These acquisition maps are produced by maximizing three di erent acquisition functions: the mean or variance alone, and the UCB, a combination of the mean and variance (see Section: 3). The accuracy of each model's drag prediction on the best design in every bin is then measured at various stages of the sampling process (Figure 5). Mean squared error (log scale) of drag prediction on optimal designs. Models constructed using designs sampled from parameter space using a Sobol sequence or selected from acquisition maps produced with the mean, variance, or the UCB of the prediction.

Accuracy on High Fitness Slice

By concentrating the sampling process on either high-performing solutions or on reducing overall uncertainty we are able to produce better performing models than evenly sampling the parameter space. When both uncertainty and performance are considered, that is, when using the UCB, SAIL produces models that are an order of magnitude more accurate than uniform sampling.

Optimality and E ciency

Though our goal is not to directly compete with algorithms designed to nd one optimal solution, to accurately portray the design space it is critical that the solutions found are representative of the best designs in their region.

We compared the designs found in each bin by SAIL after 1000PE to the best design found by CMA-ES after all 20 runs for 1000PE in each valid bin (≈ 11.5 million PE in total). Figure 6 shows the median values of the prediction map, the true performance of those median designs, the optimal performance found after 20 runs of CMA-ES, and the tness di erence between these optimal values and those found by SAIL. The tness potential of the feature space is well illuminated, with found designs performing within 5% of the optimum in nearly half of bins, and the general relationship between features and performance accurately portrayed. As we have found no similiar design space exploration algorithms beyond MAP-Elites for comparison, to judge the e ciency of our algorithm we turn to convergent search techniques. As CMA-ES was not intended for use across a multitude of subproblems the total number of PEs needed to arrive at an optimized feature map is highly dependent on the number of bins in the map. Therefore we also compare SAIL to the performance of CMA-ES in a single bin. The progress of the di erent approaches is compared in Figure 7.

X up Z up

Single bin performance is taken as the median performance over all bins. Optimization may progress faster or slower depending on the bin, and this gives us a measure of how near an average bin will be to the optimum after a given number of precise evaluations. Map performance is simply this median multiplied by the number of bins. Performance of individuals produced by SA-CMA-ES and SAIL to construct the initial models is set to 0%, with the rst valid performance indicators at 25PE/bin and 50PE (total) respectively.

With the same computational budget required by CMA-ES to nd a near optimal solution in a single bin, SAIL produces solutions of similar quality in every bin.

The acceleration a orded by surrogate modeling has an even more pronounced e ect on the divergent search techniques (MAP-Elites and SAIL) than on the convergent approaches. Incorporating surrogate-assistance into CMA-ES improves performance by an order of magnitude. MAP-Elites, even when given two orders of magnitude more precise evaluations, is still unable to compete with SAIL's performance. Surrogate-assisted optimization allows for estimations of performance to be calculated based on similarity of solutions, a technique which ts neatly into the illumination approach as solutions in close proximity on the map are also likely to perform similarly.

CONCLUSION AND DISCUSSION

The SAIL algorithm produces a model of the objective function in high-performing regions across the feature space despite a limited computational budget. With the knowledge that our models are accurate, we can be con dent in the prediction map's depiction of the feature space, even if the solutions in the map have not been precisely evaluated.

Prediction maps which illuminate di erent feature combinations of the search space can be produced quickly without additional precise evaluations or model training. This allows easy exploration and visualization of the design space through various lenses. Acceleration of the illumination process allows the exploration process to take place in an anytime fashion: as soon as new samples are evaluated, the surrogate model can be reconstructed and estimates of the entirety of the feature space can be rapidly updated.

This assumes, of course, that our models can be trained quickly. In our analysis we concentrated only on the e ciency of the algorithm with regards to precise evaluations. While appropriate in extreme cases, such as uid dynamics, in practice the cost of training surrogate models must be balanced against the savings they yield. In light of the sheer number of evaluations required by MAP-Elites the savings will typically be substantial.

While directing the sampling process with the UCB of the prediction produced more accurate models than using the mean or variance alone, the importance of this improved accuracy is unclear. More investigation is needed into the e ect of di erent acquisition functions and how best to then choose samples from the resulting acquisition map. In the most expensive cases, human-in-the-loop approaches may be appropriate, with experienced designers selecting designs from the acquisition map for evaluation.

In our experiments parameter values served as features, making search within regions of the feature space with a traditional optimizer straightforward. Features are not always so easy to compute, especially if those features are behaviors identi ed during evaluation, as in evolutionary robotics [START_REF] Mouret | Encouraging Behavioral Diversity in Evolutionary Robotics: An Empirical Study[END_REF]. In cases where classifying a solution in feature space is itself expensive, it may be necessary to also construct models to approximate this classi cation.

MAP-Elites grew out of the evolutionary robotics community where it is common to employ representations that themselves evolve and grow more complex, such as NEAT [START_REF] Stanley | Evolving neural networks through augmenting topologies[END_REF]. If SAIL is to be used with non-static representations, like those produced by NEAT, or those that are static but very high dimensional, like those produced by CPPNs [START_REF] Stanley | A hypercube-based encoding for evolving large-scale neural networks[END_REF], specialized surrogate modeling techniques must be developed. Though MAP-Elites has shown remarkable potential, the intensive computation it requires precludes its use in many domains. By pairing MAP-Elites with a surrogate modeling, a Bayesian optimization equivalent for illumination is created. By enabling illumination in computationally expensive domains SAIL opens up new avenues for experiments and applications of quality-diversity techniques.

Figure 1 :

 1 Figure 1: Surrogate-Assisted Illumination (SAIL) 0) Sample design space to produce initial solutions. 1) Construct model of objective based on samples. 2) Maximize the acquisition function, which balances exploitation and exploration, in every region of the feature space, producing an acquisition map. 3) Draw next samples to test on the objective function from the acquisition map. Repeat steps 1-3. 4) Maximize tness as predicted by the resulting model, producing a prediction map.

Algorithm 2 Surrogate-Assisted Illumination 1 : 3 : 4 : 1 )

 21341 Initialize with G solutions drawn from Sobol sequence 2: X ← Sobol 1:G , P ← PE(X) PE = precise evaluation Produce Acquisition Map 5: for iter = 1 → precise_e aluation_bud et do 6: D ← (X, P) Observation Set: Genome, Performance 7:

Figure 2 :

 2 Figure2: The ten parameters used to de ne an airfoil. Dimensions of variation (X up and Z up ) in gold.

Figure 3 :

 3 Figure 3: Design Space Overview with SAIL Prediction map produced by SAIL after 1000PE. Border: Median performing designs found by SAIL in green, best designs found by CMA-ES in black.

Figure 4 :

 4 Figure 4: Drag and Li Predictions Per Bin Predicted and true values of drag (lo (C D )) and lift (C L ) for designs in each bin after 1000PE.

Figure 5 :

 5 Figure 5: Accuracy of Sampling StrategiesMean squared error (log scale) of drag prediction on optimal designs. Models constructed using designs sampled from parameter space using a Sobol sequence or selected from acquisition maps produced with the mean, variance, or the UCB of the prediction.

Fitness

Figure 6 :

 6 Figure 6: Performance of Designs Found By SAIL Top: Median predicted and true performance of designs found by SAIL with a budget of 1000PE compared with performance of optimal designs found by 20 runs of CMA-ES per bin (≈11.5 million PE) Bo om: Optimality of SAIL designs per bin.

  Optimization E ciency in a Single Bin and Over the Entire Design Space Computational e ciency of CMA-ES, SA-CMA-ES, MAP-Elites, and SAIL in precise evaluations. Bin: median progress towards optimum in every bin. Map: performance of CMA-ES and SA-CMA-ES is median bin performance multiplied by number of bins. Performance of individuals produced to construct initial models is set to 0%. Bounds indicate one standard deviation over 20 replicates. PEs and performance in log scale.

				Optimization Performance Per Precise Evaluation
	Percentage of Optimum	10 1 99% 90% 0%	98.9% (Bin) SA-CMA-ES 10 2	98.5% (Bin) CMA-ES 97.6% (Map) SAIL 10 3	98.9% (Map) SA-CMA-ES 10 4	91.9% (Map) MAP-Elites 98.5% (Map) CMA-ES 10 5
					Precise Evaluations
			CMA-ES	SA-CMA-ES	SAIL	MAP-Elites
	Figure 7:

One replicate, including data gathering, with 8 cores of a Intel Xeon 2.6GHz processor required approximately: SA-CMA-ES:32h, CMA-ES:80h, SAIL:12h, MAP-Elites:14h

github.com/agaier/sail_gecco2017

ACKNOWLEDGMENTS

This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement number 637972, project "ResiBots") and the German Federal Ministry of Education and Research (BMBF) under the Forschung an Fachhochschulen mit Unternehmen programme (grant agreement number 03FH012PX5 project "Aeromat"). The authors would like to thank Roby Velez, Alexander Hagg, and the ResiBots team for their feedback.


Mohammed Amin Tahraoui

Eric Duchene
email: eric.duchene@univ-lyon1.fr

Hamamache Kheddouci

Mariusz Wozniak

M-A Tahraoui

E Duchêne

M Woźniak

Graph Theory xx (xxxx) 1-12 LABELED EMBEDDING OF (N, N -2)-GRAPHS IN THEIR COMPLEMENTS

Keywords: Packing of graphs, Labeled packing, Permutation. 2010 Mathematics Subject Classification: 05C70

published or not. The documents may come

Labeled embedding of (n, n -2)-graphs in their complements

Context and definitions

Graph theoretical definitions

All graphs considered in this paper are finite, undirected, without loops or multiple edges. If T is a rooted tree of order n, we define an end vertex as a vertex which does not have any son, and a leaf-parent as a vertex whose all of its sons are end vertices.

Given a positive integer n, the graphs K n , P n and C n will denote respectively the complete graph, the path and the cycle on n vertices. For a graph G, we will use V (G) and E(G) to denote its vertex and edge sets respectively. Given V ′ ⊂ V , the subgraph G[V ′ ] denotes the subgraph of G induced by V ′ , i.e., E(G[V ′ ]) contains all the edges of E which have both extremities in V ′ . If a graph G has order n and size m, we say that G is an (n, m)-graph.

An independent set of G is a subset of vertices X ⊆ V , such that no two vertices in X are adjacent. An independent set is said to be maximal if no independent set properly contains it. An independent set of maximum cardinality is called a maximum independent set. For undefined terms, we refer the reader to [START_REF] Bondy | Graph Theory with Applications[END_REF]. A permutation σ is a one-to-one mapping of {1, . . . , n} into itself. We say that a permutation σ is fixed-point-free if σ(x) ̸ = x for all x of {1, . . . , n}.

The graph packing problem

The graph packing problem was introduced by Bollobás and Eldridge [START_REF] Bollobós | Packing of graphs and applications to computational complexity[END_REF] and Sauer and Spencer [START_REF] Sauer | Edge disjoint placement of graphs[END_REF] in the late 1970s. Let G 1 , . . . , G k be k graphs of order n. We say that there is a packing of G 1 , . . . , G k (into the complete graph K n ) if there exist bijections σ i : V (G i ) -→ V (K n ), where 1 ≤ i ≤ k, such that

σ * i (E(G i )) ∩ σ * j (E(G j )) = ∅ for i ̸ = j
, and here the map

σ * i : E(G i ) -→ E(K n )
is the one induced by σ i . A packing of k copies of a graph G will be called a k-placement of G. A packing of two copies of G (i.e., a 2-placement) is also called an embedding of G (into its complement G). In other words, an embedding of a graph G is a permutation σ on V (G) such that for each edge uv belonging to

E(G), its image σ(u)σ(v) does not belong to E(G).
In the literature, the question of the existence of an embedding of a given graph received a great attention (see the survey papers [START_REF]WoźniakPacking of graphs and permutations-a surveyDiscrete[END_REF]9]). In [START_REF] Burns | Every (p, p -2) graph is contained in its complement[END_REF], full characterizations of all the (n, n -1) and (n, n) embeddable graphs are given. The case of (n, n -2)-graphs was also solved independently in [START_REF] Bollobós | Packing of graphs and applications to computational complexity[END_REF][START_REF] Burns | Every (p, p -2) graph is contained in its complement[END_REF][START_REF] Sauer | Edge disjoint placement of graphs[END_REF]. In particular, it is proved in [START_REF] Sauer | Edge disjoint placement of graphs[END_REF] that any pair of (n, n -2)-graphs can be packed into K n .

In [START_REF] Duchêne | Labeled packing of graphs[END_REF], Duchêne et al. introduced and studied the graph packing problem for a vertex labeled graph. Roughly speaking, it consists of a graph packing which preserves the labels of the vertices. We give below the formal definition of this problem.

Definition [START_REF] Duchêne | Labeled packing of graphs[END_REF]. Given a positive integer p, let G be a graph of order n and f be a mapping from V (G) to the set {1, . . . , p}. The mapping f is called a p-labeled-

packing of k copies of G into K n if there exist bijections σ i : V (G) -→ V (K n ) for 1 ≤ i ≤ k, such that: 1. σ * i (E(G)) ∩ σ * j (E(G)) = ∅ for all i ̸ = j. 2. For every vertex v of G, we have f (v)=f (σ 1 (v))=f (σ 2 (v)) = • • • =f (σ k (v)).
The maximum positive integer p for which G admits a p-labeled-packing of k copies of G is called the labeled packing number of k copies of G and is denoted by λ k (G). Throughout this paper, a labeled packing of two copies of G will be called a labeled embedding of G. It will be denoted by a pair (f, σ).

Remark that the existence of a packing of k copies of a graph G is a necessary condition for the existence of p-labeled-packing of k copies of G. Indeed, it suffices to choose p = 1. Therefore, the result of Sauer and Spencer [START_REF] Sauer | Edge disjoint placement of graphs[END_REF] ensures the existence of a p-labeled packing for (n, n -2)-graphs. An estimation of the labeled packing number of such graphs is the main issue of the current paper.

The following result was proved in [START_REF] Duchêne | Labeled packing of graphs[END_REF]. It gives an upper bound for the labeled packing number of two copies of a general graph.

Lemma 1 (Duchêne et al.,11). Let G be a graph of order n and let I be a maximum independent set of G. If there exists an embedding of G into K n , then

λ 2 (G) ≤ |I| + ⌊ n -|I| 2 ⌋
In [START_REF] Duchêne | Labeled packing of graphs[END_REF], exact values of λ 2 (G) are given when G is a cycle or a path. In almost all cases, the upper bound of the above lemma is reached. More precisely, it is shown that for all n ≥ 6,

λ 2 (P n ) ∈ {⌊ 3n 4 ⌋, ⌊ 3n 4 ⌋ + 1} λ 2 (C n ) = ⌊ 3n 4 ⌋
The case of trees is also considered [START_REF] Tahraoui | Labeled 2-packings of trees[END_REF], but only a lower bound is proposed.

Tahraoui, Duchêne, Kheddouci and Woźniak

Labeled embedding of graphs and permutations

In this section, we give a strong relationship between a labeled embedding and its permutation structure.

A permutation σ of a finite set can be written as the disjoint union of cycles (two cycles being disjoint if they do not have any common element). Here, a cycle (a 1 , . . . , a n ) is a permutation sending a i to a i+1 for 1 ≤ i ≤ n -1 and a n to a 1 . This representation is called the cyclic decomposition of σ and is denoted by

C(σ).
According to this definition, the cycles of length one correspond to fixed points of σ. For example, the cyclic decomposition of the permutation induced by the labeled embedding of T (in Figure 1) is:

{(v 1 ), (v 2 ), (v 3 ), (v 4 ), (v 5 ), (v 6 ), (v 7 ), (v 8 , v 10 ),(v 11 , v 13 ), (v 9 , v 12 )}.
We now recall a fundamental property of labeled embeddings (see [START_REF] Duchêne | Labeled packing of graphs[END_REF]). For any labeled embedding (f, σ) of a graph G, one can remark that the vertices of every cycle of C(σ) share the same label. In other words, the labeled embedding number of G exactly corresponds to the maximum number of cycles induced by an embedding of G. It means that if G admits an embedding with k cycles, then

λ 2 (G) ≥ k.
Although this correlation between labeled embeddings and the permutation's number of cycles was rencently stated, several studies can be found about the permutation structure of an embedding. In particular, the permutation structure of embeddings of (n, n -2)-graphs was investigated by Woźniak in [START_REF] Woźniak | Embedding Graphs of Small Size[END_REF]:

Theorem 2 (Woźniak,94). Let G be a graph of order n, different from

K 3 ∪ 2K 1 and K 4 ∪ 4K 1 . If |E(G)| ≤ n -2, then there exists a permutation σ on V (G) such that σ 1 , σ 2 , σ 3 define a 3-placement of G.
Moreover, σ has all its cycles of length 3, except for one of length one if n ≡ 1 mod 3 or two of length one if n ≡ 2 mod 3.

According to our previous remarks, the above theorem induces the following result in the context of labeled embeddings.

Corollary 3. Let G be a graph of order n, different from K 3 ∪2K 1 and K 4 ∪4K 1 . If |E(G)| ≤ n -2, then λ 2 (G) ≥ ⌊ n 3 ⌋ + n mod 3
In the next section, we will show that the lower bound of Corollary 3 can be improved (including for the excluded graphs).

Main result

We first define the notion of good permutation for a graph.

Definition. Given a graph G, a permutation σ on V (G) is said to be good if

• σ is an embedding of G,
• σ has at least ⌊ 2n 3 ⌋ cycles,

• every cycle of σ is of order at most 2, i.e., for every pair of distinct vertices

u, v of G, if σ(u) = v, then σ(v) = u.
The following lemma will be useful in a special case of our main result.

Lemma 4. For k > 0, the graph kC 3 ∪ 2K 1 admits a good permutation.

Proof. According to the diagram below (Figure 1), first remark that 3C 3 admits a good permutation. Indeed, the numbers inside the vertices correspond to a labeled embedding with 6 labels, with at most two vertices sharing the same label.

G ∈ {2K 1 , 3K 1 , K 1 ∪ K 2 , 2K 2 , K 1,2 ∪ K 1 }.
In each case, one can quickly check that there exist good permutations with at least two cycles. The property still holds for n = 5, where

G ∈ {K 3 ∪ 2K 1 , K 1 ∪ K 1,3 , K 2 ∪ K 1,2 , P 4 ∪ K 1 }.
Good permutations with at least three cycles can be found. Now let n ≥ 6 and assume there exists a good permutation for every (n T or H admits a leaf-parent of degree 2. By symmetry, we may assume that T admits a leaf-parent say x 1 , of degree 2. Let x 0 and x 2 be the two vertices of T such that (x 0 , x 1 , x 2 ) is an induced path of T and x 2 is an end vertex. Let y 1 be an end vertex of H and y 0 its parent. Now consider the graph

′ , n ′ - 2)-graph of order n ′ < n with n ′ ≥ 3. Since G is an (n, n -2)-graph,
G ′ = G -{x 1 , x 2 , y 1 }. Clearly, G ′ is an (n -3, n -5)-graph with n -3 ≥ 3.
Hence the induction hypothesis guarantees the existence of a good permutation σ ′ for G ′ . This permutation can be extended to a good permutation σ for G as follows:

σ(x 1 ) = { y 1 if σ ′ (x 0 ) = x 0 , x 1 otherwise. σ(x 2 ) = { x 2 if σ ′ (x 0 ) = x 0 , y 1 otherwise. σ(y 1 ) = { x 1 if σ ′ (x 0 ) = x 0 , x 2 otherwise. σ(v) = σ ′ (v) if v ∈ V (G ′ )
Since the number of cycles of σ G-G ′ equals two, and they all are of length at most 2, it ensures that σ is a good permutation for G.

Subcase 1.2:

T and H do not have any leaf-parent of degree 2 and T is a star. Let ℓ be a leaf of H and v be its neighbour (of degree at least 3). Let u be the unique vertex of degree at least 3 in T . Now consider the graph

G ′ = G -T -{ℓ, v}.
Clearly, G ′ is an (k, j)-graph with k -2 ≥ j and k ≥ 2. Hence the induction hypothesis guarantees the existence of a good permutation σ ′ for G ′ . This permutation can be extended to a good permutation σ for G by setting

σ(u) = v, σ(v) = u, σ(ℓ) = ℓ and σ(l i ) = l i for all leaf l i of T .
Subcase 1.3: T and H do not have any leaf-parent of degree 2 and T is not a star. Hence T admits two leaf parents, say u and t. Let u 1 , . . . , u k (resp. t 1 , . . . , t k ′ ) be the leaves adjacent to u (resp. t), with k, k ′ ≥ 2. Let ℓ be a leaf of H and v be its neighbour (of degree at least 3). Let x be the non-leaf neighbour of u (it exists since T is not a star and x = t if T is a bistar). See Figure 3 for a better view of these notations. Without loss of generality, we will assume k ≥ k ′ .

u 2 t t 2 t 1 x u k u 1 u t k ′ v l T H Figure 3. Subcase (1.3) If k ≥ 3, consider the graph G ′ = G -{u, u 1 , . . . , u k , t 1 , ℓ}. Hence G ′ is a (n -k -3, n -k -5)-graph satisfying |V (G ′ )| ≥ 3 (since |V (H)| ≥ 4)
, and the induction hypothesis can be applied. Let σ ′ be a good permutation for G ′ . We now define a permutation σ for G as follows:

σ(u) = { t 1 if σ ′ (x) = v, ℓ otherwise. σ(t 1 ) = { u if σ ′ (x) = v, u 1 otherwise. σ(u 1 ) = { ℓ if σ ′ (x) = v, t 1 otherwise. σ(ℓ) = { u 1 if σ ′ (x) = v, u otherwise. σ(u i ) i>1 = u i σ(v) = σ ′ (v) if v ∈ V (G ′ )
This permutation remains good as (k + 1) cycles of lengths 1 or 2 and (k + 3) vertices have been added to σ ′ et G ′ respectively.

If k = k ′ = 2, we first consider the case where x = t (i.e., T is a bistar).

Consider the graph G ′ = G -{u 1 , t 1 , ℓ} and proceed by induction as above. Since ut is an edge of T , we have σ ′ (u) ̸ = σ ′ (t) and u and t cannot be fixed points simultaneously. Therefore, at least u 1 or t 1 can be set as a fixed point in σ and the two remaining vertices in {u 1 , t 1 , ℓ} can be images of each other (thus implying that σ is good).

In the case where x ̸ = t, consider the graph G ′ = G -{u, u 1 , u 2 , t 1 , t 2 , ℓ}, and proceed as previously with the following permutation:

σ(u) =      t 1 if σ ′ (t) = t, t 1 if σ ′ (t) ̸ = t and σ ′ (x) = v, ℓ otherwise. σ(t 1 ) =      u if σ ′ (t) = t, u if σ ′ (t) ̸ = t and σ ′ (x) = v, t 1 otherwise. σ(t 2 ) =      ℓ if σ ′ (t) = t, t 2 if σ ′ (t) ̸ = t and σ ′ (x) = v, t 2 otherwise. σ(ℓ) =      t 2 if σ ′ (t) = t, ℓ if σ ′ (t) ̸ = t and σ ′ (x) = v, u otherwise. σ(u i ) i=1,2 = u i σ(v) = σ ′ (v) if v ∈ V (G ′ )
One can easily check that σ remains good since 6 vertices and at least 4 cycles have been added to σ ′ .

Case 2: |V (T )| ≥ 3 and H = K 1 .

Subcase 2.1:

There exists a leaf-parent, say x, of degree at least 3. Let ℓ be one of its leaves, and y be the unique vertex of H. Now consider the graph

G ′ = G -{x, ℓ, y}, which satisfies |E(G ′ )| ≤ |V (G ′ )| -2.
Hence it admits a good permutation σ ′ by induction hypothesis. A good permutation σ of G can thus be extended from G ′ by setting σ(x) = y, σ(y) = x and σ(ℓ) = ℓ.

Subcase 2.2: All the vertices of T which are adjacent to leaves are of degree 2.

Let x 0 be such a vertex (it exists since |V (T )| ≥ 3), let ℓ 1 be its adjacent leaf, and x 1 its second neighbor. Now let ℓ 2 be a distinct leaf from ℓ 1 in T , and x 2 be its neighbor. If

x 1 ̸ = ℓ 2 and x 1 ̸ = x 2 , then consider G ′ = G -{x 0 , ℓ 1 , ℓ 2 }, which
admits a good permutation σ ′ by induction hypothesis. Then set σ G ′ = σ ′ and

• If σ ′ (x 1 ) = x 1 : σ(x 0 ) = ℓ 2 , σ(ℓ 2 ) = x 0 , and σ(ℓ 1 ) = ℓ 1 . • If σ ′ (x 1 ) ̸ = x 1 : σ(x 0 ) = x 0 , σ(ℓ 2 ) = ℓ 1 , and σ(ℓ 1 ) = ℓ 2 .
One can now easily check that σ is good for Such a vertex exists since n ≥ 6. Consider the graph G ′ = G -{x 0 , x 1 , y}. By induction hypothesis, there exists a good permutation for G ′ , say σ ′ . We set

G. If x 1 = ℓ 2 or x 1 = x 2 ,
σ(x 0 ) = y, σ(y) = x 0 , σ(x 1 ) = x 1 and for every vertex v ∈ V (G ′ ), σ(v) = σ ′ (v),
which defines a good permutation for G. (c) G is the sum of C m (for some m ≥ 4) and two isolated vertices. If m < 8, then Figure 5 shows labeled embeddings corresponding to good permutations. x 2 , σ(x 6 ) = z, and σ(z) = x 6 . For the same reasons as in case (4.2.b), this permutation is good for G.

σ(x 1 ) =      x 1 if σ ′ (x) ̸ = x and σ ′ (y) ̸ = y, x 4 if σ ′ (x) = x, z otherwise. σ(x 2 ) = { x 4 if σ ′ (x) ̸ = x and σ ′ (y) ̸ = y, x 2 otherwise. σ(x 3 ) =      x 3 if σ ′ (x) ̸ = x and σ ′ (y) ̸ = y, z if σ ′ (x) = x, x 4 otherwise. σ(x 4 ) =      x 2 if σ ′ (x) ̸ = x and σ ′ (y) ̸ = y, x 1 if σ ′ (x) = x, x 3 otherwise. σ(z) =      z if σ ′ (x) ̸ = x and σ ′ (y) ̸ = y, x 3 if σ ′ (x) = x, x 1 otherwise. 4 C 4 ∪ 2K 1

Conclusion

Theorem 5 gives a first lower bound about the labeled embedding number of (n, n -2)-graphs. Yet, the computation of the exact value remains an open question, as this bound is not exact for many families of (n, n -2)-graphs. As an example, consider a cycle C n without two edges. Its labeled packing number is at least the one of C n , (i.e., ⌊3n/4⌋). Yet, for any large value of n, we can find an (n, n -2)-graph for which the bound is tight. Indeed, consider G as an union of k disjoint triangles with K 2 ∪ K 1 . The size of a maximum independent set for this graph equals k+2. According to Lemma 1, we have that λ 2 (G) = 2k+2 = ⌊2n/3⌋.

In addition, we mention that this result can be used to study the labeled embedding of (n, n -1)-graphs. One can show for example that the same bound is valid for the union of cycles with a single tree.

Figure 1 .

 1 Figure 1. Good permutation for 3C 3

Figure 2 .Theorem 5 . 3 ⌋

 253 Figure 2. Good permutation for 2C 3 ∪ 2K 1For k = 3, consider the permutation σ corresponding to the labeled embedding of Figure1, and extend it to G by setting σ(u) = u and σ(t) = t. Then σ remains good for G. For k > 3, we can now conclude to the existence of a good permutation for G by pairing good permutations of 3C 3 with a good permutation of rC 3 ∪ 2K 1 where r is in {1, 2, 3}.

Case 1 :

 1 at least two of its connected components are trees. Denote by T and H two trees of G of highest orders such that |V (T )| ≥ |V (H)|. In what follows, we choose to consider T and H as rooted trees. We consider the following four cases: |V (T )| ≥ 3 and |V (H)| ≥ 3. Two subcases are considered as follows: Subcase 1.1:

Case 3 :

 3 then T is either a P 3 or a P 4 . Since n ≥ 6, it implies that G admits at least another connected component which is an (n, n -1) or an (n, n) connected graph with n ≥ 2 (since |E(G)| ≥ 4, there exists at least one edge in G -T ). In other words, this component is either a tree T ′ , or a tree with an edge T ′ ∪ {e}. Let ℓ 3 be a leaf in T ′ . Note that we do not care whether ℓ 3 is adjacent to e or not. By considering G ′ = G -{x 0 , ℓ 1 , ℓ 3 } together with the above permutation where ℓ 2 is replaced by ℓ 3 , we find a good permutation for G. |V (T )| = 2. Let T = (x 0 , x 1 ) and let y be a vertex of degree 2 of G.

Case 4 :Subcase 4 . 2 :

 442 |V (T )| = 1. In this case, G contains isolated vertices (at least two) and non-tree connected components. Two subcases are considered as follows: Subcase 4.1: G has a vertex, say x, of degree at least 3. Let y and z be two isolated vertices of G. Consider the graph G ′ = G -{x, y, z}. The induction hypothesis guarantees the existence of a good permutation σ ′ for G ′ . By puttingσ(x) = y, σ(y) = x, σ(z) = z and for every vertex v ∈ V (G ′ ), σ(v) = σ ′ (v), weget a good permutation for G. The complementary subcase to (4.1), i.e., G is the sum of two isolated vertices and an union of cycles. This case is solved as follows:(a) G = kC 3 ∪ 2K 1 for some k ≥ 1. Lemma 4 allows us to conclude.(b) G has at least one cycle, say H, of order at least 4, and one cycle, say Q, of order at least 3: let (x 1 , x 2 , x 3 ) be an induced path of H, let x 4 be a vertex of Q and z, t be the two isolated vertices of G. Denote by x (resp. y) the neighbor of x 1 (resp. x 3 ) different from x 2 . Note that we may have x = y in the case H = C 4 . See Figure4for a graphical depiction of these notations.

Figure 4 .

 4 Figure 4. Case (4.2.b)

Figure 5 .

 5 Figure 5. Case C m ∪ 2K 1 for m = 4, . . . , 7


Jean-Baptiste Mouret
email: jean-baptiste.mouret@inria.fr

Konstantinos

Konstantinos Chatzilygeroudis
email: konstantinos.chatzilygeroudis@inria.fr

20 Years of Reality Gap: a few Thoughts about Simulators in Evolutionary Robotics Extended Abstract -SimER workshop

  L'archive ouverte pluridisciplinaire

INTRODUCTION

Simulators are an integral part of the Evolutionary Robotics' literature because they allow researchers to explore their ideas more easily than evolution in physical systems. Given the importance of this topic, the community would bene t from a shared culture of simulators and from sharing their toolkits. In this short paper, we describe the requirements of ER from simulators, what we tried, and how we successfully crossed the "reality gap" [START_REF] Jakobi | Evolutionary robotics and the radical envelope-of-noise hypothesis[END_REF][START_REF] Jakobi | Noise and the reality gap: e use of simulation in evolutionary robotics[END_REF] in many experiments.

WHAT SIMULATOR DO WE NEED?

Evolutionary robotics imposes speci c requirements on the simulators:

• Fast: Most ER experiments spend most of their time in the tness function, that is, in running the simulator for many time-steps; it is therefore critical that the simulators are as fast as possible.

• Short cold start: e tness function will be evaluated millions of times, therefore starting a new simulation should take as li le time as possible.

• Fully re-entrant and thread-safe: A common strategy to reduce the duration of ER experiment is to parallelize the tness evaluations (this strategy makes even more sense now that most computers have many cores); however, to do so, the simulators need to be programmed without any hidden state (e.g. static or global variables) and allow many instances to run at the same time (ideally in the same process).

• Stable: Evolution, like other optimization algorithms, tend to push the simulators to their limits and over-exploit "bugs" and instability; it is, for example, common in ER to observe simulated multi-body robots "explode" (break the joint constraint) because it allows the main body to travel a high-speed.

• Accurate: more accuracy is always be er as it reduces the reality gap (see section 4); however, the level of accuracy depends strongly on the type of study (not all arguments/demonstrations need accurate physical simulations).

• Open source: It is critical for reproducibility that (1) other researchers can con rm and expand results without needing to buy a license for some closed-source so ware, (2) that simulators can be updated to run on new systems / operating systems even if the company that made the original so ware does not exist anymore, (3) it is possible to know exactly what computation is performed; for these reasons, scienti c work in ER needs to focus on open source solutions.

• Stable API: As simulators are not the main research focus of the community, we would like to spend as li le time as possible on maintaining them.

ese requirements are o en di erent from those of "classic" robotics [START_REF] Ivaldi | Tools for simulating humanoid robot dynamics: a survey based on user feedback[END_REF]. For instance, simulators used in mechanical design do not need to be thread-safe and can take a long time to start; and they o en favor accuracy over speed. On the contrary, simulators for video-games do not need physical accuracy (only visual accuracy) and are o en primarily optimized for speed.

WHAT DID WE TRY?

Over the years, our team implemented several open source simulators by taking advantage of the available open source physics libraries; all of them have been interfaced with Sferes 2 [START_REF] Mouret | Sferes 2 : Evolvin' in the Multi-Core World[END_REF] (h p: //github.com/sferes2/sferes2), our framework for evolutionary computation:

A B C D E F G H
• Fastsim: A 2-D kinematic simulator for di erential-drive robots that is designed to be fast, simple and exible; it is programmed in C++ -h p://github.com/jbmouret/libfastsim. We tried 3 physics engines: ODE, Bullet, and Dart. ODE is welldocumented and its API (in C) did not change in many years, but it does not include recent developments like a GPU-accelerated pipeline or the Featherstone's solver for multi-body dynamics [START_REF] Featherstone | Robot dynamics algorithms[END_REF][START_REF] Vincent | Impulse-based dynamic simulation of rigid body systems[END_REF]. Bullet was not very documented when we tried it and its API was changing frequently, but it includes more modern ideas.

Dart is a new physics library designed with robots and animation in mind. It has a modern C++11 API and is faster than ODE in our experiments; however, it is under intense developments and changes o en (it seems that it is close to having a more or less xed API and the developers are very responsive).

In addition, we investigated more advanced simulation environments like Gazebo [START_REF] Koenig | Design and use paradigms for gazebo, an open-source multi-robot simulator[END_REF] (h p://www.gazebosim.org) and V-Rep [START_REF] Rohmer | V-REP: A versatile and scalable robot simulation framework[END_REF] (www.coppeliarobotics.com), but they added too much overhead to be e ciently used in our experiments (their communication layer slows the simulation down for simple robots, they can be di cult/slow to restart, and they put a lof of constraints on the interface; it is also not easy to run many instances of Gazebo in parallel).

THE REALITY GAP

In our experiments, no simulator led to any notable boost of accuracy and all of them led to reality gap issues. As a consequence, instead of working on more accurate simulation, we explored many ideas to overcome or bypass the reality gap.

One of our most promising idea is the transferability approach [6, [START_REF] Sylvain Koos | Fast Damage Recovery in Robotics with the T-Resilience Algorithm[END_REF][START_REF] Koos | Online Discovery of Locomotion Modes for Wheel-Legged Hybrid Robots: a Transferability-based Approach[END_REF][START_REF] Sylvain Koos | Crossing the Reality Gap in Evolutionary Robotics by Promoting Transferable Controllers[END_REF][START_REF] Sylvain Koos | e Transferability Approach: Crossing the Reality Gap in Evolutionary Robotics[END_REF][START_REF] Mouret | Crossing the Reality Gap: a Short Introduction to the Transferability Approach[END_REF]: instead of a empting to correct the simulator to make be er predictions, a supervised learning algorithm uses tests on the real robot to learn the limits of the simulator, that is it learns a mapping between behavior descriptors, extracted from the simulator, and a prediction of the accuracy (or con dence) of the simulation. Intuitively, learning these limits is easier than correcting the simulator because [START_REF] Cameron B Browne | A survey of monte carlo tree search methods[END_REF] there are e ects that cannot be easily modeled in a simulator (e.g. aerodynamic e ects in a simulator that does not implement aerodynamics), , especially with a few data points (e.g. 20 experiments on the robot), and (2) many non transferrable ("cheating") behaviors are easy to dismiss (e.g. exploding robots that travel at high speed, unrealistic jumps, etc.).

When the evolutionary algorithm knows the limits of the simulation, it can select individuals with good tness in simulation and high predicted transferability (e.g., with a multi-objective evolutionary algorithm or a constrained evolutionary algorithm). Put di erently, the evolutionary algorithm avoids the parts of the search space that are not accurately simulated. Using this approach, we discovered that current physics engine (e.g. ODE) can be surprisingly accurate for some behaviors. For instance, on a quadruped robot, we observed a di erence of tness of less than 10% for highperforming behaviors evolved with the transferability approach, versus a di erence of more than 100% when we transferred the best individual obtained in simulation to the real robot.

Using the same intuition but a di erent process (the MAP-Elites algorithm [START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF]), we generated thousands of potential solutions (e.g. gaits) in simulation, then selected online using a data-e cient algorithm (here Bayesian optimization [START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF]) the one that works best on the physical robot [START_REF] Cully | Robots that can adapt like animals[END_REF]. We performed our experiments with a damaged, 6-legged robot whereas the simulations were performed with an intact robot. e reality gap was, therefore very important since we added the issues of the simulator to the damage; nonetheless, we were able to nd high-performing gaits for all the damage conditions that we tested in only a dozen of tests on the real robot [START_REF] Cully | Robots that can adapt like animals[END_REF]. In this process, the online algorithm (Bayesian optimization) searches for the solutions that transfer well, that is, that are high-performing both in simulation and in reality.

Combining the transferability approach with "Novelty Search with Local Competition" [START_REF] Lehman | Evolving a diversity of virtual creatures through novelty search and local competition[END_REF] (an illumination [START_REF] Mouret | Illuminating search spaces by mapping elites[END_REF] or quality-diversity [START_REF] Pugh | ality Diversity: A New Frontier for Evolutionary Computation[END_REF] algorithm), we were able to evolve behavioral repertoires so that a 6-legged could walk in any direction [START_REF] Cully | Evolving a behavioral repertoire for a walking robot[END_REF]. Overall, the results showed that only a few short experiments on the physical robot were needed by the algorithm in order to generate a repertoire of controllers that allows the robot to reach every point in its reachable space. In more recent work, we exploited the same idea of repertoire creation (in simulation only with the intact robot) and combined it with online adaptation (via Gaussian process regression [START_REF] Edward | Gaussian processes for machine learning[END_REF]) and replanning (using Monte Carlo Tree Search [START_REF] Cameron B Browne | A survey of monte carlo tree search methods[END_REF]) to achieve fast recovery of a damaged physical 6-legged robot without any human intervention and while taking the environment into account [START_REF] Chatzilygeroudis | Reset-free Trial-and-Error Learning for Data-E cient Robot Damage Recovery[END_REF].

THE MYTH OF THE PERFECT SIMULATOR

Overall, we do not think that simulators will ever be both fast enough for evolution and accurate enough to cross the reality gap easily. e alternative would be to propose a new generation of simulators that can predict a behavior (and therefore a tness value) but also estimate their con dence in this prediction. is con dence measure is what is learned with the transferability approach, but it could be computed or estimated by the simulator1 : a er having simulated a behavior, we could retrieve the tness value (e.g. the distance covered by a walking robot) and the con dence in the tness. Typically, the con dence would be low for highly-dynamic behaviors or for behaviors that highly rely on some features of the environment (e.g. the friction coe cient), and high for more robust, static behaviors.

Such simulators would allow algorithms to avoid over ing the simulator (e.g. using multi-objective evolution) and thus to greatly increase the transferability of behaviors evolved in simulation. At least two approaches could be investigated to add con dence levels in ER simulators: (1) Monte Carlo estimates, that is, running many

Figure 1 :

 1 Figure 1: Examples of simulated and real robots we have used in our research. A. Fastsim (fast 2-D kinematic simulator) [16]. B. E-puck robot in the same maze as in A [16]. C. adruped robot simulated in RobDyn (Bullet) [14-16, 22]. D. Real quadruped robot based on Dynamixel actuators (AX-12) [14-16, 22]. E. iCub robot crawling simulated in our new Dart-based simulator [23]. F. Real iCub robot crawling. G. Hexapod robot simulated in RobDyn (ODE) [4-6, 13]. H. Real hexapod robot [3-6, 13].

is is what we did when we learned a "simulator" (a dynamic model) from scratch for data-e cient policy search in robotics[START_REF] Chatzilygeroudis | Black-Box Data-e cient Policy Search for Robotics[END_REF].simulations with variations of parameters (e.g. variations of the friction level) and measure the mean and variance, or (2) crowd-sourced estimates, that is, building an online database of experiments with a simulator and the result on the real robot (provided that there are enough users of the same simulators).Last, not all ER experiments need accurate, physical simulators. For instance, our work on behavioral diversity[START_REF] Mouret | Encouraging Behavioral Diversity in Evolutionary Robotics: an Empirical Study[END_REF] did not need such a simulator to make our point (which is, that diversity in behavior space needs to be encouraged). As a consequence, we do not always need an accurate simulation; instead, we need exible simulators that can be fast or accurate, depending on the experiment and the stage of the project.

ACKNOWLEDGMENTS

is work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement number 637972, project "ResiBots").