.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_numerical_methods/optimization/plot_ego.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_numerical_methods_optimization_plot_ego.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_numerical_methods_optimization_plot_ego.py:


EfficientGlobalOptimization examples
====================================

.. GENERATED FROM PYTHON SOURCE LINES 7-19

The EGO algorithm [jones1998]_ is an adaptative optimization method based on
Gaussian Process metamodel.

An initial design of experiment is used to build a first metamodel.
At each iteration a new point that maximizes a criterion is chosen as
optimizer candidate.

The criterion uses a tradeoff between the metamodel value and the conditional
variance.

Then the new point is evaluated using the original model and the metamodel is
relearnt on the extended design of experiments.

.. GENERATED FROM PYTHON SOURCE LINES 22-29

.. code-block:: Python

    from openturns.usecases import branin_function
    from openturns.usecases import ackley_function
    import openturns as ot
    import openturns.experimental as otexp
    import openturns.viewer as viewer









.. GENERATED FROM PYTHON SOURCE LINES 30-34

Ackley test-case
----------------

We first apply the EGO algorithm on the :ref:`Ackley function<use-case-ackley>`.

.. GENERATED FROM PYTHON SOURCE LINES 36-38

Define the problem
^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 40-41

The Ackley model is defined in the usecases module in a data class `AckleyModel` :

.. GENERATED FROM PYTHON SOURCE LINES 41-43

.. code-block:: Python

    am = ackley_function.AckleyModel()








.. GENERATED FROM PYTHON SOURCE LINES 44-45

We get the Ackley function :

.. GENERATED FROM PYTHON SOURCE LINES 45-47

.. code-block:: Python

    model = am.model








.. GENERATED FROM PYTHON SOURCE LINES 48-49

We specify the domain of the model :

.. GENERATED FROM PYTHON SOURCE LINES 49-53

.. code-block:: Python

    dim = am.dim
    lowerbound = [-4.0] * dim
    upperbound = [4.0] * dim








.. GENERATED FROM PYTHON SOURCE LINES 54-55

We know that the global minimum is at the center of the domain. It is stored in the `AckleyModel` data class.

.. GENERATED FROM PYTHON SOURCE LINES 55-57

.. code-block:: Python

    xexact = am.x0








.. GENERATED FROM PYTHON SOURCE LINES 58-59

The minimum value attained `fexact` is :

.. GENERATED FROM PYTHON SOURCE LINES 59-62

.. code-block:: Python

    fexact = model(xexact)
    fexact






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=1 values=[4.44089e-16]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 63-67

.. code-block:: Python

    graph = model.draw(lowerbound, upperbound, [100] * dim)
    graph.setTitle("Ackley function")
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_001.svg
   :alt: Ackley function
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_001.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 68-69

We see that the Ackley function has many local minimas. The global minimum, however, is unique and located at the center of the domain.

.. GENERATED FROM PYTHON SOURCE LINES 71-79

Create the initial GP
^^^^^^^^^^^^^^^^^^^^^

Before using the EGO algorithm, we must create an initial GP metamodel.
In order to do this, we must create a design of experiment which fills the space.
In this situation, the :class:`~openturns.LHSExperiment` is a good place to start (but other design of experiments may allow one to better fill the space).
We use a uniform distribution in order to create a LHS design.
The length of the first LHS is set to ten times the problem dimension as recommended in [jones1998]_.

.. GENERATED FROM PYTHON SOURCE LINES 81-90

.. code-block:: Python

    listUniformDistributions = [
        ot.Uniform(lowerbound[i], upperbound[i]) for i in range(dim)
    ]
    distribution = ot.JointDistribution(listUniformDistributions)
    sampleSize = 10 * dim
    experiment = ot.LHSExperiment(distribution, sampleSize)
    inputSample = experiment.generate()
    outputSample = model(inputSample)








.. GENERATED FROM PYTHON SOURCE LINES 91-98

.. code-block:: Python

    graph = ot.Graph(
        "Initial LHS design of experiment - n=%d" % (sampleSize), "$x_0$", "$x_1$", True
    )
    cloud = ot.Cloud(inputSample)
    graph.add(cloud)
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_002.svg
   :alt: Initial LHS design of experiment - n=20
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_002.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 99-101

We now create the GP metamodel.
We selected the :class:`~openturns.MaternModel` covariance model with a constant basis as recommended in [leriche2021]_.

.. GENERATED FROM PYTHON SOURCE LINES 103-110

.. code-block:: Python

    covarianceModel = ot.MaternModel([1.0] * dim, [0.5], 2.5)
    basis = ot.ConstantBasisFactory(dim).build()
    fitter = otexp.GaussianProcessFitter(inputSample, outputSample, covarianceModel, basis)
    fitter.run()
    gpr = otexp.GaussianProcessRegression(fitter.getResult())
    gpr.run()








.. GENERATED FROM PYTHON SOURCE LINES 111-115

Create the optimization problem
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

