Repository: OpenSourceEconomics/estimagic

Issue:
Add pygmo optimizers
### Current situation

No pygmo optimizers are supported 

### Desired Situation

- Support all pygmo optimizers that are not clearly dominated by other optimizers. 
- Improve name of options where necessary, but stay close to original names.


Discussion:
review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/OpenSourceEconomics/estimagic/pull/215"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>
codecov[bot]: # [Codecov](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics) Report
> Merging [#215](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics) (7c9353c) into [master](https://codecov.io/gh/OpenSourceEconomics/estimagic/commit/238c5e26208b59d1e0caa30a3d35d99f2241ac2a?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics) (238c5e2) will **not change** coverage.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215/graphs/tree.svg?width=650&height=150&src=pr&token=INbcX5rQK2&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics)](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics)

```diff
@@           Coverage Diff           @@
##           master     #215   +/-   ##
=======================================
  Coverage   83.54%   83.54%           
=======================================
  Files          58       58           
  Lines        3992     3992           
=======================================
  Hits         3335     3335           
  Misses        657      657           
```



------

[Continue to review full report at Codecov](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics). Last update [238c5e2...7c9353c](https://codecov.io/gh/OpenSourceEconomics/estimagic/pull/215?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=OpenSourceEconomics).


Patch:
diff --git a/.codecov.yml b/.codecov.yml
index eab9118..ddcf264 100644
--- a/.codecov.yml
+++ b/.codecov.yml
@@ -5,4 +5,4 @@ coverage:
         target: 80%
     project:
       default:
-        target: 84%
+        target: 87%
diff --git a/docs/source/getting_started/installation.rst b/docs/source/getting_started/installation.rst
index 887353c..0662c44 100644
--- a/docs/source/getting_started/installation.rst
+++ b/docs/source/getting_started/installation.rst
@@ -43,4 +43,4 @@ To enable all algorithms at once, do the following:
 
 ``conda install cyipopt``
 
-.. ``conda install pygmo``
+``conda install pygmo``
diff --git a/docs/source/how_to_guides/optimization/how_to_specify_algorithm_and_algo_options.rst b/docs/source/how_to_guides/optimization/how_to_specify_algorithm_and_algo_options.rst
index 4dfe539..d523996 100644
--- a/docs/source/how_to_guides/optimization/how_to_specify_algorithm_and_algo_options.rst
+++ b/docs/source/how_to_guides/optimization/how_to_specify_algorithm_and_algo_options.rst
@@ -188,6 +188,85 @@ Py-BOBYQA``).
 
     .. autofunction:: estimagic.optimization.nag_optimizers.nag_pybobyqa
 
+
+
+.. _pygmo_algorithms:
+
+PYGMO2 Optimizers
+------------------
+
+Please cite :cite:`Biscani2020` in addition to estimagic when using pygmo.
+estimagic supports the following `pygmo2 <https://esa.github.io/pygmo2>`_
+optimizers.
+
+.. dropdown::  pygmo_gaco
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_gaco
+
+.. dropdown::  pygmo_bee_colony
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_bee_colony
+
+.. dropdown::  pygmo_de
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_de
+
+.. dropdown::  pygmo_sea
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_sea
+
+.. dropdown::  pygmo_sga
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_sga
+
+.. dropdown::  pygmo_sade
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_sade
+
+
+.. dropdown::  pygmo_cmaes
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_cmaes
+
+.. dropdown::  pygmo_simulated_annealing
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_simulated_annealing
+
+.. dropdown::  pygmo_pso
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_pso
+
+.. dropdown::  pygmo_pso_gen
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_pso_gen
+
+.. dropdown::  pygmo_mbh
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_mbh
+
+.. dropdown::  pygmo_xnes
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_xnes
+
+.. dropdown::  pygmo_gwo
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_gwo
+
+.. dropdown::  pygmo_compass_search
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_compass_search
+
+.. dropdown::  pygmo_ihs
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_ihs
+
+.. dropdown::  pygmo_de1220
+
+    .. autofunction:: estimagic.optimization.pygmo_optimizers.pygmo_de1220
+
+
+.. _ipopt_algorithm:
+
 The Interior Point Optimizer (ipopt)
 ------------------------------------
 
@@ -282,10 +361,6 @@ using an NLOPT algorithm. To install nlopt run ``conda install nlopt``.
 **References**
 
 .. bibliography:: ../../refs.bib
+    :labelprefix: algo_
     :filter: docname in docnames
     :style: unsrt
-
-
-
-.. With ``nlopt`` installed
-.. With ``pygmo`` installed
diff --git a/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst b/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst
index 7ec414c..bb61e2d 100644
--- a/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst
+++ b/docs/source/how_to_guides/optimization/how_to_specify_constraints.rst
@@ -260,8 +260,6 @@ that your read that explanation first!
     You can specify either value or lower and upper bounds. Suppose you have the
     following params DataFrame:
 
-
-
     +-------------------+-------+
     |                   | value |
     +----------+--------+-------+
diff --git a/docs/source/refs.bib b/docs/source/refs.bib
index a56d2f8..3b851e8 100644
--- a/docs/source/refs.bib
+++ b/docs/source/refs.bib
@@ -1,6 +1,7 @@
 % Encoding: UTF-8
 
 
+
 @Book{Dennis1996,
   Title                    = {Numerical Methods for Unconstrained Optimization and Nonlinear Equations},
   Author                   = {Dennis, J.E. and Schnabel, R.B.},
@@ -13,96 +14,91 @@
   Url                      = {https://books.google.de/books?id=RtxcWd0eBD0C&redir_esc=y}
 }
 
-@Book{Hansen2019,
-  Title                    = {Econometrics},
-  Author                   = {Bruce E. Hansen},
-  Editor                   = {Bruce E. Hansen},
-  Publisher                = {Unpublished},
-  Year                     = {2019},
-
-  Url                      = {https://www.ssc.wisc.edu/~bhansen/econometrics/},
-
-  Owner                    = {janos},
-  Timestamp                = {2019.10.03}
+@book{Hansen2019,
+  title     = {Econometrics},
+  author    = {Bruce E. Hansen},
+  editor    = {Bruce E. Hansen},
+  publisher = {Unpublished},
+  year      = {2019},
+  url       = {https://www.ssc.wisc.edu/~bhansen/econometrics/},
+  owner     = {janos},
+  timestamp = {2019.10.03}
 }
 
-@Book{Hansen2020,
-  Title                    = {Econometrics},
-  Author                   = {Bruce E. Hansen},
-  Editor                   = {Bruce E. Hansen},
-  Publisher                = {Unpublished},
-  Year                     = {2020},
-
-  Address                  = {https://www.ssc.wisc.edu/~bhansen/econometrics/},
-
-  Owner                    = {janos},
-  Timestamp                = {2020.03.04}
+@book{Hansen2020,
+  title     = {Econometrics},
+  author    = {Bruce E. Hansen},
+  editor    = {Bruce E. Hansen},
+  publisher = {Unpublished},
+  year      = {2020},
+  address   = {https://www.ssc.wisc.edu/~bhansen/econometrics/},
+  owner     = {janos},
+  timestamp = {2020.03.04}
 }
 
 
-@Book{Verbeek2008,
-  Title                    = {A Guide to Modern Econometrics},
-  Author                   = {Verbeek, M.},
-  Publisher                = {Wiley},
-  Year                     = {2008},
-
-  ISBN                     = {9780470517697},
-  Lccn                     = {2007050167},
-  Url                      = {https://books.google.com/books?id=uEFm6pAJZhoC}
+@book{Verbeek2008,
+  title     = {A Guide to Modern Econometrics},
+  author    = {Verbeek, M.},
+  publisher = {Wiley},
+  year      = {2008},
+  isbn      = {9780470517697},
+  lccn      = {2007050167},
+  url       = {https://books.google.com/books?id=uEFm6pAJZhoC}
 }
 
 @book{Wassermann2006,
-  title={All of nonparametric statistics},
-  author={Wasserman, Larry},
-  year={2006},
-  publisher={Springer Science \& Business Media}
+  title     = {All of nonparametric statistics},
+  author    = {Wasserman, Larry},
+  year      = {2006},
+  publisher = {Springer Science \& Business Media}
 }
 
-@Article{Groeneveld1994,
+@article{Groeneveld1994,
   author  = {Eildert Groeneveld},
   title   = {A reparameterization to improve numerical optimization in multivariate REML (co)variance component estimation},
   journal = {Genetics, Selection, Evolution : GSE},
   year    = {1994},
   volume  = {26},
-  pages   = {537 - 545},
+  pages   = {537 - 545}
 }
 
-@Article{Pinheiro1996,
+@article{Pinheiro1996,
   author  = {José C. Pinheiro and Douglas M. Bates},
   title   = {Unconstrained Parameterizations for Variance-Covariance Matrices},
   journal = {Statistics and Computing},
   year    = {1996},
   volume  = {6},
-  pages   = {289--296},
+  pages   = {289--296}
 }
 
-@TechReport{Kraft1988,
+@techreport{Kraft1988,
   author      = {Kraft, Dieter},
   institution = {DLR German Aerospace Center – Institute for Flight Mechanics},
   title       = {A software package for sequential quadratic programming},
   year        = {1988},
   address     = {Köln, Germany},
-  url         = {http://degenerateconic.com/wp-content/uploads/2018/03/DFVLR_FB_88_28.pdf},
+  url         = {http://degenerateconic.com/wp-content/uploads/2018/03/DFVLR_FB_88_28.pdf}
 }
 
-@Book{Nocedal2006,
+@book{Nocedal2006,
   author    = {Nocedal, Jorge and Wright, Stephen},
   publisher = {Springer Science \& Business Media},
   title     = {Numerical optimization},
-  year      = {2006},
+  year      = {2006}
 }
 
-@InCollection{Conn2000,
+@incollection{Conn2000,
   author    = {Conn, AR and Gould, NI and Toint, PL},
   booktitle = {Trust region methods},
   publisher = {Siam},
   title     = {Nonlinear equations and nonlinear fitting},
   year      = {2000},
   pages     = {749--774},
-  volume    = {1},
+  volume    = {1}
 }
 
-@Article{Byrd1999,
+@article{Byrd1999,
   author    = {Byrd, Richard H and Hribar, Mary E and Nocedal, Jorge},
   journal   = {SIAM Journal on Optimization},
   title     = {An interior point algorithm for large-scale nonlinear programming},
@@ -110,10 +106,10 @@
   number    = {4},
   pages     = {877--900},
   volume    = {9},
-  publisher = {SIAM},
+  publisher = {SIAM}
 }
 
-@Article{Lalee1998,
+@article{Lalee1998,
   author    = {Lalee, Marucha and Nocedal, Jorge and Plantenga, Todd},
   journal   = {SIAM Journal on Optimization},
   title     = {On the implementation of an algorithm for large-scale equality constrained optimization},
@@ -121,10 +117,10 @@
   number    = {3},
   pages     = {682--706},
   volume    = {8},
-  publisher = {SIAM},
+  publisher = {SIAM}
 }
 
-@Article{Gao2012,
+@article{Gao2012,
   author    = {Gao, Fuchang and Han, Lixing},
   journal   = {Computational Optimization and Applications},
   title     = {Implementing the Nelder-Mead simplex algorithm with adaptive parameters},
@@ -132,28 +128,28 @@
   number    = {1},
   pages     = {259--277},
   volume    = {51},
-  publisher = {Springer},
+  publisher = {Springer}
 }
 
-@InCollection{Powell1994,
+@incollection{Powell1994,
   author    = {Powell, Michael JD},
   booktitle = {Advances in optimization and numerical analysis},
   publisher = {Springer},
   title     = {A direct search optimization method that models the objective and constraint functions by linear interpolation},
   year      = {1994},
-  pages     = {51--67},
+  pages     = {51--67}
 }
 
-@Article{Powell1998,
+@article{Powell1998,
   author    = {Powell, Michael JD},
   journal   = {Acta numerica},
   title     = {Direct search algorithms for optimization calculations},
   year      = {1998},
   pages     = {287--336},
-  publisher = {Cambridge University Press},
+  publisher = {Cambridge University Press}
 }
 
-@Article{Powell2007,
+@article{Powell2007,
   author    = {Powell, Michael JD},
   journal   = {Mathematics Today-Bulletin of the Institute of Mathematics and its Applications},
   title     = {A view of algorithms for optimization without derivatives},
@@ -161,80 +157,80 @@
   number    = {5},
   pages     = {170--174},
   volume    = {43},
-  publisher = {Citeseer},
+  publisher = {Citeseer}
 }
 
-@TechReport{Benson2017,
+@techreport{Benson2017,
   author      = {Benson, S and McInnes, LC and Mor{\'e}, JJ and Munson, T and Sarich, J},
   institution = {Technical Report ANL/MCS-TM-322, Argonne National Laboratory},
   title       = {TAO user manual (revision 3.7)},
   year        = {2017},
-  url         = {http://web.mit.edu/tao-petsc_v3.7/tao_manual.pdf},
+  url         = {http://web.mit.edu/tao-petsc_v3.7/tao_manual.pdf}
 }
 
-@TechReport{Wild2015,
+@techreport{Wild2015,
   author      = {Wild, Stefan M.},
   institution = {Argonne National Laboratory},
   title       = {Solving Derivative-Free Nonlinear Least Squares Problems with POUNDERS},
   year        = {2015},
-  url         = {https://www.mcs.anl.gov/papers/P5120-0414.pdf},
+  url         = {https://www.mcs.anl.gov/papers/P5120-0414.pdf}
 }
 
-@Misc{Cartis2018,
+@misc{Cartis2018,
   author        = {Coralia Cartis and Jan Fiala and Benjamin Marteau and Lindon Roberts},
   title         = {Improving the Flexibility and Robustness of Model-Based Derivative-Free Optimization Solvers},
   year          = {2018},
   archiveprefix = {arXiv},
   eprint        = {1804.00154},
-  primaryclass  = {math.OC},
+  primaryclass  = {math.OC}
 }
 
-@Misc{Cartis2018a,
+@misc{Cartis2018a,
   author        = {Coralia Cartis and Lindon Roberts and Oliver Sheridan-Methven},
   title         = {Escaping local minima with derivative-free methods: a numerical investigation},
   year          = {2018},
   archiveprefix = {arXiv},
   eprint        = {1812.11343},
-  primaryclass  = {math.OC},
+  primaryclass  = {math.OC}
 }
 
-@Article{Powell2009,
+@article{Powell2009,
   author  = {Powell, Michael JD},
   journal = {Cambridge NA Report NA2009/06, University of Cambridge, Cambridge},
   title   = {The BOBYQA algorithm for bound constrained optimization without derivatives},
   year    = {2009},
-  pages   = {26--46},
+  pages   = {26--46}
 }
 
-@Misc{Cartis2018b,
+@misc{Cartis2018b,
   author        = {Coralia Cartis and Jan Fiala and Benjamin Marteau and Lindon Roberts},
   title         = {Improving the Flexibility and Robustness of Model-Based Derivative-Free Optimization Solvers},
   year          = {2018},
   archiveprefix = {arXiv},
   eprint        = {1804.00154},
-  primaryclass  = {math.OC},
+  primaryclass  = {math.OC}
 }
 
-@Unpublished{Saxton2018,
+@unpublished{Saxton2018,
   title  = {resample: Randomization-based inference in Python.},
   author = {Daniel Saxton},
   year   = {2018},
-  note = {unpublished},
-  Url    = {https://github.com/dsaxton/resample},
+  note   = {unpublished},
+  url    = {https://github.com/dsaxton/resample}
 }
 
 @article{CameronMiller2015,
-  title={A practitioner’s guide to cluster-robust inference},
-  author={Cameron, A Colin and Miller, Douglas L},
-  journal={Journal of human resources},
-  volume={50},
-  number={2},
-  pages={317--372},
-  year={2015},
-  publisher={University of Wisconsin Press}
+  title     = {A practitioner’s guide to cluster-robust inference},
+  author    = {Cameron, A Colin and Miller, Douglas L},
+  journal   = {Journal of human resources},
+  volume    = {50},
+  number    = {2},
+  pages     = {317--372},
+  year      = {2015},
+  publisher = {University of Wisconsin Press}
 }
 
-@Article{Waechter2005,
+@article{Waechter2005,
   author    = {Andreas Wächter and Lorenz T. Biegler},
   journal   = {{SIAM} Journal on Optimization},
   title     = {Line Search Filter Methods for Nonlinear Programming: Local Convergence},
@@ -244,10 +240,10 @@
   pages     = {32--48},
   volume    = {16},
   doi       = {10.1137/s1052623403426544},
-  publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})},
+  publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})}
 }
 
-@Article{Waechter2005a,
+@article{Waechter2005a,
   author    = {Andreas Wächter and Lorenz T. Biegler},
   journal   = {{SIAM} Journal on Optimization},
   title     = {Line Search Filter Methods for Nonlinear Programming: Motivation and Global Convergence},
@@ -257,10 +253,10 @@
   pages     = {1--31},
   volume    = {16},
   doi       = {10.1137/s1052623403426556},
-  publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})},
+  publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})}
 }
 
-@Article{Waechter2005b,
+@article{Waechter2005b,
   author    = {Andreas Wächter and Lorenz T. Biegler},
   journal   = {Mathematical Programming},
   title     = {On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming},
@@ -270,10 +266,10 @@
   pages     = {25--57},
   volume    = {106},
   doi       = {10.1007/s10107-004-0559-y},
-  publisher = {Springer Science and Business Media {LLC}},
+  publisher = {Springer Science and Business Media {LLC}}
 }
 
-@Article{Nocedal2009,
+@article{Nocedal2009,
   author    = {Jorge Nocedal and Andreas Wächter and Richard A. Waltz},
   journal   = {{SIAM} Journal on Optimization},
   title     = {Adaptive Barrier Update Strategies for Nonlinear Interior Methods},
@@ -283,7 +279,218 @@
   pages     = {1674--1693},
   volume    = {19},
   doi       = {10.1137/060649513},
-  publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})},
+  publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})}
+}
+
+@article{Schlueter2009,
+  author    = {Martin Schlüter and Jose A. Egea and Julio R. Banga},
+  journal   = {Computers {\&} Operations Research},
+  title     = {Extended ant colony optimization for non-convex mixed integer nonlinear programming},
+  year      = {2009},
+  month     = {jul},
+  number    = {7},
+  pages     = {2217--2229},
+  volume    = {36},
+  doi       = {10.1016/j.cor.2008.08.015},
+  publisher = {Elsevier {BV}}
+}
+
+@article{Karaboga2007,
+  author    = {Dervis Karaboga and Bahriye Basturk},
+  journal   = {Journal of Global Optimization},
+  title     = {A powerful and efficient algorithm for numerical function optimization: artificial bee colony ({ABC}) algorithm},
+  year      = {2007},
+  month     = {apr},
+  number    = {3},
+  pages     = {459--471},
+  volume    = {39},
+  doi       = {10.1007/s10898-007-9149-x},
+  publisher = {Springer Science and Business Media {LLC}}
+}
+
+@article{Mernik2015,
+  author    = {Marjan Mernik and Shih-Hsi Liu and Dervis Karaboga and Matej {\v{C}}repin{\v{s}}ek},
+  journal   = {Information Sciences},
+  title     = {On clarifying misconceptions when comparing variants of the Artificial Bee Colony Algorithm by offering a new implementation},
+  year      = {2015},
+  month     = {jan},
+  pages     = {115--127},
+  volume    = {291},
+  doi       = {10.1016/j.ins.2014.08.040},
+  publisher = {Elsevier {BV}}
+}
+
+@article{Storn1997,
+  author    = {Rainer Storn and Kenneth Price},
+  journal   = {Journal of Global Optimization},
+  title     = {Differential Evolution – A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces},
+  year      = {1997},
+  number    = {4},
+  pages     = {341--359},
+  volume    = {11},
+  doi       = {10.1023/a:1008202821328},
+  publisher = {Springer Science and Business Media {LLC}}
+}
+
+@article{Oliveto2007,
+  author    = {Pietro S. Oliveto and Jun He and Xin Yao},
+  journal   = {International Journal of Automation and Computing},
+  title     = {Time complexity of evolutionary algorithms for combinatorial optimization: A decade of results},
+  year      = {2007},
+  month     = {jul},
+  number    = {3},
+  pages     = {281--293},
+  volume    = {4},
+  doi       = {10.1007/s11633-007-0281-3},
+  publisher = {Springer Science and Business Media {LLC}}
+}
+
+
+@article{Brest2006,
+  author  = {Brest, Janez and Greiner, Sao and Boskovic, Borko and Mernik, Marjan and Zumer, Viljem},
+  journal = {IEEE Transactions on Evolutionary Computation},
+  title   = {Self-Adapting Control Parameters in Differential Evolution: A Comparative Study on Numerical Benchmark Problems},
+  year    = {2006},
+  number  = {6},
+  pages   = {646-657},
+  volume  = {10},
+  doi     = {10.1109/TEVC.2006.872133}
+}
+
+@inproceedings{Elsayed2011,
+  author    = {Elsayed, Saber M. and Sarker, Ruhul A. and Essam, Daryl L.},
+  booktitle = {2011 IEEE Congress of Evolutionary Computation (CEC)},
+  title     = {Differential evolution with multiple strategies for solving CEC2011 real-world numerical optimization problems},
+  year      = {2011},
+  pages     = {1041-1048},
+  doi       = {10.1109/CEC.2011.5949732}
+}
+
+@incollection{Hansen2006,
+  author    = {Nikolaus Hansen},
+  booktitle = {Towards a New Evolutionary Computation},
+  publisher = {Springer Berlin Heidelberg},
+  title     = {The {CMA} Evolution Strategy: A Comparing Review},
+  year      = {2006},
+  pages     = {75--102},
+  doi       = {10.1007/3-540-32494-1_4}
+}
+
+@article{Corana1987,
+  author     = {Corana, A. and Marchesi, M. and Martini, C. and Ridella, S.},
+  title      = {Minimizing Multimodal Functions of Continuous Variables with the “Simulated Annealing” Algorithm—Corrigenda for This Article is Available Here},
+  year       = {1987},
+  issue_date = {Sept. 1987},
+  publisher  = {Association for Computing Machinery},
+  address    = {New York, NY, USA},
+  volume     = {13},
+  number     = {3},
+  issn       = {0098-3500},
+  url        = {https://doi.org/10.1145/29380.29864},
+  doi        = {10.1145/29380.29864},
+  abstract   = {A new global optimization algorithm for functions of continuous variables is presented,
+derived from the “Simulated Annealing” algorithm recently introduced in combinatorial
+optimization.The algorithm is essentially an iterative random search procedure with
+adaptive moves along the coordinate directions. It permits uphill moves under the
+control of a probabilistic criterion, thus tending to avoid the first local minima
+encountered.The algorithm has been tested against the Nelder and Mead simplex method
+and against a version of Adaptive Random Search. The test functions were Rosenbrock
+valleys and multiminima functions in 2,4, and 10 dimensions.The new method proved
+to be more reliable than the others, being always able to find the optimum, or at
+least a point very close to it. It is quite costly in term of function evaluations,
+but its cost can be predicted in advance, depending only slightly on the starting
+point.},
+  journal    = {ACM Trans. Math. Softw.},
+  month      = sep,
+  pages      = {262–280},
+  numpages   = {19}
+}
+
+@Article{Poli2007,
+  author    = {Riccardo Poli and James Kennedy and Tim Blackwell},
+  journal   = {Swarm Intelligence},
+  title     = {Particle swarm optimization},
+  year      = {2007},
+  month     = {aug},
+  number    = {1},
+  pages     = {33--57},
+  volume    = {1},
+  doi       = {10.1007/s11721-007-0002-0},
+  publisher = {Springer Science and Business Media {LLC}},
+}
+
+@Article{Wales1997,
+  author    = {David J. Wales and Jonathan P. K. Doye},
+  journal   = {The Journal of Physical Chemistry A},
+  title     = {Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms},
+  year      = {1997},
+  month     = {jul},
+  number    = {28},
+  pages     = {5111--5116},
+  volume    = {101},
+  doi       = {10.1021/jp970984n},
+  publisher = {American Chemical Society ({ACS})},
 }
 
+@InProceedings{Glasmachers2010,
+  author    = {Tobias Glasmachers and Tom Schaul and Sun Yi and Daan Wierstra and Jürgen Schmidhuber},
+  booktitle = {Proceedings of the 12th annual conference on Genetic and evolutionary computation - {GECCO} {\textquotesingle}10},
+  title     = {Exponential natural evolution strategies},
+  year      = {2010},
+  publisher = {{ACM} Press},
+  doi       = {10.1145/1830483.1830557},
+}
+
+@Article{Mahdavi2007,
+  author    = {M. Mahdavi and M. Fesanghary and E. Damangir},
+  journal   = {Applied Mathematics and Computation},
+  title     = {An improved harmony search algorithm for solving optimization problems},
+  year      = {2007},
+  month     = {may},
+  number    = {2},
+  pages     = {1567--1579},
+  volume    = {188},
+  doi       = {10.1016/j.amc.2006.11.033},
+  publisher = {Elsevier {BV}},
+}
+
+@Article{Mirjalili2014,
+  author    = {Seyedali Mirjalili and Seyed Mohammad Mirjalili and Andrew Lewis},
+  journal   = {Advances in Engineering Software},
+  title     = {Grey Wolf Optimizer},
+  year      = {2014},
+  month     = {mar},
+  pages     = {46--61},
+  volume    = {69},
+  doi       = {10.1016/j.advengsoft.2013.12.007},
+  publisher = {Elsevier {BV}},
+}
+
+@Article{Kolda2003,
+  author    = {Tamara G. Kolda and Robert Michael Lewis and Virginia Torczon},
+  journal   = {{SIAM} Review},
+  title     = {Optimization by Direct Search: New Perspectives on Some Classical and Modern Methods},
+  year      = {2003},
+  month     = {jan},
+  number    = {3},
+  pages     = {385--482},
+  volume    = {45},
+  doi       = {10.1137/s003614450242889},
+}
+
+@Article{Biscani2020,
+  doi = {10.21105/joss.02338},
+  url = {https://doi.org/10.21105/joss.02338},
+  year = {2020},
+  publisher = {The Open Journal},
+  volume = {5},
+  number = {53},
+  pages = {2338},
+  author = {Francesco Biscani and Dario Izzo},
+  title = {A parallel global multiobjective framework for optimization: pagmo},
+  journal = {Journal of Open Source Software}
+}
+
+
+
 @Comment{jabref-meta: databaseType:bibtex;}
diff --git a/environment.yml b/environment.yml
index dcdb6fc..6cb3575 100644
--- a/environment.yml
+++ b/environment.yml
@@ -6,7 +6,6 @@ channels:
   - opensourceeconomics
 
 dependencies:
-
   - python=3.8
   - pip
   - anaconda-client
@@ -15,6 +14,7 @@ dependencies:
   - click
   - conda-build
   - conda-verify
+  - cyipopt
   - fuzzywuzzy
   - joblib
   - cloudpickle
@@ -24,6 +24,7 @@ dependencies:
   - pandas
   - pdbpp
   - petsc4py>=3.11
+  - pygmo
   - pytest
   - pytest-cov
   - pytest-xdist
@@ -39,15 +40,14 @@ dependencies:
   - pytask>=0.0.11
   - nlopt
   - sphinx-panels
-
-  - cyipopt
+  - pygmo
 
   - pip:
-    - black
-    - blackcellmagic
-    - bump2version
-    - pre-commit
-    - Py-BOBYQA
-    - DFO-LS
-    - gif
-    - gif[matplotlib]
+      - black
+      - blackcellmagic
+      - bump2version
+      - pre-commit
+      - Py-BOBYQA
+      - DFO-LS
+      - gif
+      - gif[matplotlib]
diff --git a/estimagic/config.py b/estimagic/config.py
index 314bd81..8e2fe55 100644
--- a/estimagic/config.py
+++ b/estimagic/config.py
@@ -53,6 +53,13 @@ except ImportError:
 else:
     IS_DFOLS_INSTALLED = True
 
+try:
+    import pygmo  # noqa: F401
+except ImportError:
+    IS_PYGMO_INSTALLED = False
+else:
+    IS_PYGMO_INSTALLED = True
+
 
 try:
     import cyipopt  # noqa: F401
diff --git a/estimagic/optimization/__init__.py b/estimagic/optimization/__init__.py
index 71541b2..ec75616 100644
--- a/estimagic/optimization/__init__.py
+++ b/estimagic/optimization/__init__.py
@@ -5,9 +5,11 @@ from estimagic.config import IS_DFOLS_INSTALLED
 from estimagic.config import IS_NLOPT_INSTALLED
 from estimagic.config import IS_PETSC4PY_INSTALLED
 from estimagic.config import IS_PYBOBYQA_INSTALLED
+from estimagic.config import IS_PYGMO_INSTALLED
 from estimagic.optimization import cyipopt_optimizers
 from estimagic.optimization import nag_optimizers
 from estimagic.optimization import nlopt_optimizers
+from estimagic.optimization import pygmo_optimizers
 from estimagic.optimization import scipy_optimizers
 from estimagic.optimization import tao_optimizers
 
@@ -25,6 +27,7 @@ if IS_NLOPT_INSTALLED:
     COLLECTED_FUNCTIONS.update(
         **dict(inspect.getmembers(nlopt_optimizers, inspect.isfunction))
     )
+
 # drop private and helper functions
 AVAILABLE_ALGORITHMS = {}
 PUBLIC_HELPERS = [
@@ -42,7 +45,21 @@ if IS_PYBOBYQA_INSTALLED:
 if IS_DFOLS_INSTALLED:
     AVAILABLE_ALGORITHMS["nag_dfols"] = nag_optimizers.nag_dfols
 
+if IS_PYGMO_INSTALLED:
+    _PYGMO_FUNCTIONS = dict(inspect.getmembers(pygmo_optimizers, inspect.isfunction))
+    PYGMO_ALGORITHMS = {
+        k: v for k, v in _PYGMO_FUNCTIONS.items() if k.startswith("pygmo_")
+    }
+    AVAILABLE_ALGORITHMS.update(**PYGMO_ALGORITHMS)
+
 if IS_CYIPOPT_INSTALLED:
     AVAILABLE_ALGORITHMS["ipopt"] = cyipopt_optimizers.ipopt
 
-GLOBAL_ALGORITHMS = ["nlopt_direct", "nlopt_esch", "nlopt_isres", "nlopt_crs2_lm"]
+GLOBAL_ALGORITHMS = [
+    "nlopt_direct",
+    "nlopt_esch",
+    "nlopt_isres",
+    "nlopt_crs2_lm",
+]
+if IS_PYGMO_INSTALLED:
+    GLOBAL_ALGORITHMS += PYGMO_ALGORITHMS.keys()
diff --git a/estimagic/optimization/algo_options.py b/estimagic/optimization/algo_options.py
index 6ff1ddc..c38639d 100644
--- a/estimagic/optimization/algo_options.py
+++ b/estimagic/optimization/algo_options.py
@@ -101,7 +101,7 @@ STOPPING_MAX_CRITERION_EVALUATIONS = 1_000_000
 """
 
 
-STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL = 5_000
+STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL = 1_000
 """int:
     If the maximum number of function evaluation is reached, the optimization stops
     but we do not count this as successful convergence. The function evaluations used
diff --git a/estimagic/optimization/nlopt_optimizers.py b/estimagic/optimization/nlopt_optimizers.py
index b6981c5..f7e639e 100644
--- a/estimagic/optimization/nlopt_optimizers.py
+++ b/estimagic/optimization/nlopt_optimizers.py
@@ -962,33 +962,33 @@ def nlopt_isres(
 ):
     """Optimize a scalar function using the ISRES algorithm.
 
-    ISRES is an implementation of "Improved Stochastic Evolution Strategy" written
-    for solving optimization problems with non-linear constraints. The algorithm
-    is supposed to be a global method, in that it has heuristics to avoid local
-    minima. However, no convergence proof is available.
+    ISRES is an implementation of "Improved Stochastic Evolution Strategy"
+    written for solving optimization problems with non-linear constraints. The
+    algorithm is supposed to be a global method, in that it has heuristics to
+    avoid local minima. However, no convergence proof is available.
 
     The original method and a refined version can be found, respecively, in:
-    Thomas Philip Runarsson and Xin Yao, "Search biases in constrained evolutionary
-    optimization," IEEE Trans. on Systems, Man, and Cybernetics Part C: Applications
-    and Reviews, vol. 35 (no. 2), pp. 233-243 (2005).
-    Thomas P. Runarsson and Xin Yao, "Stochastic ranking for constrained evolutionary
-    optimization," IEEE Trans. Evolutionary Computation, vol. 4 (no. 3), pp. 284-294
-    (2000).
+    Thomas Philip Runarsson and Xin Yao, "Search biases in constrained
+    evolutionary optimization," IEEE Trans. on Systems, Man, and Cybernetics
+    Part C: Applications and Reviews, vol. 35 (no. 2), pp. 233-243 (2005).
+    Thomas P. Runarsson and Xin Yao, "Stochastic ranking for constrained
+    evolutionary optimization," IEEE Trans. Evolutionary Computation, vol. 4
+    (no. 3), pp. 284-294 (2000).
 
 
     ``nlopt_isres`` supports the following ``algo_options``:
 
-    - convergence.relative_params_tolerance (float):  Stop when the relative movement
-      between parameter vectors is smaller than this.
-    - convergence.absolute_params_tolerance (float): Stop when the absolute movement
-      between parameter vectors is smaller than this.
+    - convergence.relative_params_tolerance (float):  Stop when the relative
+      movement between parameter vectors is smaller than this.
+    - convergence.absolute_params_tolerance (float): Stop when the absolute
+      movement between parameter vectors is smaller than this.
     - convergence.relative_criterion_tolerance (float): Stop when the relative
       improvement between two iterations is smaller than this.
-    - convergence.absolute_criterion_tolerance (float): Stop when the change of the
-      criterion function between two iterations is smaller than this.
-    - stopping.max_criterion_evaluations (int): If the maximum number of function
-      evaluation is reached, the optimization stops but we do not count this
-      as convergence.
+    - convergence.absolute_criterion_tolerance (float): Stop when the change of
+      the criterion function between two iterations is smaller than this.
+    - stopping.max_criterion_evaluations (int): If the maximum number of
+      function evaluation is reached, the optimization stops but we do not count
+      this as convergence.
 
     """
     out = _minimize_nlopt(
@@ -1018,26 +1018,24 @@ def nlopt_crs2_lm(
     convergence_relative_criterion_tolerance=CONVERGENCE_RELATIVE_CRITERION_TOLERANCE,
     convergence_absolute_criterion_tolerance=CONVERGENCE_ABSOLUTE_CRITERION_TOLERANCE,
     stopping_max_criterion_evaluations=STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL,
-    genetic_popsize=None,
+    population_size=None,
 ):
     """Optimize a scalar function using the CRS2_LM algorithm.
 
-    This implementation of controlled random search method with local mutation is
-    based on:
-    P. Kaelo and M. M. Ali, "Some variants of the controlled random search algorithm
+    This implementation of controlled random search method with local mutation is based
+    on: P. Kaelo and M. M. Ali, "Some variants of the controlled random search algorithm
     for global optimization," J. Optim. Theory Appl. 130 (2), 253-264 (2006).
 
-    The original CRS method is described in:
-    W. L. Price, "A controlled random search procedure for global optimization,"
-    in Towards Global Optimization 2, p. 71-84 edited by L. C. W. Dixon and G. P.
-    Szego (North-Holland Press, Amsterdam, 1978).
-    W. L. Price, "Global optimization by controlled random search," J. Optim. Theory
-    Appl. 40 (3), p. 333-348 (1983).
+    The original CRS method is described in: W. L. Price, "A controlled random search
+    procedure for global optimization," in Towards Global Optimization 2, p. 71-84
+    edited by L. C. W. Dixon and G. P. Szego (North-Holland Press, Amsterdam, 1978). W.
+    L. Price, "Global optimization by controlled random search," J. Optim. Theory Appl.
+    40 (3), p. 333-348 (1983).
 
     CRS class of algorithms starts with random population of points and evolves the
     points "randomly". The size of the initial population can be set via the param-
-    meter genetic_popsize. If the user doesn't specify a value, it is
-    set to the nlopt default of 10*(n+1).
+    meter population_size. If the user doesn't specify a value, it is set to the nlopt
+    default of 10*(n+1).
 
     ``nlopt_isres`` supports the following ``algo_options``:
 
@@ -1050,14 +1048,14 @@ def nlopt_crs2_lm(
     - convergence.absolute_criterion_tolerance (float): Stop when the change of the
       criterion function between two iterations is smaller than this.
     - stopping.max_criterion_evaluations (int): If the maximum number of function
-      evaluation is reached, the optimization stops but we do not count this
-      as convergence.
-    - genetic_popsize(int): The size of the population of the starting
-      points.
+      evaluation is reached, the optimization stops but we do not count this as
+      convergence.
+    - population_size (int): Size of the population. If None, it's set to be
+      10 * (number of parameters + 1).
 
     """
-    if not genetic_popsize:
-        genetic_popsize = 10 * (len(x) + 1)
+    if population_size is None:
+        population_size = 10 * (len(x) + 1)
     out = _minimize_nlopt(
         criterion_and_derivative,
         x,
@@ -1070,7 +1068,7 @@ def nlopt_crs2_lm(
         convergence_ftol_rel=convergence_relative_criterion_tolerance,
         convergence_ftol_abs=convergence_absolute_criterion_tolerance,
         stopping_max_eval=stopping_max_criterion_evaluations,
-        population_size=genetic_popsize,
+        population_size=population_size,
     )
     return out
 
diff --git a/estimagic/optimization/pygmo_optimizers.py b/estimagic/optimization/pygmo_optimizers.py
new file mode 100644
index 0000000..cf2df43
--- /dev/null
+++ b/estimagic/optimization/pygmo_optimizers.py
@@ -0,0 +1,1900 @@
+"""Implement pygmo optimizers."""
+import functools
+import warnings
+
+import numpy as np
+
+from estimagic import batch_evaluators
+from estimagic.config import IS_PYGMO_INSTALLED
+from estimagic.optimization.algo_options import CONVERGENCE_RELATIVE_PARAMS_TOLERANCE
+from estimagic.optimization.algo_options import (
+    STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL,
+)
+
+STOPPING_MAX_ITERATIONS_GENETIC = 250
+
+try:
+    import pygmo as pg
+except ImportError:
+    pass
+
+
+def pygmo_gaco(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    kernel_size=63,
+    speed_parameter_q=1.0,
+    oracle=0.0,
+    accuracy=0.01,
+    threshold=1,
+    speed_of_std_values_convergence=7,
+    stopping_max_n_without_improvements=100000,
+    stopping_max_criterion_evaluations=STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL,
+    focus=0.0,
+    cache=False,
+):
+    """Minimize a scalar function using the generalized ant colony algorithm.
+
+    The version available through pygmo is an generalized version of the
+    original ant colony algorithm proposed by :cite:`Schlueter2009`.
+
+    This algorithm can be applied to box-bounded problems.
+
+    Ant colony optimization is a class of optimization algorithms modeled on the
+    actions of an ant colony. Artificial "ants" (e.g. simulation agents) locate
+    optimal solutions by moving through a parameter space representing all
+    possible solutions. Real ants lay down pheromones directing each other to
+    resources while exploring their environment. The simulated "ants" similarly
+    record their positions and the quality of their solutions, so that in later
+    simulation iterations more ants locate better solutions.
+
+    The generalized ant colony algorithm generates future generations of ants by
+    using a multi-kernel gaussian distribution based on three parameters (i.e.,
+    pheromone values) which are computed depending on the quality of each
+    previous solution. The solutions are ranked through an oracle penalty
+    method.
+
+    - population_size (int): Size of the population. If None, it's twice the
+      number of parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch
+      evaluator (currently 'joblib' and 'pathos_mp') or Callable with the same
+      interface as the estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed
+      to be part of the initial population. This saves one criterion function
+      evaluation that cannot be done in parallel with other evaluations. Default
+      False.
+
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - kernel_size (int): Number of solutions stored in the solution archive.
+    - speed_parameter_q (float): This parameter manages the convergence speed
+      towards the found minima (the smaller the faster). In the pygmo
+      documentation it is referred to as $q$. It must be positive and can be
+      larger than 1. The default is 1.0 until **threshold** is reached. Then it
+      is set to 0.01.
+    - oracle (float): oracle parameter used in the penalty method.
+    - accuracy (float): accuracy parameter for maintaining a minimum penalty
+      function's values distances.
+    - threshold (int): when the iteration counter reaches the threshold the
+      convergence speed is set to 0.01 automatically. To deactivate this effect
+      set the threshold to stopping.max_iterations which is the largest allowed
+      value.
+    - speed_of_std_values_convergence (int): parameter that determines the
+      convergence speed of the standard deviations. This must be an integer
+      (`n_gen_mark` in pygmo and pagmo).
+    - stopping.max_n_without_improvements (int): if a positive integer is
+      assigned here, the algorithm will count the runs without improvements, if
+      this number exceeds the given value, the algorithm will be stopped.
+    - stopping.max_criterion_evaluations (int): maximum number of function
+      evaluations.
+    - focus (float): this parameter makes the search for the optimum greedier
+      and more focused on local improvements (the higher the greedier). If the
+      value is very high, the search is more focused around the current best
+      solutions. Values larger than 1 are allowed.
+    - cache (bool): if True, memory is activated in the algorithm for multiple calls.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+
+    if isinstance(speed_of_std_values_convergence, float):
+        if not speed_of_std_values_convergence.is_integer():
+            raise ValueError(
+                "The speed_of_std_values_convergence parameter must be an integer. "
+                f"You specified {speed_of_std_values_convergence}."
+            )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "ker": kernel_size,
+        "q": speed_parameter_q,
+        "oracle": oracle,
+        "acc": accuracy,
+        "threshold": threshold,
+        "n_gen_mark": int(speed_of_std_values_convergence),
+        "impstop": stopping_max_n_without_improvements,
+        "evalstop": stopping_max_criterion_evaluations,
+        "focus": focus,
+        "memory": cache,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="gaco",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_bee_colony(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    max_n_trials=1,
+    population_size=None,
+):
+    """Minimize a scalar function using the artifical bee colony algorithm.
+
+    The Artificial Bee Colony Algorithm was originally proposed by
+    :cite:`Karaboga2007`. The implemented version of the algorithm is proposed
+    in :cite:`Mernik2015`. The algorithm is only suited for bounded parameter
+    spaces.
+
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch
+      evaluator (currently 'joblib' and 'pathos_mp') or Callable with the same
+      interface as the estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed
+      to be part of the initial population. This saves one criterion function
+      evaluation that cannot be done in parallel with other evaluations. Default
+      False.
+    - max_n_trials (int): Maximum number of trials for abandoning a source.
+      Default is 1.
+    - population_size (int): Size of the population. If None, it's twice the
+      number of parameters but at least 20.
+
+    """
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=20
+    )
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options={
+            "limit": max_n_trials,
+            "gen": int(stopping_max_iterations),
+        },
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="bee_colony",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_de(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    weight_coefficient=0.8,
+    crossover_probability=0.9,
+    mutation_variant="rand/1/exp",
+    convergence_criterion_tolerance=1e-6,
+    convergence_relative_params_tolerance=CONVERGENCE_RELATIVE_PARAMS_TOLERANCE,
+):
+    """Minimize a scalar function using the differential evolution algorithm.
+
+    Differential Evolution is a heuristic optimizer originally presented in
+    :cite:`Storn1997`. The algorithm is only suited for bounded parameter
+    spaces.
+
+    - population_size (int): Size of the population. If None, it's twice the
+      number of parameters but at least 10.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch
+      evaluator (currently 'joblib' and 'pathos_mp') or Callable with the same
+      interface as the estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed
+      to be part of the initial population. This saves one criterion function
+      evaluation that cannot be done in parallel with other evaluations. Default
+      False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - weight_coefficient (float): Weight coefficient. It is denoted by $F$ in
+      the main paper and must lie in [0, 2]. It controls the amplification of
+      the differential variation $(x_{r_2, G} - x_{r_3, G})$.
+    - crossover_probability (float): Crossover probability.
+    - mutation_variant (str or int): code for the mutation variant to create a
+      new candidate individual. The default is . The following are available:
+        - "best/1/exp" (1, when specified as int)
+        - "rand/1/exp" (2, when specified as int)
+        - "rand-to-best/1/exp" (3, when specified as int)
+        - "best/2/exp" (4, when specified as int)
+        - "rand/2/exp" (5, when specified as int)
+        - "best/1/bin" (6, when specified as int)
+        - "rand/1/bin" (7, when specified as int)
+        - "rand-to-best/1/bin" (8, when specified as int)
+        - "best/2/bin" (9, when specified as int)
+        - "rand/2/bin" (10, when specified as int)
+    - convergence.criterion_tolerance: stopping criteria on the criterion
+      tolerance. Default is 1e-6. It is not clear whether this is the absolute
+      or relative criterion tolerance.
+    - convergence.relative_params_tolerance: stopping criteria on the x
+      tolerance. In pygmo the default is 1e-6 but we use our default value of
+      1e-5.
+
+    """
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=10
+    )
+
+    # support both integer and string specification of the mutation variant
+    mutation_variant_str_to_int = {
+        "best/1/exp": 1,
+        "rand/1/exp": 2,
+        "rand-to-best/1/exp": 3,
+        "best/2/exp": 4,
+        "rand/2/exp": 5,
+        "best/1/bin": 6,
+        "rand/1/bin": 7,
+        "rand-to-best/1/bin": 8,
+        "best/2/bin": 9,
+        "rand/2/bin": 10,
+    }
+    mutation_variant = _convert_str_to_int(
+        str_to_int=mutation_variant_str_to_int, value=mutation_variant
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "F": weight_coefficient,
+        "CR": crossover_probability,
+        "variant": mutation_variant,
+        "ftol": convergence_criterion_tolerance,
+        "xtol": convergence_relative_params_tolerance,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="de",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_sea(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    stopping_max_iterations=10_000,  # Each generation will compute the objective once
+):
+    r"""Minimize a scalar function using the (N+1)-ES simple evolutionary algorithm.
+
+    This algorithm represents the simplest evolutionary strategy, where a population of
+    $\lambda$ individuals at each generation produces one offspring by mutating its best
+    individual uniformly at random within the bounds. Should the offspring be better
+    than the worst individual in the population it will substitute it.
+
+    See :cite:`Oliveto2007`.
+
+    The algorithm is only suited for bounded parameter spaces.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 10.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping_max_iterations (int): number of generations to consider. Each generation
+      will compute the objective function once.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=10
+    )
+
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options={"gen": int(stopping_max_iterations)},
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="sea",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_sga(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    crossover_probability=0.9,
+    crossover_strategy="exponential",
+    eta_c=None,
+    mutation_probability=0.02,
+    mutation_strategy="polynomial",
+    mutation_polynomial_distribution_index=None,
+    mutation_gaussian_width=None,
+    selection_strategy="tournament",
+    selection_truncated_n_best=None,
+    selection_tournament_size=None,
+):
+    """Minimize a scalar function using a simple genetic algorithm.
+
+    A detailed description of the algorithm can be found `in the pagmo2 documentation
+    <https://esa.github.io/pagmo2/docs/cpp/algorithms/sga.html>`_.
+
+    See also: :cite:`Oliveto2007`.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - crossover_probability (float): Crossover probability.
+    - crossover_strategy (str): the crossover strategy. One of “exponential”,“binomial”,
+      “single” or “sbx”. Default is "exponential".
+    - eta_c (float): distribution index for “sbx” crossover. This is an inactive
+      parameter if other types of crossovers are selected. Can be in [1, 100].
+    - mutation_probability (float): Mutation probability.
+    - mutation_strategy (str): Mutation strategy. Must be "gaussian", "polynomial" or
+      "uniform". Default is "polynomial".
+    - mutation_polynomial_distribution_index (float): Must be in [0, 1]. Default is 1.
+    - mutation_gaussian_width (float): Must be in [0, 1]. Default is 1.
+    - selection_strategy (str): Selection strategy. Must be "tournament" or "truncated".
+    - selection_truncated_n_best (int): number of best individuals to use in the
+      "truncated" selection mechanism.
+    - selection_tournament_size (int): size of the tournament in the "tournament"
+      selection mechanism. Default is 1.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+
+    if eta_c is not None and crossover_strategy != "sbx":
+        warnings.warn(
+            f"You specified crossover strategy {crossover_strategy} and eta_c. "
+            "However, eta_c is ignored because it is only used when the "
+            "crossover_strategy is set to sbx."
+        )
+    eta_c = 1.0 if eta_c is None else eta_c
+
+    if (
+        mutation_polynomial_distribution_index is not None
+    ) and mutation_strategy != "polynomial":
+        warnings.warn(
+            "You specified a mutation_polynomial_distribution_index but did not choose "
+            "polynomial as your mutation_strategy. Thus, "
+            "mutation_polynomial_distribution_index will be ignored."
+        )
+    if mutation_gaussian_width is not None and mutation_strategy != "gaussian":
+        warnings.warn(
+            "You specified a mutation_gaussian_width but did not choose gaussion as "
+            "your mutation_strategy. Thus, mutation_gaussian_width will be ignored."
+        )
+    if selection_strategy != "truncated" and selection_truncated_n_best is not None:
+        warnings.warn(
+            "You specified selection_truncated_n_best but did not specify truncated as "
+            "your selection strategy. Therefore, selection_truncated_n_best is ignored."
+        )
+    if selection_strategy != "tournament" and selection_tournament_size is not None:
+        warnings.warn(
+            "You specified selection_tournament_size but did not specify tournament as "
+            "your selection strategy. Therefore, selection_tournament_size is ignored."
+        )
+
+    if mutation_strategy == "gaussian" and mutation_gaussian_width is not None:
+        param_m = mutation_gaussian_width
+    elif (
+        mutation_strategy == "polynomial"
+        and mutation_polynomial_distribution_index is not None
+    ):
+        param_m = mutation_polynomial_distribution_index
+    else:
+        param_m = 1.0
+
+    if selection_strategy == "truncated" and selection_truncated_n_best is not None:
+        param_s = selection_truncated_n_best
+    elif selection_strategy == "tournament" and selection_tournament_size is not None:
+        param_s = selection_tournament_size
+    else:
+        param_s = 2
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "cr": crossover_probability,
+        "eta_c": eta_c,
+        "m": mutation_probability,
+        "param_m": param_m,
+        "crossover": crossover_strategy,
+        "mutation": mutation_strategy,
+        "selection": selection_strategy,
+        "param_s": param_s,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="sga",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_sade(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    jde=True,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    mutation_variant="rand/1/exp",
+    keep_adapted_params=False,
+    ftol=1e-6,
+    xtol=1e-6,
+):
+    """Minimize a scalar function using Self-adaptive Differential Evolution.
+
+    The original Differential Evolution algorithm (pygmo_de) can be significantly
+    improved introducing the idea of parameter self-adaptation.
+
+    Many different proposals have been made to self-adapt both the crossover and the
+    F parameters of the original differential evolution algorithm. pygmo's
+    implementation supports two different mechanisms. The first one, proposed by
+    :cite:`Brest2006`, does not make use of the differential evolution operators to
+    produce new values for the weight coefficient $F$ and the crossover probability
+    $CR$ and, strictly speaking, is thus not self-adaptation, rather parameter control.
+    The resulting differential evolution variant is often referred to as jDE.
+    The second variant is inspired by the ideas introduced by :cite:`Elsayed2011` and
+    uses a variaton of the selected DE operator to produce new $CR$ anf $F$ parameters
+    for each individual. This variant is referred to iDE.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - jde (bool): Whether to use the jDE self-adaptation variant to control the $F$ and
+      $CR$ parameter. If True jDE is used, else iDE.
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - mutation_variant (int or str): code for the mutation variant to create a new
+      candidate individual. The default is "rand/1/exp". The first ten are the
+      classical mutation variants introduced in the orginal DE algorithm, the remaining
+      ones are, instead, considered in the work by :cite:`Elsayed2011`.
+      The following are available:
+
+        - "best/1/exp" or 1
+        - "rand/1/exp" or 2
+        - "rand-to-best/1/exp" or 3
+        - "best/2/exp" or 4
+        - "rand/2/exp" or 5
+        - "best/1/bin" or 6
+        - "rand/1/bin" or 7
+        - "rand-to-best/1/bin" or 8
+        - "best/2/bin" or 9
+        - "rand/2/bin" or 10
+        - "rand/3/exp" or 11
+        - "rand/3/bin" or 12
+        - "best/3/exp" or 13
+        - "best/3/bin" or 14
+        - "rand-to-current/2/exp" or 15
+        - "rand-to-current/2/bin" or 16
+        - "rand-to-best-and-current/2/exp" or 17
+        - "rand-to-best-and-current/2/bin" or 18
+
+    - keep_adapted_params (bool):  when true the adapted parameters $CR$ anf $F$ are
+      not reset between successive calls to the evolve method. Default is False.
+    - ftol (float): stopping criteria on the x tolerance.
+    - xtol (float): stopping criteria on the f tolerance.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+    mutation_variant_str_to_int = {
+        "best/1/exp": 1,
+        "rand/1/exp": 2,
+        "rand-to-best/1/exp": 3,
+        "best/2/exp": 4,
+        "rand/2/exp": 5,
+        "best/1/bin": 6,
+        "rand/1/bin": 7,
+        "rand-to-best/1/bin": 8,
+        "best/2/bin": 9,
+        "rand/2/bin": 10,
+        "rand/3/exp": 11,
+        "rand/3/bin": 12,
+        "best/3/exp": 13,
+        "best/3/bin": 14,
+        "rand-to-current/2/exp": 15,
+        "rand-to-current/2/bin": 16,
+        "rand-to-best-and-current/2/exp": 17,
+        "rand-to-best-and-current/2/bin": 18,
+    }
+    mutation_variant = _convert_str_to_int(
+        str_to_int=mutation_variant_str_to_int, value=mutation_variant
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "variant": mutation_variant,
+        "variant_adptv": 1 if jde else 2,
+        "ftol": ftol,
+        "xtol": xtol,
+        "memory": keep_adapted_params,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="sade",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_cmaes(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    backward_horizon=None,
+    variance_loss_compensation=None,
+    learning_rate_rank_one_update=None,
+    learning_rate_rank_mu_update=None,
+    initial_step_size=0.5,
+    ftol=1e-6,
+    xtol=1e-6,
+    keep_adapted_params=False,
+):
+    r"""Minimize a scalar function using the Covariance Matrix Evolutionary Strategy.
+
+    CMA-ES is one of the most successful algorithm, classified as an Evolutionary
+    Strategy, for derivative-free global optimization. The version supported by
+    estimagic is the version described in :cite:`Hansen2006`.
+
+    In contrast to the pygmo version, estimagic always sets force_bounds to True. This
+    avoids that ill defined parameter values are evaluated.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - backward_horizon (float): backward time horizon for the evolution path. It must
+      lie betwen 0 and 1.
+    - variance_loss_compensation (float): makes partly up for the small variance loss in
+      case the indicator is zero. `cs` in the MATLAB Code of :cite:`Hansen2006`. It must
+      lie between 0 and 1.
+    - learning_rate_rank_one_update (float): learning rate for the rank-one update of
+      the covariance matrix. `c1` in the pygmo and pagmo documentation. It must lie
+      between 0 and 1.
+    - learning_rate_rank_mu_update (float): learning rate for the rank-mu update of the
+      covariance matrix. `cmu` in the pygmo and pagmo documentation. It must lie between
+      0 and 1.
+    - initial_step_size (float): initial step size, :math:`\sigma^0` in the original
+      paper.
+    - ftol (float): stopping criteria on the x tolerance.
+    - xtol (float): stopping criteria on the f tolerance.
+    - keep_adapted_params (bool):  when true the adapted parameters are not reset
+      between successive calls to the evolve method. Default is False.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "cc": _replace_none(var=backward_horizon, none_value=-1.0),
+        "cs": _replace_none(var=variance_loss_compensation, none_value=-1.0),
+        "c1": _replace_none(var=learning_rate_rank_one_update, none_value=-1.0),
+        "cmu": _replace_none(var=learning_rate_rank_mu_update, none_value=-1.0),
+        "sigma0": initial_step_size,
+        "ftol": ftol,
+        "xtol": xtol,
+        "memory": keep_adapted_params,
+        "force_bounds": True,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="cmaes",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_simulated_annealing(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    start_temperature=10.0,
+    end_temperature=0.01,
+    n_temp_adjustments=10,
+    n_range_adjustments=10,
+    bin_size=10,
+    start_range=1.0,
+):
+    """Minimize a function with the simulated annealing algorithm.
+
+    This version of the simulated annealing algorithm is, essentially, an iterative
+    random search procedure with adaptive moves along the coordinate directions. It
+    permits uphill moves under the control of metropolis criterion, in the hope to avoid
+    the first local minima encountered. This version is the one proposed in
+    :cite:`Corana1987`.
+
+    .. note: When selecting the starting and final temperature values it helps to think
+        about the tempertaure as the deterioration in the objective function value that
+        still has a 37% chance of being accepted.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - start_temperature (float): starting temperature. Must be > 0.
+    - end_temperature (float): final temperature. Our default (0.01) is lower than in
+      pygmo and pagmo. The final temperature must be positive.
+    - n_temp_adjustments (int): number of temperature adjustments in the annealing
+      schedule.
+    - n_range_adjustments (int): number of adjustments of the search range performed at
+      a constant temperature.
+    - bin_size (int): number of mutations that are used to compute the acceptance rate.
+    - start_range (float): starting range for mutating the decision vector. It must lie
+      between 0 and 1.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+
+    algo_specific_options = {
+        "Ts": start_temperature,
+        "Tf": end_temperature,
+        "n_T_adj": int(n_temp_adjustments),
+        "n_range_adj": int(n_range_adjustments),
+        "bin_size": bin_size,
+        "start_range": start_range,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="simulated_annealing",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_pso(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    omega=0.7298,
+    force_of_previous_best=2.05,
+    force_of_best_in_neighborhood=2.05,
+    max_velocity=0.5,
+    algo_variant=5,
+    neighbor_definition="lbest",
+    neighbor_param=None,
+    keep_velocities=False,
+):
+    r"""Minimize a scalar function using Particle Swarm Optimization.
+
+    Particle swarm optimization (PSO) is a population based algorithm inspired by the
+    foraging behaviour of swarms. In PSO each point has memory of the position where it
+    achieved the best performance xli (local memory) and of the best decision vector
+    :math:`x^g` in a certain neighbourhood, and uses this information to update its
+    position.
+
+    For a survey on particle swarm optimization algorithms, see :cite:`Poli2007`.
+
+    Each particle determines its future position :math:`x_{i+1} = x_i + v_i` where
+
+    .. math:: v_{i+1} = \omega (v_i + \eta_1 \cdot \mathbf{r}_1 \cdot (x_i - x^{l}_i) +
+        \eta_2 \cdot \mathbf{r}_2 \cdot (x_i - x^g))
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 10.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+
+    - omega (float): depending on the variant chosen, :math:`\omega` is the particles'
+      inertia weight or the construction coefficient. It must lie between 0 and 1.
+    - force_of_previous_best (float): :math:`\eta_1` in the equation above. It's the
+      magnitude of the force, applied to the particle’s velocity, in the direction of
+      its previous best position. It must lie between 0 and 4.
+    - force_of_best_in_neighborhood (float): :math:`\eta_2` in the equation above. It's
+      the magnitude of the force, applied to the particle’s velocity, in the direction
+      of the best position in its neighborhood. It must lie between 0 and 4.
+    - max_velocity (float): maximum allowed particle velocity as fraction of the box
+      bounds. It must lie between 0 and 1.
+    - algo_variant (int or str): algorithm variant to be used:
+        - 1 or "canonical_inertia": Canonical (with inertia weight)
+        - 2 or "social_and_cog_rand": Same social and cognitive rand.
+        - 3 or "all_components_rand": Same rand. for all components
+        - 4 or "one_rand": Only one rand.
+        - 5 or "canonical_constriction": Canonical (with constriction fact.)
+        - 6 or "fips": Fully Informed (FIPS)
+
+    - neighbor_definition (int or str): swarm topology that defines each particle's
+      neighbors that is to be used:
+        - 1 or "gbest"
+        - 2 or "lbest"
+        - 3 or "Von Neumann"
+        - 4 or "Adaptive random"
+
+    - neighbor_param (int): the neighbourhood parameter. If the lbest topology is
+      selected (neighbor_definition=2), it represents each particle's indegree (also
+      outdegree) in the swarm topology. Particles have neighbours up to a radius of k =
+      neighbor_param / 2 in the ring. If the Randomly-varying neighbourhood topology is
+      selected (neighbor_definition=4), it represents each particle’s maximum outdegree
+      in the swarm topology. The minimum outdegree is 1 (the particle always connects
+      back to itself). If neighbor_definition is 1 or 3 this parameter is ignored.
+    - keep_velocities (bool): when true the particle velocities are not reset between
+      successive calls to `evolve`.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    if neighbor_definition in [1, 3] and neighbor_param is not None:
+        warnings.warn(
+            "You gave a neighbor parameter but selected a neighbor_definition "
+            "that ignores this parameter."
+        )
+
+    neighbor_param = _replace_none(neighbor_param, 4)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=10
+    )
+
+    neighbor_definition_str_to_int = {
+        "gbest": 1,
+        "lbest": 2,
+        "Von Neumann": 3,
+        "Adaptive random": 4,
+    }
+    algo_variant_str_to_int = {
+        "canonical_inertia": 1,
+        "social_and_cog_rand": 2,
+        "all_components_rand": 3,
+        "one_rand": 4,
+        "canonical_constriction": 5,
+        "fips": 6,
+    }
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "omega": omega,
+        "eta1": force_of_previous_best,
+        "eta2": force_of_best_in_neighborhood,
+        "max_vel": max_velocity,
+        "variant": _convert_str_to_int(algo_variant_str_to_int, algo_variant),
+        "neighb_type": _convert_str_to_int(
+            neighbor_definition_str_to_int, neighbor_definition
+        ),
+        "neighb_param": neighbor_param,
+        "memory": keep_velocities,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="pso",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_pso_gen(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    omega=0.7298,
+    force_of_previous_best=2.05,
+    force_of_best_in_neighborhood=2.05,
+    max_velocity=0.5,
+    algo_variant="canonical_constriction",
+    neighbor_definition=2,
+    neighbor_param=None,
+    keep_velocities=False,
+):
+    r"""Minimize a scalar function with generational Particle Swarm Optimization.
+
+    Particle Swarm Optimization (generational) is identical to pso, but does update the
+    velocities of each particle before new particle positions are computed (taking into
+    consideration all updated particle velocities). Each particle is thus evaluated on
+    the same seed within a generation as opposed to the standard PSO which evaluates
+    single particle at a time. Consequently, the generational PSO algorithm is suited
+    for stochastic optimization problems.
+
+    For a survey on particle swarm optimization algorithms, see :cite:`Poli2007`.
+
+    Each particle determines its future position :math:`x_{i+1} = x_i + v_i` where
+
+    .. math:: v_{i+1} = \omega (v_i + \eta_1 \cdot \mathbf{r}_1 \cdot (x_i - x^{l}_i) +
+        \eta_2 \cdot \mathbf{r}_2 \cdot (x_i - x^g))
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 10.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+
+    - omega (float): depending on the variant chosen, :math:`\omega` is the particles'
+      inertia weight or the constructuion coefficient. It must lie between 0 and 1.
+    - force_of_previous_best (float): :math:`\eta_1` in the equation above. It's the
+      magnitude of the force, applied to the particle’s velocity, in the direction of
+      its previous best position. It must lie between 0 and 4.
+    - force_of_best_in_neighborhood (float): :math:`\eta_2` in the equation above. It's
+      the magnitude of the force, applied to the particle’s velocity, in the direction
+      of the best position in its neighborhood. It must lie between 0 and 4.
+    - max_velocity (float): maximum allowed particle velocity as fraction of the box
+      bounds. It must lie between 0 and 1.
+    - algo_variant (int): code of the algorithm's variant to be used:
+
+        - 1 or "canonical_inertia": Canonical (with inertia weight)
+        - 2 or "social_and_cog_rand": Same social and cognitive rand.
+        - 3 or "all_components_rand": Same rand. for all components
+        - 4 or "one_rand": Only one rand.
+        - 5 or "canonical_constriction": Canonical (with constriction fact.)
+        - 6 or "fips": Fully Informed (FIPS)
+
+    - neighbor_definition (int): code for the swarm topology that defines each
+      particle's neighbors that is to be used:
+
+        - 1 or "gbest"
+        - 2 or "lbest"
+        - 3 or "Von Neumann"
+        - 4 or "Adaptive random"
+
+    - neighbor_param (int): the neighbourhood parameter. If the lbest topology is
+      selected (neighbor_definition=2), it represents each particle's indegree (also
+      outdegree) in the swarm topology. Particles have neighbours up to a radius of k =
+      neighbor_param / 2 in the ring. If the Randomly-varying neighbourhood topology is
+      selected (neighbor_definition=4), it represents each particle’s maximum outdegree
+      in the swarm topology. The minimum outdegree is 1 (the particle always connects
+      back to itself). If neighbor_definition is 1 or 3 this parameter is ignored.
+    - keep_velocities (bool): when true the particle velocities are not reset between
+      successive calls to `evolve`.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    if neighbor_definition in [1, 3] and neighbor_param is not None:
+        warnings.warn(
+            "You gave a neighbor parameter but selected a neighbor_definition "
+            "that ignores this parameter."
+        )
+    neighbor_param = _replace_none(neighbor_param, 4)
+    neighbor_str_to_int = {
+        "gbest": 1,
+        "lbest": 2,
+        "Von Neumann": 3,
+        "Adaptive random": 4,
+    }
+    neighbor_param = _convert_str_to_int(neighbor_str_to_int, neighbor_param)
+
+    algo_variant_str_to_int = {
+        "canonical_inertia": 1,
+        "social_and_cog_rand": 2,
+        "all_components_rand": 3,
+        "one_rand": 4,
+        "canonical_constriction": 5,
+        "fips": 6,
+    }
+    algo_variant = _convert_str_to_int(algo_variant_str_to_int, algo_variant)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=10
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "omega": omega,
+        "eta1": force_of_previous_best,
+        "eta2": force_of_best_in_neighborhood,
+        "max_vel": max_velocity,
+        "variant": algo_variant,
+        "neighb_type": neighbor_definition,
+        "neighb_param": neighbor_param,
+        "memory": keep_velocities,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="pso_gen",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_mbh(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    inner_algorithm=None,
+    # this is 30 instead of 5 in pygmo for our sum of squares test to pass
+    stopping_max_inner_runs_without_improvement=30,
+    perturbation=0.01,
+):
+    """Minimize a scalar function using generalized Monotonic Basin Hopping.
+
+    Monotonic basin hopping, or simply, basin hopping, is an algorithm rooted in the
+    idea of mapping the objective function $f(x_0)$ into the local minima found starting
+    from $x_0$. This simple idea allows a substantial increase of efficiency in solving
+    problems, such as the Lennard-Jones cluster or the MGA-1DSM interplanetary
+    trajectory problem that are conjectured to have a so-called funnel structure.
+
+    See :cite:`Wales1997` for the paper introducing the basin hopping idea for a
+    Lennard-Jones cluster optimization.
+
+    pygmo provides an original generalization of this concept resulting in a
+    meta-algorithm that operates on a population. When a population containing a single
+    individual is used the original method is recovered.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 250.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - inner_algorithm (pygmo.algorithm): an pygmo algorithm or a user-defined algorithm,
+      either C++ or Python. If None the `pygmo.compass_search` algorithm will be used.
+    - stopping.max_inner_runs_without_improvement (int): consecutive runs of the inner
+      algorithm that need to result in no improvement for mbh to stop.
+    - perturbation (float): the perturbation to be applied to each component.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    # the min default population size is this large to pass our sum of squares tests.
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=250
+    )
+
+    algo_specific_options = {
+        "algo": inner_algorithm,
+        "stop": stopping_max_inner_runs_without_improvement,
+        "perturb": perturbation,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="mbh",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_xnes(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    learning_rate_mean_update=1.0,
+    learning_rate_step_size_update=None,
+    learning_rate_cov_matrix_update=None,
+    initial_search_share=1.0,
+    ftol=1e-6,
+    xtol=1e-6,
+    keep_adapted_params=False,
+):
+    r"""Minimize a scalar function using Exponential Evolution Strategies.
+
+    Exponential Natural Evolution Strategies is an algorithm closely related to CMAES
+    and based on the adaptation of a gaussian sampling distribution via the so-called
+    natural gradient. Like CMAES it is based on the idea of sampling new trial vectors
+    from a multivariate distribution and using the new sampled points to update the
+    distribution parameters. Naively this could be done following the gradient of the
+    expected fitness as approximated by a finite number of sampled points. While this
+    idea offers a powerful lead on algorithmic construction it has some major drawbacks
+    that are solved in the so-called Natural Evolution Strategies class of algorithms by
+    adopting, instead, the natural gradient. xNES is one of the most performing variants
+    in this class.
+
+    See :cite:`Glasmachers2010` and the `pagmo documentation on xNES
+    <https://esa.github.io/pagmo2/docs/cpp/algorithms/xnes.html#_CPPv4N5pagmo4xnesE>`_
+    for details.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+
+    - learning_rate_mean_update (float): learning rate for the mean update
+      (:math:`\eta_\mu`). It must be between 0 and 1 or None.
+    - learning_rate_step_size_update (float): learning rate for the step-size update. It
+      must be between 0 and 1 or None.
+    - learning_rate_cov_matrix_update (float): learning rate for the covariance matrix
+      update. It must be between 0 and 1 or None.
+    - initial_search_share (float): share of the given search space that will be
+      initally searched. It must be between 0 and 1. Default is 1.
+    - ftol (float): stopping criteria on the x tolerance.
+    - xtol (float): stopping criteria on the f tolerance.
+    - keep_adapted_params (bool): when true the adapted parameters are not reset between
+      successive calls to the evolve method. Default is False.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "eta_mu": _replace_none(learning_rate_mean_update, -1),
+        "eta_sigma": _replace_none(learning_rate_step_size_update, -1),
+        "eta_b": _replace_none(learning_rate_cov_matrix_update, -1),
+        "sigma0": initial_search_share,
+        "ftol": ftol,
+        "xtol": xtol,
+        "memory": keep_adapted_params,
+        "force_bounds": True,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="xnes",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_gwo(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+):
+    """Minimize a scalar function usinng the Grey Wolf Optimizer.
+
+    The grey wolf optimizer was proposed by :cite:`Mirjalili2014`. The pygmo
+    implementation that is wrapped by estimagic is pased on the pseudo code provided in
+    that paper.
+
+    This algorithm is a classic example of a highly criticizable line of search that led
+    in the first decades of our millenia to the development of an entire zoo of
+    metaphors inspiring optimzation heuristics. In our opinion they, as is the case for
+    the grey wolf optimizer, are often but small variations of already existing
+    heuristics rebranded with unnecessray and convoluted biological metaphors. In the
+    case of GWO this is particularly evident as the position update rule is shokingly
+    trivial and can also be easily seen as a product of an evolutionary metaphor or a
+    particle swarm one. Such an update rule is also not particulary effective and
+    results in a rather poor performance most of times.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options={"gen": int(stopping_max_iterations)},
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="gwo",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_compass_search(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    stopping_max_criterion_evaluations=STOPPING_MAX_CRITERION_EVALUATIONS_GLOBAL,
+    start_range=0.1,
+    stop_range=0.01,
+    reduction_coeff=0.5,
+):
+    """Minimize a scalar function using compass search.
+
+    The algorithm is described in :cite:`Kolda2003`.
+
+    It is considered slow but reliable. It should not be used for stochastic problems.
+
+    - population_size (int): Size of the population. Even though the algorithm is not
+      population based the population size does affect the results of the algorithm.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_criterion_evaluations (int): maximum number of function evaluations.
+    - start_range (float): the start range. Must be in (0, 1].
+    - stop_range (float): the stop range. Must be in (0, start_range].
+    - reduction_coeff (float): the range reduction coefficient. Must be in (0, 1).
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+    if population_size is not None:
+        warnings.warn(
+            f"You specified population size {population_size}. "
+            "compass_search does not have a population so this argument is ignored."
+        )
+    else:
+        # if discard_start_params is False population_size - 1 must still be positive
+        population_size = 100
+
+    algo_specific_options = {
+        "max_fevals": stopping_max_criterion_evaluations,
+        "start_range": start_range,
+        "stop_range": stop_range,
+        "reduction_coeff": reduction_coeff,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="compass_search",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_ihs(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    choose_from_memory_probability=0.85,
+    min_pitch_adjustment_rate=0.35,
+    max_pitch_adjustment_rate=0.99,
+    min_distance_bandwidth=1e-5,
+    max_distance_bandwidth=1.0,
+):
+    """Minimize a scalar function using the improved harmony search algorithm.
+
+    Improved harmony search (IHS) was introduced by :cite:`Mahdavi2007`.
+    IHS supports stochastic problems.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - choose_from_memory_probability (float): probability of choosing from memory
+      (similar to a crossover probability).
+    - min_pitch_adjustment_rate (float): minimum pitch adjustment rate. (similar to a
+      mutation rate). It must be between 0 and 1.
+    - max_pitch_adjustment_rate (float): maximum pitch adjustment rate. (similar to a
+      mutation rate). It must be between 0 and 1.
+    - min_distance_bandwidth (float): minimum distance bandwidth. (similar to a mutation
+      width). It must be positive.
+    - max_distance_bandwidth (float): maximum distance bandwidth. (similar to a mutation
+      width).
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+
+    if population_size is not None:
+        warnings.warn("The population size has no effect on IHS' performance.")
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=1
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "phmcr": choose_from_memory_probability,
+        "ppar_min": min_pitch_adjustment_rate,
+        "ppar_max": max_pitch_adjustment_rate,
+        "bw_min": min_distance_bandwidth,
+        "bw_max": max_distance_bandwidth,
+    }
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="ihs",
+        algo_options=algo_options,
+    )
+    return res
+
+
+def pygmo_de1220(
+    criterion_and_derivative,
+    x,
+    lower_bounds,
+    upper_bounds,
+    *,
+    population_size=None,
+    batch_evaluator=None,
+    n_cores=1,
+    seed=None,
+    discard_start_params=False,
+    #
+    jde=True,
+    stopping_max_iterations=STOPPING_MAX_ITERATIONS_GENETIC,
+    allowed_variants=None,
+    keep_adapted_params=False,
+    ftol=1e-6,
+    xtol=1e-6,
+):
+    """Minimize a scalar function using Self-adaptive Differential Evolution, pygmo flavor.
+
+    See `the PAGMO documentation for details
+    <https://esa.github.io/pagmo2/docs/cpp/algorithms/de1220.html>`_.
+
+    - population_size (int): Size of the population. If None, it's twice the number of
+      parameters but at least 64.
+    - batch_evaluator (str or Callable): Name of a pre-implemented batch evaluator
+      (currently 'joblib' and 'pathos_mp') or Callable with the same interface as the
+      estimagic batch_evaluators. See :ref:`batch_evaluators`.
+    - n_cores (int): Number of cores to use.
+    - seed (int): seed used by the internal random number generator.
+    - discard_start_params (bool): If True, the start params are not guaranteed to be
+      part of the initial population. This saves one criterion function evaluation that
+      cannot be done in parallel with other evaluations. Default False.
+    - jde (bool): Whether to use the jDE self-adaptation variant to control the $F$ and
+      $CR$ parameter. If True jDE is used, else iDE.
+    - stopping.max_iterations (int): Number of generations to evolve.
+    - allowed_variants (array-like object): allowed mutation variants (can be codes
+      or strings). Each code refers to one mutation variant to create a new candidate
+      individual. The first ten refer to the classical mutation variants introduced in
+      the original DE algorithm, the remaining ones are, instead, considered in the work
+      by :cite:`Elsayed2011`. The default is ["rand/1/exp", "rand-to-best/1/exp",
+      "rand/1/bin", "rand/2/bin", "best/3/exp", "best/3/bin", "rand-to-current/2/exp",
+      "rand-to-current/2/bin"]. The following are available:
+        - 1 or "best/1/exp"
+        - 2 or "rand/1/exp"
+        - 3 or "rand-to-best/1/exp"
+        - 4 or "best/2/exp"
+        - 5 or "rand/2/exp"
+        - 6 or "best/1/bin"
+        - 7 or "rand/1/bin"
+        - 8 or "rand-to-best/1/bin"
+        - 9 or "best/2/bin"
+        - 10 or "rand/2/bin"
+        - 11 or "rand/3/exp"
+        - 12 or "rand/3/bin"
+        - 13 or "best/3/exp"
+        - 14 or "best/3/bin"
+        - 15 or "rand-to-current/2/exp"
+        - 16 or "rand-to-current/2/bin"
+        - 17 or "rand-to-best-and-current/2/exp"
+        - 18 or "rand-to-best-and-current/2/bin"
+
+    - keep_adapted_params (bool):  when true the adapted parameters $CR$ anf $F$ are not
+      reset between successive calls to the evolve method. Default is False.
+    - ftol (float): stopping criteria on the x tolerance.
+    - xtol (float): stopping criteria on the f tolerance.
+
+    """
+    _check_that_every_param_is_bounded(lower_bounds, upper_bounds)
+    variant_str_to_int = {
+        "best/1/exp": 1,
+        "rand/1/exp": 2,
+        "rand-to-best/1/exp": 3,
+        "best/2/exp": 4,
+        "rand/2/exp": 5,
+        "best/1/bin": 6,
+        "rand/1/bin": 7,
+        "rand-to-best/1/bin": 8,
+        "best/2/bin": 9,
+        "rand/2/bin": 10,
+        "rand/3/exp": 11,
+        "rand/3/bin": 12,
+        "best/3/exp": 13,
+        "best/3/bin": 14,
+        "rand-to-current/2/exp": 15,
+        "rand-to-current/2/bin": 16,
+        "rand-to-best-and-current/2/exp": 17,
+        "rand-to-best-and-current/2/bin": 18,
+    }
+    if allowed_variants is None:
+        allowed_variant_codes = [2, 3, 7, 10, 13, 14, 15, 16]
+    else:
+        allowed_variant_codes = []
+        for variant in allowed_variants:
+            allowed_variant_codes.append(
+                _convert_str_to_int(variant_str_to_int, variant)
+            )
+
+    population_size = _determine_population_size(
+        population_size=population_size, x=x, lower_bound=64
+    )
+
+    algo_specific_options = {
+        "gen": int(stopping_max_iterations),
+        "variant_adptv": 1 if jde else 2,
+        "ftol": ftol,
+        "xtol": xtol,
+        "memory": keep_adapted_params,
+        "allowed_variants": allowed_variant_codes,
+    }
+
+    algo_options = _create_algo_options(
+        population_size=population_size,
+        n_cores=n_cores,
+        seed=seed,
+        discard_start_params=discard_start_params,
+        batch_evaluator=batch_evaluator,
+        algo_specific_options=algo_specific_options,
+    )
+
+    res = _minimize_pygmo(
+        criterion_and_derivative=criterion_and_derivative,
+        x=x,
+        lower_bounds=lower_bounds,
+        upper_bounds=upper_bounds,
+        method="de1220",
+        algo_options=algo_options,
+    )
+    return res
+
+
+# ====================================================================================
+
+
+def _minimize_pygmo(
+    criterion_and_derivative, x, lower_bounds, upper_bounds, method, algo_options=None
+):
+    """Minimize a function with pygmo.
+
+    Args:
+        criterion_and_derivative (callable):
+        x (np.ndarray): Starting values of the parameters.
+        lower_bounds (np.ndarray):
+        upper_bounds (np.ndarray):
+        method (str): One of the optimizers of the pygmo package.
+        algo_options (dict): Options for the optimizer. In addition to
+            the algo options that will be passed directly to the pygmo
+            algorithms we have the following entries:
+            - population_size (int): Population size for genetic algorithms.
+            - batch_evaluator (str or callable): An estimagic batch evaluator,
+                default joblib batch evaluator.
+            - n_cores (int): Number of cores used for parallel evaluation of
+                the criterion function. Default 1.
+            - seed (int or None): Random seed for drawing the initial
+                population.
+            - discard_start_params (bool): If True, the start params are not
+                guaranteed to be part of the initial population. This saves one
+                criterion function evaluation that cannot be done in parallel
+                with other evaluations. Default False.
+
+    Returns:
+        results (dict): Dictionary with optimization results.
+
+    """
+    if not IS_PYGMO_INSTALLED:
+        raise NotImplementedError(
+            f"The pygmo package is not installed and required for '{method}'."
+        )
+
+    algo_options = {} if algo_options is None else algo_options.copy()
+    population_size = algo_options.pop("population_size", 1)
+
+    batch_evaluator = algo_options.pop("batch_evaluator", "joblib_batch_evaluator")
+    if isinstance(batch_evaluator, str):
+        batch_evaluator = getattr(batch_evaluators, batch_evaluator)
+    n_cores = algo_options.pop("n_cores", 1)
+    seed = algo_options.pop("seed", None)
+    discard_start_params = algo_options.pop("discard_start_params", False)
+
+    algo_info = {
+        "parallelizes": n_cores != 1,
+        "name": f"pygmo_{method}",
+        "needs_scaling": True,
+        "primary_criterion_entry": "value",
+    }
+    func = functools.partial(
+        criterion_and_derivative, task="criterion", algorithm_info=algo_info
+    )
+    gradient = functools.partial(
+        criterion_and_derivative, task="derivative", algorithm_info=algo_info
+    )
+
+    bounds = (lower_bounds, upper_bounds)
+    prob = _create_problem(
+        func=func,
+        bounds=bounds,
+        gradient_=gradient,
+        dim=len(x),
+        batch_evaluator=batch_evaluator,
+        n_cores=n_cores,
+    )
+    algo = _create_algorithm(method, algo_options, n_cores)
+    pop = _create_population(
+        prob, population_size, x, seed=seed, discard_start_params=discard_start_params
+    )
+    evolved = algo.evolve(pop)
+    result = _process_pygmo_results(evolved)
+
+    return result
+
+
+def _create_problem(func, bounds, gradient_, dim, batch_evaluator, n_cores):
+    class Problem:
+        def fitness(self, x):
+            return [func(x)]
+
+        def get_bounds(self):
+            return bounds
+
+        def gradient(self, dv):
+            return gradient_(dv)
+
+        def batch_fitness(self, dvs):
+            dv_list = dvs.reshape(-1, dim)
+            eval_list = batch_evaluator(
+                func=func,
+                arguments=dv_list,
+                n_cores=n_cores,
+                # Error handling is done on a higher level
+                error_handling="raise",
+            )
+            evals = np.array(eval_list)
+            return evals
+
+    problem = pg.problem(Problem())
+    return problem
+
+
+def _create_algorithm(method, algo_options, n_cores):
+    """Create a pygmo algorithm."""
+    pygmo_uda = getattr(pg, method)
+    algo = pygmo_uda(**algo_options)
+    try:
+        algo.set_bfe(pg.bfe())
+    except AttributeError:
+        if n_cores >= 2:
+            warnings.warn(
+                f"Your specified algorithm {method} does not support parallelization. "
+                "Choose another algorithm such as pygmo_gaco to parallelize."
+            )
+    out = pg.algorithm(algo)
+    return out
+
+
+def _create_population(problem, population_size, x, seed, discard_start_params):
+    """Create a pygmo population object.
+    Args:
+        problem (pygmo.Problem)
+        algo_options (dict)
+        x (np.ndarray)
+    Todo:
+        - constrain random initial values to be in some bounds
+    """
+    if not discard_start_params:
+        population_size = population_size - 1
+
+    pop = pg.population(
+        problem,
+        size=population_size,
+        seed=seed,
+        b=pg.bfe(),
+    )
+    if not discard_start_params:
+        pop.push_back(x)
+    return pop
+
+
+def _process_pygmo_results(evolved):
+    results = {
+        # Harmonized results.
+        "solution_x": evolved.champion_x,
+        "solution_criterion": evolved.champion_f[0],
+        "solution_derivative": None,
+        "solution_hessian": None,
+        "n_criterion_evaluations": evolved.problem.get_fevals(),
+        "n_derivative_evaluations": evolved.problem.get_gevals(),
+        "n_iterations": None,
+        "success": True,
+        "reached_convergence_criterion": "Number of generations reached.",
+        "message": None,
+    }
+
+    return results
+
+
+def _create_algo_options(
+    population_size,
+    n_cores,
+    seed,
+    discard_start_params,
+    batch_evaluator,
+    algo_specific_options,
+):
+    algo_options = {
+        "population_size": _replace_none(population_size, -1),
+        "n_cores": n_cores,
+        "seed": seed,
+        "discard_start_params": discard_start_params,
+    }
+    if batch_evaluator is not None:
+        algo_options["batch_evaluator"] = batch_evaluator
+
+    algo_options.update(**algo_specific_options)
+    return algo_options
+
+
+def _check_that_every_param_is_bounded(lower_bounds, upper_bounds):
+    assert np.isfinite(lower_bounds).all(), "The lower bounds must all be finite."
+    assert np.isfinite(upper_bounds).all(), "The upper bounds must all be finite."
+
+
+def _determine_population_size(population_size, x, lower_bound):
+    if population_size is None:
+        population_size = int(np.clip(10 * (len(x) + 1), lower_bound, np.inf))
+    else:
+        population_size = int(population_size)
+    return population_size
+
+
+def _convert_str_to_int(str_to_int, value):
+    if value in str_to_int.keys():
+        out = str_to_int[value]
+    elif value not in str_to_int.values():
+        raise ValueError(
+            f"You specified {value} as value. "
+            f"It must be one of {', '.join(str_to_int.keys())}"
+        )
+    else:
+        out = value
+    return out
+
+
+def _replace_none(var, none_value):
+    out = var if var is not None else none_value
+    return out
diff --git a/tox.ini b/tox.ini
index cc7371b..4b4f516 100644
--- a/tox.ini
+++ b/tox.ini
@@ -36,6 +36,7 @@ conda_deps =
     dill
     cyipopt
     nlopt
+    pygmo
 commands = pytest {posargs}
 
 
@@ -95,6 +96,7 @@ filterwarnings =
     ignore: the imp module is deprecated
     ignore: indexing past lexsort depth may impact performance.
     ignore: Method .ptp is deprecated and will be removed in a future version. Use numpy.ptp instead.
+    ignore: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only
 markers =
     wip: Tests that are work-in-progress.
     slow: Tests that take a long time to run and are skipped in continuous integration.


Test patch:
diff --git a/estimagic/tests/optimization/test_many_algorithms.py b/estimagic/tests/optimization/test_many_algorithms.py
index 7421659..ad03d7b 100644
--- a/estimagic/tests/optimization/test_many_algorithms.py
+++ b/estimagic/tests/optimization/test_many_algorithms.py
@@ -5,6 +5,7 @@
 
 """
 import inspect
+import sys
 
 import numpy as np
 import pandas as pd
@@ -64,17 +65,22 @@ def test_algorithm_on_sum_of_squares_with_binding_bounds(algorithm):
     aaae(res["solution_params"]["value"].to_numpy(), np.array([1, 0, -1]), decimal=3)
 
 
+skip_msg = (
+    "The very slow tests of global algorithms are only run on linux which always "
+    "runs much faster in continuous integration."
+)
+
+
+@pytest.mark.skipif(sys.platform != "linux", reason=skip_msg)
 @pytest.mark.parametrize("algorithm", GLOBAL_ALGORITHMS_AVAILABLE)
 def test_global_algorithms_on_sum_of_squares(algorithm):
     params = pd.DataFrame()
-    params["value"] = [2.5, 2, -1]
-    params["lower_bound"] = [0.5, -1, -1.5]
-    params["upper_bound"] = [3, 3, -0.5]
+    params["value"] = [0.35, 0.35]
+    params["lower_bound"] = [0.2, -0.5]
+    params["upper_bound"] = [1, 0.5]
     res = minimize(
         criterion=sos_dict_criterion,
         params=params,
         algorithm=algorithm,
     )
-    aaae(
-        res["solution_params"]["value"].to_numpy(), np.array([0.5, 0, -0.5]), decimal=1
-    )
+    aaae(res["solution_params"]["value"].to_numpy(), np.array([0.2, 0]), decimal=1)
diff --git a/estimagic/tests/optimization/test_optimization_helpers.py b/estimagic/tests/optimization/test_optimization_helpers.py
new file mode 100644
index 0000000..4509e6d
--- /dev/null
+++ b/estimagic/tests/optimization/test_optimization_helpers.py
@@ -0,0 +1,49 @@
+"""Test optimization helper functions."""
+import numpy as np
+import pytest
+
+from estimagic.optimization.pygmo_optimizers import _check_that_every_param_is_bounded
+from estimagic.optimization.pygmo_optimizers import _convert_str_to_int
+from estimagic.optimization.pygmo_optimizers import _determine_population_size
+from estimagic.optimization.pygmo_optimizers import _replace_none
+
+
+def test_check_that_every_param_is_bounded():
+    _check_that_every_param_is_bounded([1, 2, 3], [4, 5, 6])
+    with pytest.raises(AssertionError):
+        _check_that_every_param_is_bounded([-np.inf, 2, 3], [4, 5, 6])
+    with pytest.raises(AssertionError):
+        _check_that_every_param_is_bounded([np.nan, 2, 3], [4, 5, 6])
+    with pytest.raises(AssertionError):
+        _check_that_every_param_is_bounded([np.inf, 2, 3], [4, 5, 6])
+
+
+test_cases = [
+    # popsize, x, lower_bound, expected
+    (55.3, None, None, 55),
+    (None, np.ones(5), 500, 500),
+    (None, np.ones(5), 4, 60),
+]
+
+
+@pytest.mark.parametrize("popsize, x, lower_bound, expected", test_cases)
+def test_determine_population_size(popsize, x, lower_bound, expected):
+    res = _determine_population_size(
+        population_size=popsize, x=x, lower_bound=lower_bound
+    )
+    assert res == expected
+
+
+def test_convert_str_to_int():
+    d = {"a": 1, "b": 3}
+    assert _convert_str_to_int(d, "a") == 1
+    assert _convert_str_to_int(d, 1) == 1
+    with pytest.raises(ValueError):
+        _convert_str_to_int(d, 5)
+    with pytest.raises(ValueError):
+        _convert_str_to_int(d, "hello")
+
+
+def test_replace_none():
+    assert _replace_none(None, 5) == 5
+    assert _replace_none("a", 10) == "a"
diff --git a/estimagic/tests/test_utilities.py b/estimagic/tests/test_utilities.py
index f0d0a63..0316ddb 100644
--- a/estimagic/tests/test_utilities.py
+++ b/estimagic/tests/test_utilities.py
@@ -7,6 +7,7 @@ from estimagic.config import IS_CYIPOPT_INSTALLED
 from estimagic.config import IS_DFOLS_INSTALLED
 from estimagic.config import IS_PETSC4PY_INSTALLED
 from estimagic.config import IS_PYBOBYQA_INSTALLED
+from estimagic.config import IS_PYGMO_INSTALLED
 from estimagic.optimization import AVAILABLE_ALGORITHMS
 from estimagic.utilities import calculate_trustregion_initial_radius
 from estimagic.utilities import chol_params_to_lower_triangular_matrix
@@ -180,5 +181,6 @@ def test_available_algorithms():
     assert ("nag_dfols" in present_algo_names) is IS_DFOLS_INSTALLED
     assert ("tao_pounders" in present_algo_names) is IS_PETSC4PY_INSTALLED
     assert ("nag_pybobyqa" in present_algo_names) is IS_PYBOBYQA_INSTALLED
+    assert ("pygmo_gaco" in present_algo_names) is IS_PYGMO_INSTALLED
     assert ("ipopt" in present_algo_names) is IS_CYIPOPT_INSTALLED
     assert "get_scipy_bounds" not in present_algo_names