Repository: scipy/scipy
Pull request: ENH: x0 for differential_evolution
closes #13300 

Problem statement:
Add a single individual to a latinhypercube initial population in differential_evolution
**Is your feature request related to a problem? Please describe.**
I'd like to pass a single individual to `differential_evolution`, where it still starts with a `latinhypercube`, but adds the individual. Currently `init=array` requires an array of size `(M, len(x))` where `M>=4`. This seems also to skip the `latinhypercube`, and I have only a single individual.

Context is that this individual is the optimal solution from a related optimization. I don't want the search space constrained, but I want this related optimum to be part of the population.

**Describe the solution you'd like**
New arg that defines additional individual(s) to add to the population, while still respecting the `init` arg and allowing fewer than 4 individuals.

Hints from issue discussion:
I've thought about this kind of thing before. What can happen is that the solve process spends a lot of time converging around the supplied `x0` to start with, unless `x0` is a really bad solution. If `x0` is a pretty good global solution then it's not clear to me that `differential_evolution` will necessarily improve on that. Moreover one wouldn't want to use any `best*` strategies.

Commit message: MAINT: ensure parameters are within bounds

Patch to review:
diff --git a/scipy/optimize/_differentialevolution.py b/scipy/optimize/_differentialevolution.py
index 34448100beea..e54babb1f42a 100644
--- a/scipy/optimize/_differentialevolution.py
+++ b/scipy/optimize/_differentialevolution.py
@@ -24,7 +24,7 @@ def differential_evolution(func, bounds, args=(), strategy='best1bin',
                            mutation=(0.5, 1), recombination=0.7, seed=None,
                            callback=None, disp=False, polish=True,
                            init='latinhypercube', atol=0, updating='immediate',
-                           workers=1, constraints=()):
+                           workers=1, constraints=(), x0=None):
     """Finds the global minimum of a multivariate function.
 
     Differential Evolution is stochastic in nature (does not use gradient
@@ -173,6 +173,13 @@ def differential_evolution(func, bounds, args=(), strategy='best1bin',
 
         .. versionadded:: 1.4.0
 
+    x0 : None or array-like, optional
+        Provides an initial guess to the minimization. Once the population has
+        been initialized this vector replaces the first (best) member. This
+        replacement is done even if `init` is given an initial population.
+
+        .. versionadded:: 1.7.0
+
     Returns
     -------
     res : OptimizeResult
@@ -304,7 +311,8 @@ def differential_evolution(func, bounds, args=(), strategy='best1bin',
                                      disp=disp, init=init, atol=atol,
                                      updating=updating,
                                      workers=workers,
-                                     constraints=constraints) as solver:
+                                     constraints=constraints,
+                                     x0=x0) as solver:
         ret = solver.solve()
 
     return ret
@@ -448,6 +456,10 @@ class DifferentialEvolutionSolver(object):
     constraints : {NonLinearConstraint, LinearConstraint, Bounds}
         Constraints on the solver, over and above those applied by the `bounds`
         kwd. Uses the approach by Lampinen.
+    x0 : None or array-like, optional
+        Provides an initial guess to the minimization. Once the population has
+        been initialized this vector replaces the first (best) member. This
+        replacement is done even if `init` is given an initial population.
     """
 
     # Dispatch of mutation strategy method (binomial or exponential).
@@ -473,7 +485,7 @@ def __init__(self, func, bounds, args=(),
                  tol=0.01, mutation=(0.5, 1), recombination=0.7, seed=None,
                  maxfun=np.inf, callback=None, disp=False, polish=True,
                  init='latinhypercube', atol=0, updating='immediate',
-                 workers=1, constraints=()):
+                 workers=1, constraints=(), x0=None):
 
         if strategy in self._binomial:
             self.mutation_func = getattr(self, self._binomial[strategy])
@@ -580,11 +592,15 @@ def __init__(self, func, bounds, args=(),
         else:
             self.init_population_array(init)
 
-        # infrastructure for constraints
-        # dummy parameter vector for preparing constraints, this is required so
-        # that the number of constraints is known.
-        x0 = self._scale_parameters(self.population[0])
+        if x0 is not None:
+            # scale to within unit interval and
+            # ensure parameters are within bounds.
+            x0_scaled = np.clip(
+                self._unscale_parameters(np.asarray(x0)), 0.0, 1.0
+            )
+            self.population[0] = x0_scaled
 
+        # infrastructure for constraints
         self.constraints = constraints
         self._wrapped_constraints = []
 
@@ -592,9 +608,13 @@ def __init__(self, func, bounds, args=(),
             # sequence of constraints, this will also deal with default
             # keyword parameter
             for c in constraints:
-                self._wrapped_constraints.append(_ConstraintWrapper(c, x0))
+                self._wrapped_constraints.append(
+                    _ConstraintWrapper(c, self.x)
+                )
         else:
-            self._wrapped_constraints = [_ConstraintWrapper(constraints, x0)]
+            self._wrapped_constraints = [
+                _ConstraintWrapper(constraints, self.x)
+            ]
 
         self.constraint_violation = np.zeros((self.num_population_members, 1))
         self.feasible = np.ones(self.num_population_members, bool)
diff --git a/scipy/optimize/tests/test__differential_evolution.py b/scipy/optimize/tests/test__differential_evolution.py
index f40a27b923e5..a826f427bb94 100644
--- a/scipy/optimize/tests/test__differential_evolution.py
+++ b/scipy/optimize/tests/test__differential_evolution.py
@@ -536,6 +536,15 @@ def test_population_initiation(self):
                       *(rosen, self.bounds),
                       **{'init': population})
 
+        # provide an initial solution
+        # bounds are [(0, 2), (0, 2)]
+        x0 = np.random.uniform(low=0.0, high=2.0, size=2)
+        solver = DifferentialEvolutionSolver(
+            rosen, self.bounds, x0=x0
+        )
+        # parameters are scaled to unit interval
+        assert_allclose(solver.population[0], x0 / 2.0)
+
     def test_infinite_objective_function(self):
         # Test that there are no problems if the objective function
         # returns inf on some runs


Code review comments:
diff_hunk: @@ -580,21 +592,29 @@ def __init__(self, func, bounds, args=(),
         else:
             self.init_population_array(init)
 
-        # infrastructure for constraints
-        # dummy parameter vector for preparing constraints, this is required so
-        # that the number of constraints is known.
-        x0 = self._scale_parameters(self.population[0])
+        if x0 is not None:
+            # scale to within unit interval and
+            # ensure parameters are within bounds.
+            x0_scaled = np.clip(
+                self._unscale_parameters(np.asarray(x0)), 0.0, 1.0
path: scipy/optimize/_differentialevolution.py
text: @user1:
Shouldn't `x0` ouside of specified `bounds` just raise a `ValueError`? Both are user-provided, so I'd consider the inconsistency to be user error. We typically prefer raising rather than silently changing the value.

@user1:
Also, can you add a test? There's still no test with `x0` being passed to the function that's public API.

@author:
> Shouldn't x0 ouside of specified bounds just raise a ValueError? 

If the user supplies an initial population to `init` the existing behaviour is to [do the same thing](https://github.com/scipy/scipy/blob/master/scipy/optimize/_differentialevolution.py#L678) as this PR. I think it definitely makes sense to check `x0` behaviour for the deterministic minimizers, but it's less important here. The solver success isn't degraded by a specific `x0`, so is a check worthwhile adding?

As the function exists it's quite fault tolerant, making the check is going to make it more likely to raise Errors that the user has to deal with, and even if they do fix the error I don't think the eventual solve trajectory will change all that much.