We finally create the :class:`~openturns.OptimizationProblem` and solve it with class:`~openturns.EfficientGlobalOptimization`.

.. GENERATED FROM PYTHON SOURCE LINES 117-122

.. code-block:: Python

    problem = ot.OptimizationProblem()
    problem.setObjective(model)
    bounds = ot.Interval(lowerbound, upperbound)
    problem.setBounds(bounds)








.. GENERATED FROM PYTHON SOURCE LINES 123-127

In order to show the various options, we configure them all, even if most could be left to default settings in this case.

The most important method is :class:`~openturns.experimental.EfficientGlobalOptimization.setMaximumCallsNumber` which limits the number of iterations that the algorithm can perform.
In the Ackley example, we choose to perform 30 iterations of the algorithm.

.. GENERATED FROM PYTHON SOURCE LINES 129-134

.. code-block:: Python

    algo = otexp.EfficientGlobalOptimization(problem, gpr.getResult())
    algo.setMaximumCallsNumber(30)
    algo.run()
    result = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 135-137

.. code-block:: Python

    result.getIterationNumber()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    30



.. GENERATED FROM PYTHON SOURCE LINES 138-140

.. code-block:: Python

    result.getOptimalPoint()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=2 values=[-0.00627899,-0.00528803]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 141-143

.. code-block:: Python

    result.getOptimalValue()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=1 values=[0.0250124]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 144-146

.. code-block:: Python

    fexact






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=1 values=[4.44089e-16]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 147-153

Compared to the minimum function value, we see that the EGO algorithm
provides solution which is accurate.
Indeed, the optimal point is in the neighbourhood of the exact solution,
and this is quite an impressive success given the limited amount of function
evaluations: only 20 evaluations for the initial DOE and 30 iterations of
the EGO algorithm, for a total equal to 50 function evaluations.

.. GENERATED FROM PYTHON SOURCE LINES 155-160

.. code-block:: Python

    graph = result.drawOptimalValueHistory()
    optimum_curve = ot.Curve(ot.Sample([[0, fexact[0]], [29, fexact[0]]]))
    graph.add(optimum_curve)
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_003.svg
   :alt: Optimal value history
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_003.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 161-163

.. code-block:: Python

    inputHistory = result.getInputSample()








.. GENERATED FROM PYTHON SOURCE LINES 164-176

.. code-block:: Python

    graph = model.draw(lowerbound, upperbound, [100] * dim)
    graph.setTitle("Ackley function")
    cloud = ot.Cloud(inputSample, "initial")
    cloud.setPointStyle("bullet")
    cloud.setColor("black")
    graph.add(cloud)
    cloud = ot.Cloud(inputHistory, "solution")
    cloud.setPointStyle("diamond")
    cloud.setColor("forestgreen")
    graph.add(cloud)
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_004.svg
   :alt: Ackley function
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_004.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 177-179

We see that the initial (black) points are dispersed in the whole domain,
while the solution points are much closer to the solution.

.. GENERATED FROM PYTHON SOURCE LINES 181-185

Branin test-case
----------------

We now take a look at the :ref:`Branin-Hoo<use-case-branin>` function.

.. GENERATED FROM PYTHON SOURCE LINES 187-189

Define the problem
^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 191-192

The Branin model is defined in the usecases module in a data class `BraninModel` :

.. GENERATED FROM PYTHON SOURCE LINES 192-194

.. code-block:: Python

    bm = branin_function.BraninModel()








.. GENERATED FROM PYTHON SOURCE LINES 195-196

We load the dimension,

.. GENERATED FROM PYTHON SOURCE LINES 196-198

.. code-block:: Python

    dim = bm.dim








.. GENERATED FROM PYTHON SOURCE LINES 199-200

the domain boundaries,

.. GENERATED FROM PYTHON SOURCE LINES 200-203

.. code-block:: Python

    lowerbound = bm.lowerbound
    upperbound = bm.upperbound








.. GENERATED FROM PYTHON SOURCE LINES 204-205

and we load the model function and its noise:

.. GENERATED FROM PYTHON SOURCE LINES 205-208

.. code-block:: Python

    objectiveFunction = bm.model
    noise = bm.trueNoiseFunction








.. GENERATED FROM PYTHON SOURCE LINES 209-210

We build a sample out of the three minima :

.. GENERATED FROM PYTHON SOURCE LINES 210-212

.. code-block:: Python

    xexact = ot.Sample([bm.xexact1, bm.xexact2, bm.xexact3])








.. GENERATED FROM PYTHON SOURCE LINES 213-214

The minimum value attained `fexact` is :

.. GENERATED FROM PYTHON SOURCE LINES 214-217

.. code-block:: Python

    fexact = objectiveFunction(xexact)
    fexact






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr><td></td><th>y0</th></tr>
      <tr><th>0</th><td>-1.04741</td></tr>
      <tr><th>1</th><td>-1.04741</td></tr>
      <tr><th>2</th><td>-1.04741</td></tr>
    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 218-222

