Aurélie Du

Céline Girard

Jacques Morel

Camille Roubille

Toward a holistic management of skin and joint disorders in PsA patients


Alexandre Legay

Thomas Tiennot

Jean-François Gelly

Maxime Sermesant
email: maxime.sermesant@inria.fr

Jean Bulté
email: jean.bulte@ge.com

End-to-end Cardiac Ultrasound Simulation for a Better Understanding of Image Quality

Keywords: Ultrasound, Cardiac Modelling, Probe Design, Image Quality

Ultrasound imaging is a very versatile and fast medical imaging modality, however it can suffer from serious image quality degradation. The origin of such loss of image quality is often difficult to identify in detail, therefore it makes it difficult to design probes and tools that are less impacted. The objective of this manuscript is to present an endto-end simulation pipeline that makes it possible to generate synthetic ultrasound images while controlling every step of the pipeline, from the simulated cardiac function, to the torso anatomy, probe parameters, and reconstruction process. Such a pipeline enables to vary every parameter in order to quantitatively evaluate its impact on the final image quality. We present here first results on classical ultrasound phantoms and a digital heart. The utility of this pipeline is exemplified with the impact of ribs on the resulting cardiac ultrasound image.

Introduction

Simulation of medical images has been an active research area for several years. Synthetic images are often associated with ground-truth information on the underlying anatomy and function used to generate the images [START_REF] Alessandrini | A Pipeline for the Generation of Realistic 3D Synthetic Echocardiographic Sequences: Methodology and Open-access Database[END_REF].

Therefore, it enables to evaluate image processing algorithms [START_REF] De Craene | 3D strain assessment in ultrasound (straus): A synthetic comparison of five tracking methodologies[END_REF] as well as generating data for machine learning. Approaches vary from simulating the whole physics of medical image acquisition to warping an existing image, including all the intermediate combinations.

Physiology: Heart

Acquisition: Probe

Reconstruction: Image

Anatomy: Torso Fig. 1. End-to-end modelling pipeline for echocardiography simulation.

In this work, we developed an original approach integrating a full simulation pipeline in order to leverage image simulation for image quality understanding (see Fig. 1). It includes detailed modelling of both human anatomy and ultrasound probe. This is the first time that such simulations are directly linked to probe parameters in order to understand better the relationship between these parameters and image quality. Moreover, some modifications of our anatomy numerical model such as addition, removal and resizing of organs, tissues and bones can be easily done and their impacts studied.

Cardiac Model

Modelling cardiac electromechanical activity has been an active research area for the last decades, as can be seen in references from recent reviews like [START_REF] Chabiniok | Multiphysics and multiscale modelling, data-model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics[END_REF]. Simulating 3D cardiac function enables to completely control the multi-physics phenomena including electrophysiology, biomechanics and blood flows. In the context of simulating ultrasound images, cardiac modelling enables to generate ground-truth data on position, motion, strain, velocities, etc. for the subsequently generated images.

In this work, we used a detailed segmentation of the heart available from the Visible Human Project [START_REF] Ackerman | The visible human project[END_REF]. We only simulated a static ultrasound image as a proof-of-concept of the whole pipeline. See Fig. 2 for the anatomy of the heart and the imaging slice used in this manuscript.

Torso Model

Once cardiac activity has been simulated, it is necessary to include it within a digital representation of the human torso, as the heterogeneities within the thorax are responsible for an important part of artifacts and echogenicity problems. We used a detailed segmentation of the human torso from the Visible Human project. It enables to introduce various structures of the human body and to see their influence on the simulated images.

Probe Model

Numerical simulations were performed using k-Wave, a k-space pseudospectral method-based solver [START_REF] Treeby | Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method[END_REF]. The k-Wave toolbox solves coupled first-order acoustic equations. An ultrasonic source, a medium (heterogeneous or not) and a sensor should be defined within these equations. In our case, the ultrasonic source and the sensor are modeled by the 3Sc-RS sector probe designed for cardiac applications (General Electric Healthcare, Illinois, USA.) The k-Wave solution enables the definition of every physical property of the probe: the number of piezoelectric elements and their width, height and spatial positioning. It also allows the complete tuning of the emission and reception beam-forming with the definition of the focal distance in azimuth, the focal distance in elevation as well as transmitted and received apodizations.

We could have chosen Field II [START_REF] Jensen | Field: A program for simulating ultrasound systems[END_REF] [6] but k-Wave was selected for two main reasons. Firstly, it defines the medium of propagation as a grid where the physical parameters influencing the propagation of the ultrasonic wave (speed of sound, density, absorption, acoustic non-linearity parameter B/A) can be defined inhomogeneously. In contrast, Field II defines it as a collection of scatterers defined by their positions and amplitudes. The final objective of this project was to study the anatomical causes of the image degradation so we decided to keep a simulation tool that allowed us to model these physical parameters precisely. Secondly, k-Wave offers a free parallel version on GPU whereas the parallel version of Field II is not open-source. Considering that we can use a GPU cluster for the project, k-Wave was chosen.

The input signal used to drive the piezoelectric elements was made of n sinusoidal periods (0.5 < n < 3) at a central frequency f 0 = 1.7 MHz. Results are presented with harmonic imaging (reception at 2f 0 = 3.4 MHz).

Reconstruction

B-mode imaging consists in a gray-scaled representation of the echogenicity of the tissue. A mechanical wave is sent in the body and is back-scattered because of the speed of sound and density heterogeneities. k-Wave models the medium of propagation as a speed of sound and a density distribution mapped on a Cartesian grid. Other parameters like absorption or the non-linearity coefficient B/A can also be defined on this grid. This medium of propagation can be defined in 2D or 3D. The choice made should be studied case by case.

To do 2D modeling is equivalent to ignore the influence of the planes parallel to the imaging plane. Clearly, this hypothesis is directly satisfied if the medium of propagation is invariant along the elevation direction of the probe. The 3D case can be considered to verify this hypothesis or in all other situations where the parallel planes are important. 3D simulations can be made using an optimized version of k-Wave running on Graphics Processing Unit (GPU) which drastically reduces the computation time (up to 10 times faster than the original code on Matlab). This version is currently unavailable in 2D and 2D simulations can only be run in Matlab (The MathWorks Inc., Massachusetts, USA.) Nonetheless, the memory available on current GPU limits the size of the Cartesian grid that has to be modeled. This essential drawback will be discussed in Section 7.

