.. 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

        Click :ref:`here <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 6-18

The EGO algorithm (Jones, 1998) is an adaptative optimization method based on
kriging.

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 experiment.

.. GENERATED FROM PYTHON SOURCE LINES 21-32

.. code-block:: default

    from openturns.usecases import branin_function as branin_function
    from openturns.usecases import ackley_function as ackley_function
    import openturns as ot
    import openturns.viewer as viewer
    from matplotlib import pylab as plt
    import math as m
    ot.RandomGenerator.SetSeed(0)
    ot.ResourceMap.SetAsString("KrigingAlgorithm-LinearAlgebra",  "LAPACK")
    ot.Log.Show(ot.Log.NONE)









.. GENERATED FROM PYTHON SOURCE LINES 33-37

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

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

.. GENERATED FROM PYTHON SOURCE LINES 39-41

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

.. GENERATED FROM PYTHON SOURCE LINES 43-44

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

.. GENERATED FROM PYTHON SOURCE LINES 44-46

.. code-block:: default

    am = ackley_function.AckleyModel()








.. GENERATED FROM PYTHON SOURCE LINES 47-48

We get the Ackley function :

.. GENERATED FROM PYTHON SOURCE LINES 48-50

.. code-block:: default

    model = am.model








.. GENERATED FROM PYTHON SOURCE LINES 51-52

We specify the domain of the model :

.. GENERATED FROM PYTHON SOURCE LINES 52-56

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 57-58

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 58-60

.. code-block:: default

    xexact = am.x0








.. GENERATED FROM PYTHON SOURCE LINES 61-62

The minimum value attained `fexact` is :

.. GENERATED FROM PYTHON SOURCE LINES 62-65

.. code-block:: default

    fexact = model(xexact)
    fexact






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[4.44089e-16]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 66-70

.. code-block:: default

    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.png
   :alt: Ackley function
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 71-72

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 74-78

Create the initial kriging
^^^^^^^^^^^^^^^^^^^^^^^^^^

Before using the EGO algorithm, we must create an initial kriging. In order to do this, we must create a design of experiment which fills the space. In this situation, the `LHSExperiment` is a good place to start (but other design of experiments may allow to better fill the space). We use a uniform distribution in order to create a LHS design with 50 points.

.. GENERATED FROM PYTHON SOURCE LINES 80-88

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 89-95

.. code-block:: default

    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.png
   :alt: Initial LHS design of experiment - n=50
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 96-97

We now create the kriging metamodel. We selected the `SquaredExponential` covariance model with a constant basis (the `MaternModel` may perform better in this case). We use default settings (1.0) for the scale parameters of the covariance model, but configure the amplitude to 0.1, which better corresponds to the properties of the Ackley function.

.. GENERATED FROM PYTHON SOURCE LINES 99-105

.. code-block:: default

    covarianceModel = ot.SquaredExponential([1.0] * dim, [0.5])
    basis = ot.ConstantBasisFactory(dim).build()
    kriging = ot.KrigingAlgorithm(
        inputSample, outputSample, covarianceModel, basis)
    kriging.run()








.. GENERATED FROM PYTHON SOURCE LINES 106-110

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

We finally create the `OptimizationProblem` and solve it with `EfficientGlobalOptimization`.

.. GENERATED FROM PYTHON SOURCE LINES 112-117

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 118-121

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 `setMaximumEvaluationNumber` which limits the number of iterations that the algorithm can perform. In the Ackley example, we choose to perform 10 iterations of the algorithm.

.. GENERATED FROM PYTHON SOURCE LINES 123-128

.. code-block:: default

    algo = ot.EfficientGlobalOptimization(problem, kriging.getResult())
    algo.setMaximumEvaluationNumber(10)
    algo.run()
    result = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 129-131

.. code-block:: default

    result.getIterationNumber()





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

 Out:

 .. code-block:: none


    10



.. GENERATED FROM PYTHON SOURCE LINES 132-134

.. code-block:: default

    result.getOptimalPoint()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[-0.0975796,0.839969]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 135-137

.. code-block:: default

    result.getOptimalValue()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[3.00508]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 138-140

.. code-block:: default

    fexact






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[4.44089e-16]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 141-142

Compared to the minimum function value, we see that the EGO algorithm provides solution which is not very accurate. However, 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 60 evaluations for the initial DOE and 10 iterations of the EGO algorithm, for a total equal to 70 function evaluations.

.. GENERATED FROM PYTHON SOURCE LINES 144-147

.. code-block:: default

    graph = result.drawOptimalValueHistory()
    view = viewer.View(graph)




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





.. GENERATED FROM PYTHON SOURCE LINES 148-150

.. code-block:: default

    inputHistory = result.getInputSample()








.. GENERATED FROM PYTHON SOURCE LINES 151-165

.. code-block:: default

    graph = model.draw(lowerbound, upperbound, [100]*dim)
    graph.setLegends([""])
    graph.setTitle(
        "Ackley function. Initial : black bullet. Solution : green diamond.")
    cloud = ot.Cloud(inputSample)
    cloud.setPointStyle("bullet")
    cloud.setColor("black")
    graph.add(cloud)
    cloud = ot.Cloud(inputHistory)
    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.png
   :alt: Ackley function. Initial : black bullet. Solution : green diamond.
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 166-169

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

However, the final solution produced by the EGO algorithm is not very accurate. This is why we finalize the process by adding a local optimization step.

.. GENERATED FROM PYTHON SOURCE LINES 171-176

.. code-block:: default

    algo2 = ot.NLopt(problem, 'LD_LBFGS')
    algo2.setStartingPoint(result.getOptimalPoint())
    algo2.run()
    result = algo2.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 177-179

.. code-block:: default

    result.getOptimalPoint()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[4.5981e-07,0.952166]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 180-181

The corrected solution is now extremely accurate.

.. GENERATED FROM PYTHON SOURCE LINES 183-186

.. code-block:: default

    graph = result.drawOptimalValueHistory()
    view = viewer.View(graph)




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





.. GENERATED FROM PYTHON SOURCE LINES 187-191

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

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

.. GENERATED FROM PYTHON SOURCE LINES 193-195

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

.. GENERATED FROM PYTHON SOURCE LINES 197-198

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

.. GENERATED FROM PYTHON SOURCE LINES 198-200

.. code-block:: default

    bm = branin_function.BraninModel()








.. GENERATED FROM PYTHON SOURCE LINES 201-202

We load the dimension,

.. GENERATED FROM PYTHON SOURCE LINES 202-204

.. code-block:: default

    dim = bm.dim








.. GENERATED FROM PYTHON SOURCE LINES 205-206

the domain boundaries,

.. GENERATED FROM PYTHON SOURCE LINES 206-209

.. code-block:: default

    lowerbound = bm.lowerbound
    upperbound = bm.upperbound








.. GENERATED FROM PYTHON SOURCE LINES 210-211

and we load the model function :

.. GENERATED FROM PYTHON SOURCE LINES 211-214

.. code-block:: default

    model = bm.model
    objectiveFunction = model.getMarginal(0)








.. GENERATED FROM PYTHON SOURCE LINES 215-216

We build a sample out of the three minima :

.. GENERATED FROM PYTHON SOURCE LINES 216-218

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 219-220

The minimum value attained `fexact` is :

.. GENERATED FROM PYTHON SOURCE LINES 220-223

.. code-block:: default

    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 224-228

.. code-block:: default

    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_006.png
   :alt: Branin function
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 229-230

The Branin function has three local minimas.

.. GENERATED FROM PYTHON SOURCE LINES 232-234

Create the initial kriging
^^^^^^^^^^^^^^^^^^^^^^^^^^

.. GENERATED FROM PYTHON SOURCE LINES 236-243

.. code-block:: default

    distribution = ot.ComposedDistribution([ot.Uniform(0.0, 1.0)] * dim)
    sampleSize = 50
    experiment = ot.LHSExperiment(distribution, sampleSize)
    inputSample = experiment.generate()
    modelEval = model(inputSample)
    outputSample = modelEval.getMarginal(0)








.. GENERATED FROM PYTHON SOURCE LINES 244-250

.. code-block:: default

    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_007.png
   :alt: Initial LHS design of experiment - n=50
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_007.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 251-256

.. code-block:: default

    covarianceModel = ot.SquaredExponential([1.0] * dim, [1.0])
    basis = ot.ConstantBasisFactory(dim).build()
    kriging = ot.KrigingAlgorithm(
        inputSample, outputSample, covarianceModel, basis)








.. GENERATED FROM PYTHON SOURCE LINES 257-261

.. code-block:: default

    noise = [x[1] for x in modelEval]
    kriging.setNoise(noise)
    kriging.run()








.. GENERATED FROM PYTHON SOURCE LINES 262-264

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

.. GENERATED FROM PYTHON SOURCE LINES 266-267

We define the problem :

.. GENERATED FROM PYTHON SOURCE LINES 267-272

.. code-block:: default

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








.. GENERATED FROM PYTHON SOURCE LINES 273-274

We configure the maximum number of function evaluations to 20. We assume that the function is noisy, with a constant variance.

.. GENERATED FROM PYTHON SOURCE LINES 276-277

We configure the algorithm :

.. GENERATED FROM PYTHON SOURCE LINES 277-286

.. code-block:: default

    algo = ot.EfficientGlobalOptimization(problem, kriging.getResult())
    # assume constant noise var
    guessedNoiseFunction = 0.1
    noiseModel = ot.SymbolicFunction(['x1', 'x2'], [str(guessedNoiseFunction)])
    algo.setNoiseModel(noiseModel)
    algo.setMaximumEvaluationNumber(20)
    algo.run()
    result = algo.getResult()








.. GENERATED FROM PYTHON SOURCE LINES 287-289

.. code-block:: default

    result.getIterationNumber()





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

 Out:

 .. code-block:: none


    20



.. GENERATED FROM PYTHON SOURCE LINES 290-292

.. code-block:: default

    result.getOptimalPoint()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[0.124453,0.803571]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 293-295

.. code-block:: default

    result.getOptimalValue()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <p>[-1.04662]</p>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 296-298

.. code-block:: default

    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 299-301

.. code-block:: default

    inputHistory = result.getInputSample()








.. GENERATED FROM PYTHON SOURCE LINES 302-316

.. code-block:: default

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




.. image-sg:: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_008.png
   :alt: Branin function. Initial : black bullet. Solution : green diamond.
   :srcset: /auto_numerical_methods/optimization/images/sphx_glr_plot_ego_008.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 317-318

We see that the EGO algorithm found the second local minimum. Given the limited number of function evaluations, the other local minimas have not been explored.

.. GENERATED FROM PYTHON SOURCE LINES 320-325

.. code-block:: default

    graph = result.drawOptimalValueHistory()
    view = viewer.View(graph, axes_kw={"xticks": range(
        0, result.getIterationNumber(), 5)})

    plt.show()



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






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  5.802 seconds)


.. _sphx_glr_download_auto_numerical_methods_optimization_plot_ego.py:


.. only :: html

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



  .. 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-jupyter

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


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_