.. code-block:: Python

    graph = objectiveFunction.draw(lowerbound, upperbound, [100] * dim)
    graph.setTitle("Branin function")
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_005.svg
   :alt: Branin function
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_005.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 223-224

The Branin function has three local minimas.

.. GENERATED FROM PYTHON SOURCE LINES 226-228

Create the initial GP
^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 230-236

.. code-block:: Python

    distribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dim)
    sampleSize = 10 * dim
    experiment = ot.LHSExperiment(distribution, sampleSize)
    inputSample = experiment.generate()
    outputSample = objectiveFunction(inputSample)








.. GENERATED FROM PYTHON SOURCE LINES 237-244

.. code-block:: Python

    graph = ot.Graph(
        "Initial LHS design of experiment - n=%d" % (sampleSize), "$x_0$", "$x_1$", True
    )
    cloud = ot.Cloud(inputSample)
    graph.add(cloud)
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_006.svg
   :alt: Initial LHS design of experiment - n=20
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_006.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 245-246

Configure the covariance model with the noise

.. GENERATED FROM PYTHON SOURCE LINES 246-249

.. code-block:: Python

    covarianceModel = ot.MaternModel([1.0] * dim, [0.5], 2.5)
    covarianceModel.setNuggetFactor(noise)








.. GENERATED FROM PYTHON SOURCE LINES 250-251

Build the initial GP

.. GENERATED FROM PYTHON SOURCE LINES 251-257

.. code-block:: Python

    basis = ot.ConstantBasisFactory(dim).build()
    fitter = otexp.GaussianProcessFitter(inputSample, outputSample, covarianceModel, basis)
    fitter.run()
    gpr = otexp.GaussianProcessRegression(fitter.getResult())
    gpr.run()








.. GENERATED FROM PYTHON SOURCE LINES 258-260

Create and solve the problem
^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 262-263

We define the problem :

.. GENERATED FROM PYTHON SOURCE LINES 263-268

.. code-block:: Python

    problem = ot.OptimizationProblem()
    problem.setObjective(objectiveFunction)
    bounds = ot.Interval(lowerbound, upperbound)
    problem.setBounds(bounds)








.. GENERATED FROM PYTHON SOURCE LINES 269-271

We configure the algorithm
the nugget factor set in the covariance model with enable the AEI formulation

.. GENERATED FROM PYTHON SOURCE LINES 271-274

.. code-block:: Python

    algo = otexp.EfficientGlobalOptimization(problem, gpr.getResult())
    algo.setMaximumCallsNumber(30)








.. GENERATED FROM PYTHON SOURCE LINES 275-276

We run the algorithm and get the result:

.. GENERATED FROM PYTHON SOURCE LINES 276-279

.. code-block:: Python

    algo.run()
    result = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 280-282

.. code-block:: Python

    result.getIterationNumber()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    30



.. GENERATED FROM PYTHON SOURCE LINES 283-285

.. code-block:: Python

    result.getOptimalPoint()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=2 values=[0.541991,0.151992]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 286-288

.. code-block:: Python

    result.getOptimalValue()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=1 values=[-1.0474]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 289-291

.. code-block:: Python

    fexact






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr><td></td><th>y0</th></tr>
      <tr><th>0</th><td>-1.04741</td></tr>
      <tr><th>1</th><td>-1.04741</td></tr>
      <tr><th>2</th><td>-1.04741</td></tr>
    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 292-294

.. code-block:: Python

    inputHistory = result.getInputSample()








.. GENERATED FROM PYTHON SOURCE LINES 295-307

.. code-block:: Python

    graph = objectiveFunction.draw(lowerbound, upperbound, [100] * dim)
    graph.setTitle("Branin function")
    cloud = ot.Cloud(inputSample, "initial")
    cloud.setPointStyle("bullet")
    cloud.setColor("black")
    graph.add(cloud)
    cloud = ot.Cloud(inputHistory, "solution")
    cloud.setPointStyle("diamond")
    cloud.setColor("forestgreen")
    graph.add(cloud)
    view = viewer.View(graph)




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_007.svg
   :alt: Branin function
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_007.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 308-309

We see that the EGO algorithm reached the different optima locations.

.. GENERATED FROM PYTHON SOURCE LINES 311-317

.. code-block:: Python

    graph = result.drawOptimalValueHistory()
    optimum_curve = ot.Curve(ot.Sample([[0, fexact[0][0]], [29, fexact[0][0]]]))
    graph.add(optimum_curve)
    view = viewer.View(graph, axes_kw={"xticks": range(0, result.getIterationNumber(), 5)})





.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_008.svg
   :alt: Optimal value history
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_008.svg
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 318-319

.. code-block:: Python

    viewer.View.ShowAll()








.. _sphx_glr_download_auto_numerical_methods_optimization_plot_ego.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_ego.ipynb <plot_ego.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_ego.py <plot_ego.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_ego.zip <plot_ego.zip>`