Results

Theoretical validation of the Point Spread Function (PSF)

The Point Spread Function (PSF) of any imaging system is its response to a point source. For B-mode imaging, such a source consists in an infinitesimally small region where the speed of sound and the density distributions slightly differ from the rest of the medium. Recovering the physical PSF of the 3Sc sector probe within our pipeline is a key step towards its validation.

We here present this validation for the 2D case (see the Discussion section for its extension in 3D). Nonetheless, the spatial steps and the number of points used for the simulation grid have a direct impact on the quality of the PSF recovering. In fact, the maximum frequency that can propagate and be solved in the numerical simulations is inversely proportional to the smallest spatial step used. Here we chose a 2008 x 2008 2D grid with respectively 90 µm and 75 µm spatial steps. The corresponding dimensions of the modeled medium are 180 mm (depth) and 150 mm (width). Thus the maximum frequency respecting the Nyquist theorem of two points-per-wavelength is f max = 8.6 MHz.

The emitted signal is made of 2.56 sinusoidal cycles at f 0 = 1.7 MHz. No apodization or filtering is applied to this signal. The imaging sequence is composed of 120 angles between -30 • and 30 • . The modeled medium consists in an homogeneous speed of sound of c 0 = 1540 m.s -1 and an homogeneous density of d 0 = 1035 kg.m -3 , apart from the point source where c 1 = 1600 m.s -1 and d 1 = 1040 kg.m -3 . This point source is made of one single point of the simulation grid aligned with the middle of the simulated probe and placed at 100 mm. An intensity profile is extracted from the 100 mm-radius circle centered in the middle of the 1D-transducer. A log compression is applied and relative intensities are displayed with respect to the orientation of the point with the probe (Fig. 4). Both profiles are normalized to the maximum intensity obtained in the simulation. The simulated profile is in excellent agreement with the theoretical expectation. The agreement (in particular for the five main lobes) is nonetheless sufficient to validate the PSF recovery within our simulations. This result also confirms the relevance of the spatial steps chosen in our study.

Experimental validation on an acoustic phantom

The theoretical validation of the PSF has been completed by an experimental validation of the simulations. B-mode images of an acoustic phantom 403GS LE (Gammex Inc., Middleton, Wisconsin, USA.) have been acquired with the 3Sc probe connected to a commercial console Vivid S70 (General Electric Healthcare, Illinois, USA.) Contrary to more complex modeling (like in Section 6.3), simulating the phantom in 2D is still relevant since it is invariant in one direction (height direction). Elevation imaging plane is thus not impacted by structural heterogeneities. The grid parameters used for the theoretical validation are conserved (90 µm and 75 µm spatial steps, 2008 x 2008 points, 180 x 150 mm dimensions). Here we chose to use an emission frequency at 2.3 MHz and to focus on harmonic imaging (reception at 4.6 MHz). The spatial steps upper defined guarantee 3 points-per-wavelength at 4.6 MHz. Figure 5 provides both the experimental and simulated B-mode images. The agreement between the two images is fairly good. Both pin targets and anechoic cysts' acoustic behaviours are caught up within our simulation (respectively higher and lower reflection than the background). Speckle pattern of the phantom background is also well recovered.

Numerical simulation of a cardiac image in apical view

An echographic image of the numeric heart phantom described in Section 2 has been simulated. In that case, the elevation dimension has to be modeled leading to the realization of simulations in 3D. The imaging slice used here is shown in Fig. 2(b). The grid dimensions are 1024 x 512 x 64 (depth x width x heigth) with respectively 189 µm, 219 µm and 200 µm spatial steps. These grid dimensions correspond to the largest grid that can fit the memory of GPU currently available and used in this study. The corresponding dimensions of the modeled medium are 193 mm (depth), 112 mm (width) and 15 mm (height). Thus, the maximum frequency respecting the Nyquist theorem of two points-per-wavelength is f max = 4.1 MHz in depth, f max = 3.4 MHz in width and f max = 2.3 MHz in height. In that case we study fundamental imaging at f 0 = 2.8 MHz. The simulated image is provided in Fig. 6 (a). Myocardium as well as both ventricles and atria are well recovered. Speckle pattern is coherent near the probe but spatial resolution is slightly too high further away from it (see Discussion). Such an image can be obtained in about 30 minutes.

Finally we provide preliminary results showing how our pipeline can be used to study the impact of anatomical structures on image quality. We have added a modelling of the rib cage on our previous imaging slice (see Fig. 2(c)). The corresponding simulated image is provided in Fig. 6(b). One can see that the rib cage causes a shadowing of the right ventricle that makes disappear its walls. Ribs also cause a blurring of left atrium and left ventricle. This leads to an apparent thickening of the mitral valve as well as a distortion of the atria's shape. These results are coherent with physicians experience during clinical exams.

Nonetheless both images of Fig. 6 display a numerical artefact because of the limited number of points along the height direction. As mentioned above the maximum frequency that satisfies the Nyquist theorem in that direction is f max = 2.3 MHz which is lower than our fundamental frequency at emission (f 0 = 2.8 MHz). Thus a bright point can be seen around the coordinate (100,60) on both images. This is due to aliasing and the non-propagation of frequencies above f max = 2.3 MHz.

Discussion

The accuracy of our simulations has been demonstrated both theoretically and experimentally in the 2D case where smaller spatial steps can be used with reasonable computational costs. The extension to 3D simulations is still limited by the memory size currently available on GPU. Video cards used in this study (GeForce GTX 1080 Ti, NVIDIA Corporation, California, USA) have a memory of 11 GB which can contain a grid made of 2 25 points. Cardiac simulations presented in this study are designed to meet the best trade-off between this memory constraint, the spatial steps needed to propagate relevant frequencies and the actual size of the medium to simulate. To do so, we had to increase the step size in the height direction leading to aliasing and poor spatial resolution far from the probe.

Nonetheless we believe these limitations will soon disappear. First, the memory available on GPU is rapidly increasing. Second, a discussion with the creators of k-Wave lets us think that a new version of the optimized code, distributed over several GPUs, will soon be available. This will lead to smaller computation times and will also allow running simulations on a far larger grid than today, therefore improving precision and complexity. For now, absorption and non-linearities are not taken into account. These two features can directly be set in k-Wave when defining the medium of propagation. Their addition will increase the accuracy of our pipeline.

Moreover, the propagation of shear-waves in bony structures like ribs would also have to be studied. The k-Wave toolbox provides a wave-equation solver based on a visco-elastic definition of the medium of propagation [START_REF] Treeby | Modelling elastic wave propagation using the k-wave matlab toolbox[END_REF]. This solver has not been used in this study but its implementation should be straightforward.

Finally, the cardiac modeling could benefit from many add-ons related to many anatomic or functional diseases that have to be studied. Important computational improvements can even lead to a 4D (3 spatial dimensions and time) modeling of the heart. In such a case, complex dysfunctions like arrhythmia and their translation into an echographic image could be developed.

Conclusion

An end-to-end simulation pipeline including realistic cardiac, torso and ultrasonic probe modelling has been presented. The control over every simulation steps offers a great freedom for investigating the echographic image quality. The simulations' accuracy has been confirmed by recovering a theoretical PSF and the comparison of 2D simulations with experimental B-mode images of an acoustic phantom. Finally, simulated B-mode images of the heart model with and without the rib cage have been created to investigate its impact on the image quality. The inclusion of more advanced anatomical and physiological features together with current computational improvements, making this simulation pipeline a useful tool to investigate the key aspects of image quality.

Fig. 2 .

 2 Fig. 2. (a) 3D cardiac model and probe location. Corresponding 2D imaging slice without (b) and with the rib cage (c).

Fig. 3 .

 3 Fig. 3. 3D torso model (a) and imaging plan (b).

Fig. 4 .

 4 Fig. 4. Simulated (blue line) and theoretical (red line) PSF.

Fig. 5 .

 5 Fig. 5. (a) Simulated B-mode image of the Gammex phantom. (b) Corresponding Experimental B-mode image of the Gammex phantom acquired with the 3Sc probe on the Vivid S70 console.

Fig. 6 .

 6 Fig. 6. Simulated B-mode image of the heart model without (a) and with (b) ribs within the imaging field.


Julien Villemonteix

Emmanuel Vazquez
email: emmanuel.vazquez@supelec.fr

Maryan Sidorkiewicz

Eric Walter

Global optimization of expensive-to-evaluate functions: an empirical comparison of two sampling criteria

Keywords: expected improvement, global optimization, Kriging

In many global optimization problems motivated by engineering applications, the number of function evaluations is severely limited by time or cost. To ensure that each of these evaluations usefully contributes to the localization of good candidates for the role of global minimizer, a stochastic model of the function can be built to conduct a sequential choice of evaluation points. Based on Gaussian processes and Kriging, the authors have recently introduced the informational approach to global optimization (IAGO) which provides a one-step optimal choice of evaluation points in terms of reduction of uncertainty on the location of the minimizers. To do so, the probability density of the minimizers is approximated using conditional simulations of the Gaussian process model behind Kriging. In this paper, an empirical comparison between the underlying sampling criterion called conditional minimizer entropy (CME) and the standard expected improvement sampling criterion (EI) is presented. Classical test functions are used as well as sample paths of the Gaussian model and an industrial application. They show the interest of the CME sampling criterion in terms of evaluation savings.

Introduction

To minimize an expensive-to-evaluate function f , a common approach is to use a cheap approximation of this function, which can lead to significant savings over traditional methods. In this context, global optimization techniques based on Gaussian processes and Kriging (see, e.g., [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty[END_REF]) are often preferred, for they provide an appealing probabilistic framework to account for the uncertainty on the function approximation. Expensive-to-evaluate functions are often encountered in industrial optimization problems, where the function value may be the output of complex computer simulations, or the result of costly measurements on prototypes.

Most Kriging-based strategies proposed in the past few years (see, e.g., [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] and the references therein) implicitly seek a likely value for a global minimizer, and then assume it to be a suitable location for the next evaluation of f . Yet, making full use of Kriging, it is possible to explicitly account for the uncertainty on the global minimizers, and the most likely location of a global optimizer is not necessarily a good evaluation point to improve the accumulated knowledge on the global minimizers.

Based on these considerations, the Informational Approach to Global Optimization (IAGO) strategy recently proposed in [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] evaluates f where the potential for reduction of the uncertainty on the location of the minimizers is deemed to be highest. The entropy of the conditional distribution of the global minimizers is taken as the uncertainty measure, and is approximated using conditional simulations of the Gaussian process modeling f . This approach has two main advantages over classical Kriging-based global optimization methods, such as the Efficient Global Optimization (EGO) algorithm (see [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]). First, it should lead to significant savings on the number of evaluations of f . Second, results under the form of probability distributions are particularly attractive. The purpose of this paper is to evidence the evaluation savings that can be obtained via the use of IAGO.

EGO and IAGO differ only by the sampling criterion used for choosing the next evaluation point. These two criteria, namely expected improvement (EI) for EGO and conditional minimizer entropy (CME) for IAGO, undergo a series of numerical experiments. The first experiments are conducted on four classical test functions. Later on, empirical convergence rates are estimated using sample paths of a Gaussian process. A final comparison is performed on a real-case application to the design of intake ports in the automotive industry, for which a single evaluation of the function to be optimized requires about ten hours of computer time.

The Kriging framework is briefly recalled in Section 2, as well as the definitions of the EI and CME criteria. A brief description of computational aspects of the IAGO approach is also presented. Section 3 reports the empirical comparison of these two criteria. Finally, Section 4 presents conclusions and offers perspectives for future work.

Kriging-based global optimization

Let X, the factor space, be a compact subset of R d and f : X → R be the function to be minimized. The objective is to find x * a global minimizer of f over X when the evaluation of f is expensive. To do so, a cheap model of f (also known as surrogate approximation) based on previous evaluations will be used. Even if deterministic models have been discussed (as in the response surface methodology, see, e.g., [START_REF] Myers | Response Surface Methodology: Process and Product Optimization Using Designed Experiments[END_REF]), it is stochastic models that will retain our attention, and more precisely the Bayesian approach to global optimization (see, e.g., [START_REF] Mockus | Bayesian Approach to Global Optimization: Theory and Applications[END_REF]). In this framework, f is viewed as a realisation (or sample path) of a stochastic process F (F can also be viewed as a Bayesian prior on f ). The distribution of F conditionally to past evaluation results for f is used to design a sampling criterion to be optimized to choose an additional evaluation point for f . When F is Gaussian (we make this assumption in the rest of the paper), the conditional distribution of F at an untried point is also Gaussian with mean and variance that can be obtained analytically using Kriging (prediction based on Gaussian processes has been known for more than 50 years as Kriging in geostatistics and we shall keep to this terminology). Gaussian models and Kriging have been introduced in the field of Bayesian optimization in [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], through the Efficient Global Optimization (EGO) algorithm. Since then (see [START_REF] Mockus | Bayesian Approach to Global Optimization: Theory and Applications[END_REF] for an overview of previous work in the field), Gaussian processes and Kriging have been the object of most publications in the field of Bayesian global optimization, with improvements of the EGO algorithm (see, e.g., [START_REF] Williams | Sequential design of computer experiments to minimize integrated response functions[END_REF] or [START_REF] Huang | Global optimization of stochastic black-box systems via sequential kriging meta-models[END_REF]) and comparative studies (see, e.g. [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] or [START_REF] Sasena | Exploration of metamodeling sampling criteria for constrained global optimization[END_REF]). Our contribution to the field is also based on Kriging.

Linear prediction

In this section, we recall some well-known facts about Kriging on which the rest of the paper is based (for more details, see [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty[END_REF][START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] and the references therein).

Let k(., .) be the covariance function of F, and x be a point in X where F is to be predicted. The mean of F(x) is assumed to be a finite linear combination of known functions p i of x, m(x) = β β β T p(x), where β β β is a vector of fixed coefficients to be computed, and p(x) = [p 1 (x), ..., p l (x)] T . Usually the functions p i are monomials of low degree in the components of x (in practice, their degrees do not exceed two).

Given the vector f n = [ f (x 1 ), ..., f (x n )] T of past evaluations at points in S n = {x 1 , ..., x n } ∈ X n (a sample value of F n = [F(x 1 ), ..., F(x n )] T ), the Kriging predictor F(x) of F(x) is the minimum-variance unbiased linear predictor in the vector space span{F(x 1 ), ..., F(x n )}. It can be written as

F(x) = λ λ λ(x) T F n , (1)
with F n = [F(x 1 ), ..., F(x n )] T , and λ λ λ(x) the vector of Kriging coefficients for the prediction at x. The vector of coefficients λ λ λ(x) is solution of the linear system of equations

K P P T 0 λ λ λ(x) µ µ µ(x) = k(x) p(x) , (2)
with 0 a matrix of zeros, K = (k(x i , x j )) 1≤i, j≤n the covariance matrix of F at all evaluation points in S n , k(x) = [k(x 1 , x), ..., k(x n , x)] T , the vector of covariances between F(x) and F n , and

P =    p(x 1 ) T . . . p(x n ) T    .
The Kriging coefficients at x can thus be computed without evaluating f (x), along with the variance of the prediction error σ2

(x) = k(x, x) -λ λ λ(x) T k(x) -p(x) T µ µ µ(x) , (3)
as these quantities only depend on the covariance of F. Once f has been evaluated at all x i in S n , the prediction of f (x) is the conditional mean of F, given by

f (x) = E[ F(x)|F n ] = λ λ λ(x) T f n , with F n = {F n = f n } the evaluation results. When f is evaluated exactly, Kriging is an interpolation (∀ x i ∈ S n , F(x i ) = F(x i )).
Although noise on the evaluation results could easily be taken into account in the prediction, in what follows, the evaluations are assumed to be noise-free (see [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] for the noisy case). As advocated in [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF], the covariance of F is chosen within the Matérn class of covariance functions (cf. [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] and the reference therein for more details on the choice of a covariance), and the covariance parameters are either set a priori or estimated from the data using the maximum-likelihood method.

After the evaluations in S n , f (x) is viewed as a sample path of F that interpolates the data f n . Such sample paths, known as conditional sample paths, are realizations of F conditionally to F n and are essential to the IAGO approach. They represent all the behaviors that are deemed possible for f given the results of evaluations in S n . Figure 1(a) illustrates the relationships between f , f , σ and the conditional sample paths.

Kriging-based sampling criteria

Among the many sampling criterion available in the literature, we feel that expected improvement (EI), which has been the object of most publications in the field for the last ten years, is the most suited for a comparison with the one we proposed in [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF].

Expected improvement

This sampling criterion corresponds to a one-step optimal strategy given the Gaussian prior F on the unknown function f . Let f * = min x∈X f (x) be the global minimum of f , S n be a set of n evaluation points in X, and consider M n = min x i ∈S n F(x i ) an estimator for f * . For the loss function

L(S n , F) = M n -f * ,
the risk, or expected loss for a candidate point c for the evaluation of f , given the evaluation results f n , is given by

E(L(S n ∪ {c}, F)|F n ) = E(min{M n , F(c)}|F n ) -f * . (4)
One can show that minimizing (4) is equivalent to maximizing the EI criterion as presented for example in [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF], i.e., EI(c

) = E [I(c)|F n ] , (5)
with

I(c) = 0 if F(c) ≥ M n M n -F(c) otherwise .
One can easily rewrite [START_REF] Cover | Elements of Information Theory[END_REF] as

EI(c) = σ(c) uΦ(u) + Φ ′ (u) , (6)
with

u = m n -f (c) σ(c) , m n = E[M n |F n ] = min x i ∈S n f (x i
) the current estimation of the minimum, and Φ the normal cumulative distribution. The new evaluation point is then chosen as a global maximizer of EI(c).

Conditional minimizer entropy

The IAGO approach is based on two complementary principles, that set it apart from previous work in Bayesian global optimization. First, a one-step optimal sampling criterion for the reduction of the uncertainty on the minimizers. Second, the use of Kriging to evaluate this sampling criterion by approximating the distribution of the minimizers conditionally to past evaluations. We now briefly present our sampling criterion, and refer to [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] for computational details.

In [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF], conditional entropy has been introduced to measure the information gained on the minimizers by an additional evaluation of f . This Stepwise Uncertainty Reduction (SUR) strategy [START_REF] Geman | An active testing model for tracking roads in satellite images[END_REF], chooses the point that potentially brings the largest reduction in entropy (seen as a measure of uncertainty).

More formally, given our Gaussian prior F on the function f to be minimized, the uncertainty on the minimizer x * can be measured by the entropy of the global minimizers

H(X * ) = -∑ x∈G p X * (x)log(p X * (x)),
with X * a random vector uniformly distributed in the set of the global minimizers of F over a discrete approximation G of X, and p X * the point mass density of X * . Now, given a vector f n of evaluation results, the uncertainty left on x * is the entropy of p X * (•|F n )

the point mass density of X * conditionally to the evaluation results F n (or in short conditional minimizer density),

H(X * |F n ) = -∑ x∈G p X * (x|F n )log(p X * (x|F n )).
The idea of the IAGO strategy is iteratively to ensure a one-step optimal reduction of the entropy of this distribution.

The risk associated with a candidate evaluation at c ∈ X is then chosen as the differential entropy of the global minimizers conditionally to the potential result of an evaluation at c (in short CME for conditional minimizer entropy)

H n (c) = H(X * |F n , F(c)),
and the evaluation is performed at

x n+1 = arg min c∈X H n (c).
From the definition of conditional entropy [START_REF] Cover | Elements of Information Theory[END_REF], we can write

H n (c) = Z y∈R p F(c) (y|F n ) -∑ x∈G p X * (x|F n , F(c) = y) log(p X * (x|F n , F(c) = y)) dy, ( 7
)
with p F(c) (•|F n ) the distribution of F(c) and p X * (•|F n , F(c) = y) the distribution of X * conditionally to F n and {F(c) = y}.
The CME H n (c), as written in [START_REF] Hartman | Some experiments in global optimization[END_REF], can be viewed as an expected loss, the loss function being the entropy of p X * (•|F n , F(c) = y) the conditional minimizer density after n + 1 evaluations.

Practical aspects

The distribution p F(c) (•|F n ) is Gaussian, with mean and variance simply obtained by Kriging. There is, however, no result in the literature that we can use to describe analytically any useful property of the conditional minimizer density. To compute [START_REF] Hartman | Some experiments in global optimization[END_REF], we resort to an approximation that is conducted via Monte-Carlo simulations of F conditionally to available evaluation results F n and to a potential evaluation result y at c (this approximation as well as recommendations for the choice of G are described in details in [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF]). This approximation leads to a complexity in O(N) for the computation of H n (c), with N the size of the discrete approximation of X. Note that in IAGO the conditional minimizer density is thus available at each step and provides (at least for low-dimensional problems) a clear view of the progress achieved in the optimization process (cf. Figure 1(b)).

In the Bayesian optimization framework, the expensive-to-evaluate function is replaced by a cheap criterion, updated after each evaluation, which has to be optimized for a new evaluation point to be chosen. Up to now, we have focused on the choice of criterion, but no attention has been paid to the entire procedure for global optimization, including for example an update process for the Kriging prediction. To keep this paper focused on a comparison between sampling criteria, we shall only mention the classical framework of the Efficient Global Optimization (EGO) (Algorithm 1).

EGO (see, e.g. [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]) starts with a small initial design used to get a first estimate of the parameters of the covariance and to compute a first Kriging model. Based on this model, an additional point is selected in the design space to be the location of the next evaluation of f in order to maximize the EI criterion. The parameters of the covariance are then re-estimated, the Kriging model is re-computed, and the process of choosing new points continues until the improvement expected from sampling additional points has become sufficiently small. The CME criterion can easily be inserted in a similar algorithm in place of EI to transform EGO into IAGO.

Algorithm 1: Efficient global optimization framework

Input: Initial design of evaluation points and corresponding values of f Output: Additional evaluations 1. while the evaluation budget is not exhausted or some other convergence condition is not satisfied 2.

do Estimate the parameters of the covariance 3.

Compute the Kriging model 4.

Optimize the sampling criterion (EI or CME here) 5.

Evaluate f

Empirical comparison between EI and CME

As presented in the previous section, EI and CME are both Kriging-based sampling criteria and both onestep optimal in some sense. CME should lead to faster convergence rates, and this for three major reasons. First, EI aims at estimating the minimum, while CME concentrates on the minimizers. The search is therefore likely to be more global when based on the latter. Second, EI aims at improving the estimation of the minimum by sampling where its appearance is most probable. It seems more reasonable to try diminishing the uncertainty associated with its position. For example, it might be excessively costly to refine the estimation in a small neighborhood of a potential minimum, which may only be local, while evaluations using H n could show that a large part of the search space has a very low probability of containing the global minimum (this idea will be confirmed in Section 3.3). Third, the computation of CME involves the statistical properties of the sample paths of F, while, by contrast the computation of EI involves only the conditional mean and variance of F at c. A more thorough use of the available information on the function is indeed appealing in this context of expensive, and therefore sparse, evaluations.

To substantiate these intuitions, a comparison of EI and CME is in order.

Experimental conditions

To make this comparison fair, we propose to study the behaviors of EI and CME independently Algorithm 1) and from the optimization method to be used to optimize the sampling criteria (Step 4 in Algorithm 1). These aspects are quite complex, and ad-hoc strategies have been proposed in the literature (see [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] for an example of optimization method for the sampling criteria). However, our first objective here is to motivate the choice of a sampling criterion. Therefore, we conducted our experiments using the same Matérn covariance with the same values, fixed a priori, for the covariance parameters. The set G of potential evaluation points was identical for both criteria, and the choice of the next evaluation point was carried out via an exhaustive computation of the relevant sampling criterion over this set. The question to be addressed in what follows then boils down to the following interrogation: Given the same prior information on the function, which sampling criterion chooses the best point (in a sense to be discussed later) amongst a finite set of possible evaluations points?

Tests on classical benchmarks

The four test functions used in this section are taken from [START_REF] Huang | Global optimization of stochastic black-box systems via sequential kriging meta-models[END_REF], where a comparison was conducted between EI and classical global optimization schemes such as DIRECT (see, e.g., [START_REF] Perttunen | A computational geometric approach to feasible region division in constrained global optimization[END_REF]). The problem dimensions range from two to five, and all functions present several local minimizers (see Table 2 in Appendix 2). The covariance parameters are estimated beforehand on the results of 200 evaluations randomly chosen in search space (using a latin hyper cube sampler), and the two criteria are optimized over a latin hyper cube design containing 1000 points randomly re-sampled after every evaluation.

A single point x 1 is randomly chosen in search space as a common starting point for both criteria, and 50 runs are conducted for each function to reduce the dependency on the starting point. After the i-th evaluation of f , the efficiency of each criteria is measured by

G i = f (x 1 ) -m i f (x 1 ) -f * ,
with m i = min x∈{x 1 ,...x i } f (x) the current estimate of the global minimum. G i (a modified version of the quality measure used in [START_REF] Barton | Minimization algorithms for functions with random noise[END_REF]) thus describes the reduction, after i iterations of the optimization process, of the initial estimation error for the global minimum f (x 1 )f * . Table 1 presents, for each criterion, the averaged efficiency after 20, 50 and 100 evaluations. EI beats CME for the Ackley function when i = 40, but for the other three test functions CME converges faster towards the optimum than EI, and significantly so for the Hartman 3 function.

Tests on Gaussian processes simulations

Even if a comparison on classical test functions gives some perspectives on the qualities of each of the criteria, the variability of the results from one test problem to the next may be significant, so one can hardly use them to decide beforehand which sampling criterion to use on a specific problem. It therefore would be best to derive some analytical convergence rates for both criteria under reasonable hypotheses on the function to be optimized. In our context of expensive-to-evaluate functions, these convergence rates would have to be non-asymptotic, and we do not know of any such results in the literature. However, the probabilistic framework considered here makes it possible to estimate empirical convergence rates. Since the function to be optimized is assumed to be a sample path of a Gaussian process, we can estimate the convergence rates with both criteria when optimizing sample paths of a Gaussian process whose covariance is the same as that chosen for the optimization algorithm.

For the sake of brevity, we shall limit our presentation to two Gaussian processes, one with very smooth sample paths, and the other with irregular sample paths.

Two sets of 1000 sample paths were generated over a regular grid of 1500 points in [0, 1] 2 . 15 evaluations are then performed on each sample paths using both criteria. After each new evaluation, and for each criterion, estimation errors are computed for the global minimum and the minimizer, as well as the entropy of the conditional minimizer density. Two estimators of the global minimum are considered here, namely m n = min x∈{x 1 ,...x n } f (x), the best evaluation result obtained so far, and

mn = f (arg max x∈G p X * (x|F n )),
the predicted value associated with the point where the conditional minimizer density is the highest. The average convergence rates for irregular sample paths are presented on Figure 2 As expected, since the entropy of the conditional minimizer is the loss function behind CME, CME performs significantly better than EI in terms of the entropy of the conditional minimizer density and, in average, the uncertainty on the positions of the global minimizers diminishes faster if points are chosen Fig. 3 Convergence rates using EI (dotted line) and IAGO (bold line), when convergence is measured by H(X * |F n ) (a), and when convergence is measured by the estimation error for the minimum with mn as estimator (b). The dashed line represents, as reference, the convergence rate for a random choice of evaluation points. The sample paths used here are smoother than those used for Figure 2 (the parameter for the Matérn covariance are ν = 5, ρ = 0.3 and σ = 1).

using CME (cf. Figure 2(a)). This fact was guaranteed for the first evaluation since CME is one-step optimal for this loss function, but it had to be checked for several evaluations. Similarly, if the convergence is measured by the estimation error m nf * , EI is bound to perform better if we consider only the first evaluation, since the convergence measure is the loss function behind EI. However, it appears that after 15 evaluations, the performance of EI and CME are similar (cf. Figure 2(b)), suggesting that even if EI is one-step optimal, in the long run, CME will bring the largest reduction for m nf * (this is confirmed by computations, not presented here, with a larger number of evaluations).

EI would thus seem to be a better criterion in a context where very few evaluations are allowed. However, m n is estimator actually a rather poor estimator of the global minimum, and it appears that when a faster-toconverge estimator is used instead of m n , CME performs significantly better than EI (Figure 2(c)), and this right form the start. This estimator is mn , whose interest is apparent for the three search strategy considered here, for which m nf * is significantly bigger than mnf * , and more than three times so after the first evaluation (Figure 2(b) and Figure 2(c)). CME should therefore be preferred to EI when one is confronted with irregular sample paths, since it allows a better estimation of f * .

If we look at what happens on a typical sample path (see Figure 4), the drawbacks of EI are clearly evidenced. As intuitively stated at the beginning of the section, EI stalls on a local optimum because with (4) as a loss function, it is better to ensure a small improvement near a minimum already found than to check that it is effectively a global minimum. In the case of irregular sample paths this might be highly dangerous, and IAGO performs better simply because it first addresses the question of whether a minimizer is global before improving the precision on its exact position. When the sample paths are more regular, this advantage diminishes (cf. Figure 3), as the local optima are scarcer. A significant problem is left aside here, namely what happens in practice if the parameters of the covariance are poorly estimated? Does the optimization strategy still perform well? Robustness to a poor choice of covariance parameters is of course a major issue but it is not considered here.The EI and CME criteria should have similar robustness properties and may both be deceived by a poor choice of covariance as demonstrated in [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]. We feel that this problem should be tackled from a Bayesian point of view, with some prior on the covariance parameters. This will be done in future work, where we shall compare an extended version of IAGO to the methods in [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] designed to be robust to a poor choice of covariance.

Test on an industrial application: intake port design

This section presents an industrial optimization problem in the automotive field, also used for the comparison of CME and EI.

Problem description

Intake ports (Figure 5) are engine components that convey a mixture of air and fuel to the combustion chambers. The importance of this type of component lies in the properties of the flow it induces in the combustion chamber, which has a direct impact on both the performance and the emissions of pollutant by the engine. To comply with new emission standards (Euro V and Euro VI), while satisfying the ever increasing need for engine performance, the shape of intake ports has to be carefully optimized. Two oftenconflicting objectives have to be maximized simultaneously, namely the flow rate and a scalar characteristic of the turbulent flow known as tumble [START_REF] Lumley | Engines: An Introduction[END_REF]. Physics tells us that the higher the flow rate, the larger the amount of fuel that can be burnt, and consequently the larger the power delivered by the engine. Similarly, pollutants as nitrogen oxides (NO x ) and carbon monoxide (CO) are, to a large extent, created when the air/fuel mix is not homogeneous. Therefore, the larger the turbulence (and tumble accounts for the relevant properties of it), the smaller the pollution.

The specifications for these two objectives are liable to change during conception. Therefore, it is important to determine not only an optimal geometry for a given set of preferences but rather the full Pareto front. However, building prototypes for tests is exceedingly expensive, and each flow simulation by finite-element methods takes about ten hours on powerful servers. The approach advocated in this paper is therefore particularly attractive given the general will for reduction of duration and cost associated with development.

Computational issues

To extend our sampling criteria to a multi-objective problem, we use a standard procedure and consider several linear combinations of the objective functions (a.k.a. aggregations), each accounting for a different zone of the Pareto front. In [START_REF] Knowles | Parego: A hybrid algorithm with on-line landscape approxiamtion for expensive multiobjective optimization problems[END_REF], this approach has been used to extend the EI criterion to a multi-objective problem by randomly selecting a new aggregation after each evaluation of f . In this paper, we follow the same route, but use the IAGO framework to compute the entropy of the conditional density of the minimizers for the mono-objective optimization problem corresponding to each aggregation in a given set. The search can thus be directed towards the most uncertain regions of the Pareto front.

The resulting multi objective extensions of CME and EI have been applied to the optimization of six shape parameters of an intake port (these parameters are not detailed here for confidentiality reasons). To improve the number of simulations achievable in a given time, the geometry and mesh are automatically generated for each finite-element simulation. The optimization algorithm is then directly interfaced with the solver, limiting human intervention to the initialization of the procedure.

The initial value for the parameters of the Matérn covariance are estimated on simulations that have been collected during the design of previous intake ports. It thus becomes possible to initialize the algorithms with very few randomly chosen points (here with five points).

One thousand candidate points are used (N = 1000) and the parameters of the Matèrn covariance are fixed.

Results

For comparison purposes, simulations were conduced on twenty intake ports whose parameters were chosen using EI, CME, or a Latin Hyper Cube (LHC) as a reference. The Pareto-optimal points within each of the three sets of evaluation results are presented in Figure 6. Comparison between sets estimates of Pareto fronts is a tricky process, which may involve various quality measures [START_REF] Knowles | A tutorial on the performance assessment of stochastive multiobjective optimizers[END_REF]. Here however, the comparison is clearly in favor of CME, as Fig. 6 Results of 20 evaluations randomly chosen using an LHC (dots), or optimally chosen by CME (crosses) or by EI (squares). For CME and EI, only the Pareto-optimal points are presented. For LHC, a larger dot size indicates points that are, within the considered set of points, Pareto-optimal. Tumble and flow rate have been rescaled for confidentiality reasons. The dashed lines delimit the set of points that are dominated by the Pareto-optimal points obtained using CME and which dominate the "worst point", i.e. the point with worst values attained for both objectives .

-among all evaluation results, the point closest to an "ideal" solution (i.e. with best value yet obtained for both objectives, here [-0.97, -1.97] T ) has been chosen by CME; -all but one point chosen by EI are dominated by points chosen by CME. In other words, almost any good solution found by EI is bettered by a solution found by CME; -the volume of the set of points dominated by the Pareto-optimal points (cf. [START_REF] Knowles | A tutorial on the performance assessment of stochastive multiobjective optimizers[END_REF] for details on this quality measure) is 0.31 for CME (this volume is represented on Figure 6), while it is only 0.26 for EI (the reference for the computation is the point with worst values attained for both objectives as coordinates).

This test case confirms that Bayesian global optimization is easily applicable to an industrial problem, even with a very small evaluation budget. The interest of CME is apparent after only a few evaluations, as predicted by the convergence rates of Section 3.3.

Computational burden

The comparison made so far dealt only with convergence rates. The superiority of CME over EI was demonstrated for sample paths of Gaussian random processes, at least for two very different regularities of the sample paths (cf. Section 3.3), while the convergence rates are generally in favor of CME when applied to some classical test functions (cf. Section 3.2) or to an actual industrial problem. However, EI is easier to compute since it only requires the mean and variance of the prediction at the candidate point, while the complexity of the computation of IAGO is in O(N), with N the size of a discrete approximation of X used for the estimation of the conditional density.

In practice, with our implementation of IAGO (cf. [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF]) around 40s are required on an AMD opteron 285 server to choose an additional evaluation point for the sample paths of Section 3.3 (by extensive computation of H n over 1500 candidate points, which is enough in practice since this set is randomly re-sampled after each evaluation). By comparison, choosing a point with EI takes less than half a second under the same conditions. To broaden the range of potential applications, we tried to limit the computational expense by testing other approximations for the conditional minimizer density (since the manipulation of sample paths, necessary for the approximation proposed here, is responsible for most of the computational burden). We proposed for example, to estimate the derivatives of f by Kriging (as in, e.g., [START_REF] Vazquez | Estimating derivatives and integrals with kriging[END_REF]) and to compute the probability for a given point to be a local optimum and under a certain threshold. It was then easy to build a relatively accurate approximation of the conditional minimizer density, but the approximation had a detrimental effect on the convergence rate, so that EI then became more efficient. In fact, the quality of the approximation of the conditional minimizer density is important for CME to perform better than EI.

IAGO therefore remains destined to the optimization of functions that require a large amount of computer time (or more generally a significant expense) to be evaluated, which is after all what it was designed for and is the case in many applications in the industrial world including the one we presented here.

Conclusions and perspectives

In this paper, we have evidenced a clear superiority of CME over EI, especially when the function to be optimized is irregular. The comparison has been conducted using classical test functions and an actual industrial application, but above all using sample paths of the model behind Kriging. The use of sample paths indeed allows the computation of empirical convergence rates that can also be useful to tune other components of any Kriging-based algorithms (e.g., the optimization of the sampling criterion). Now that the interest of the CME criterion has been demonstrated, attention should turn to other crucial aspects of any Kriging-based optimization algorithm. One of these aspects is the improvement of robustness against a bad choice of covariance for the Gaussian process model. We feel that the Bayesian framework that we have used until now should also be useful in this respect.

Appendices

Appendix 1: Matèrn covariance

In this paper, we follow [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF] and use of the isotropic Matérn covariance:

k(x, y) = k(h) = σ 2 2 ν-1 Γ(ν) 2ν 1/2 h ρ ν K ν 2ν 1/2 h ρ ∀ (x, y) ∈ X 2 , (8)
with h the Euclidean distance between x and y, and K ν the modified Bessel function of the second kind [START_REF] Yaglom | Correlation Theory of Stationary and Related Random Functions I: Basic Results[END_REF].

The parameters of this covariance are easy to interpret, as ν controls regularity, σ 2 is the variance (k(0) = σ 2 ), and ρ represents the range of the covariance, i.e., the characteristic correlation distance. They can either be fixed using prior knowledge on the system, or be estimated from experimental data. In geostatistics, estimation is carried out using the adequacy between the empirical and model covariance (see, e.g., [START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty[END_REF]).

In other areas, cross validation (cf. [START_REF] Wahba | Support vector machines, reproducing kernel Hilbert spaces, and randomized GACV[END_REF]) and maximum likelihood (cf. [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]) are mostly employed. For simplicity and generality reasons (cf. [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]), the maximum-likelihood method is preferred here.

Appendix 2: Test functions

Fig. 1

 1 Fig. 1 (a): Conditional sample paths of F, and corresponding Kriging prediction. The squares represent available values of f , the bold line is the conditional mean f as computed by the Kriging predictor, the dotted lines provide 95% confidence intervals for the prediction ( f ± 1.96 σ) and the thin lines are conditional sample paths. (b): Estimated conditional minimizer density (p X * (•|F n )) associated with the Kriging prediction.

  (a) in terms of the entropy of the conditional minimizer density, on Figure 2(b) in terms of the estimation error for m n , and on Figure 2(c) in terms of the estimation error for mn .

Fig. 2

 2 Fig.2Convergence rates using EI (dotted line) and IAGO (bold line), when convergence is measured by H(X * |F n ) (a), by the estimation error for the global minimum with the best value obtained so far (m n ) as an estimator (b), by the estimation error for the global minimum with mn as an estimator (c). The convergence measures are averaged over 1000 sample paths of a Gaussian process with a Matérn covariance with parameters ν = 1, ρ = 0.3 and σ = 1 (see Appendix 5.1). The dashed line represents, as a reference, the convergence rate for a random choice of evaluation points.

Fig. 4

 4 Fig.[START_REF] Chilès | Geostatistics, Modeling Spatial Uncertainty[END_REF] Minimization of a sample path from the Gaussian process used to evaluate the convergence rates of Figure2. The dots indicate the evaluation points chosen by EI. The crosses indicate the evaluation points chosen by CME. The order in which the evaluations are carried out is also indicated for CME.

Fig. 5

 5 Fig. 5 Intake port. The component itself is in the middle. Below is the combustion chamber. The upper cylinder is a tranquilizing volume necessary for the convergence of finite-element simulations.

d = 2 f (x) = x 2 -π x 1 -6 2 +

 2212 10 1 -1 8π cos x 1 + 10 + 0.5x 1 -5 ≤ x 1 ≤ 10,0 ≤ x 2 ≤ 15 N local = 3,N global = 1 x * = [-3.2,12.3] T f * = -1.17

Hartman 3 [ 7 ] d = 3 f0Ackley 5 [ 1 ] d = 5 f

 73515 (x) = -∑ 4 i=1 d i exp -∑ 3 j=1 α i j (x jp i j ) , ≤ x i ≤ 1, i = 1,2,3 N local > 1,N global = 1 x * = (0.114,0.556,0.852) T , f * = -3.86 (x) = -20 exp[-0.2 1 d ∑ d i=1 x 2 i ]exp[ 1 d ∑ d i=1 cos(2πx i )] + 20 + e ∀i ∈ 1,3 -32.8 ≤ x i ≤ 32.8 N local > 1,N global = 1 x * = 0, f * = 0

Table 1

 1 Comparison of expected improvement and conditional minimizer entropy on four test functions taken from[8] * . For each criterion, the convergence measure G i is averaged over 50 runs (the estimated standard error for the estimation of these figures is always smaller than 0.01).

			G i when points are chosen using EI		G i when points are chosen using CME
		i = 20	i = 50	i = 100	i = 20	i = 50	i = 100
	Six-Hump Camel Back	0.65	1	1	0.76	1	1
	Tilted Branin	0.83	0.92	0.98	0.89	0.95	0.97
	Hartman 3	0.64	0.98	1	0.82	0.99	1
	Ackley 5	0.36	0.75	0.73	0.34	0.59	0.72

*

Table 2

 2 Test functions to be minimized [8] * . -1.6 ≤ x 1 ≤ 2.4,-0.8 ≤ x 2 ≤ 1.2 N local = 6,N global = 2 x * = [0.089,0.713] T and [0.089,0.713] T , f * = -1.03

	Name	Description
		d = 2
		f (x) = 4x 2 1 -2.1x 4 1 + 1/3x 6 1 + x 1 x 2 -4x 2 2 + 4x 4 2
	Six-Hump Camel Back [3]
	Tilted Branin [8]