Repository: scipy/scipy
Pull request: MAINT: Update optimize.shgo
Reopening of previously closed issue #14013, also addresses the reviewer comments by @tupui discussed there (also now using from scipy._lib._util.MapWrapper to hopefully lower future maintenance).

Reference issues/fixes
===
Fixes for GitHub reported issues (with tests in scipy/optimize/tests/test__shgo.py):
- Closes gh-12114 (tested by `test_21_arg_tuple_sobol`). 
    - Random and Sobol sampling with constraint functions can now handle `args`.
- Closes gh-13547 (tested by `test_21_1_jac_true`).
   - Correctly handle objective functions that return the gradient alongside the objective value 
- Closes gh-13974 (tested by `test_10_finite_time`).
    - Fixed cases where shgo was not using all possible stopping criteria.
    - Fixed an issue where an extra iteration is run before stopping the routine.
    - The simplicial sampling now generates the correct number of specified points when `n` is specified with the non-default value.


Other fixes without unittests or linked issues:
- Fixed numerous typos in the documentation (https://github.com/Stefan-Endres/shgo/issues/29).
- Resolves some of the issues in gh-13813, but this is not fully closed (will have to be done in future documentation improvements).

Improved features
===

- Parallelization is now implemented with the `workers` argument.
- Added symmetry arguments which can improve computational performance up to a maximum of O(dim!) under ideal conditions.
- Class based usage (similar to https://github.com/scipy/scipy/pull/8414) to improve usability. Example calling `SHGO.iterate()` until a "satisfactory" solution is found when no other stopping criteria is available (commonly requested feature).
- Performance is improved due to various bug fixes. See the preliminary results below.
- By default shgo should now return the minimum value or its round-off (related to https://github.com/scipy/scipy/issues/10538)
- Improved speed of SHGO class initiation (serial).

Code streamlining
====

- Use a single caching library for all sampling methods. Reduces number of code lines in main `_shgo.py` file to improve maintainability and readability.
--  (The main reason this PR for _shgo_lib is unfortunately so large is because adopting a new cache turned out to be the only way to solve the issue with the `simplicial` sampling method not iterating correctly).
- Removed extra sobol point generations.
- Fixed anti-patterns:
-- `if a is 'teststring':`
-- [type(thing)](http://docs.quantifiedcode.com/python-anti-patterns/readability/do_not_compare_types_use_isinstance.html)


Benchmark profiles
===
As a preliminary test I reran the GO benchmarks for the `simplicial` method and compared it to the previous results quickly by overlaying a transparent window performance profiles over the previous results. As can be seen below the new `simplicial` method (darker blue solid line overlay) is noticeably faster than the previous results (lighter blue solid line).

![shgo_new bench](https://user-images.githubusercontent.com/8235072/117587965-e42a9d80-b120-11eb-9fbe-a8d725434fb4.png)

The increased performance largely due to fixing the issue that required the sampling method to refine the entire complex instead of a fixed number of points.

Still to do is to rerun all the results and new sampling methods together with the other GO methods.

TODO:
===

- [x] ResourceWarning causing unittest test_18_parallelization to fail.
- [x] Reimplement adding local minimizer results to the simplicial complex.
- [x] Fix importing error in scipy/stats/_hypotests.py#L7.
- [x] Possibly force using a fixed seed in `Sobol` and `Halton` sequences to improve deterministic performance.
- [ ] Rerun benchmarks.

Proposed reviewers
===

@tupui @andyfaff 

Problem statement:
scipy.optimize.shgo(): 'args' is incorrectly passed to constraint function when 'sobol' sampling is used
Hello!

I've discovered an issue in shgo when the 'args' argument is used together with constraint functions. 
It seems like the 'args' argument to shgo() is incorrectly passed on the constraint functions when sampling_method='sobol' is used.

I'm using scipy 1.4.1.

Code sample to demonstrate the issue:
```python
from scipy.optimize import shgo
res = shgo(
    lambda x, k: x[0]**k,
    bounds=[(0, 10)],
    args=(1,),
    constraints=({'type': 'ineq', 'fun': lambda x: x[0] - 1}),
    sampling_method='sobol',
)
```
This results in the following error:
```
/usr/lib64/python3.8/site-packages/scipy/optimize/_shgo.py in fun_ref(self)
   1442             if self.g_cons is not None:
   1443                 for g in self.g_cons:
-> 1444                     if g(self.C[i, :], *self.args) < 0.0:
   1445                         eval_f = False
   1446                         break  # Breaks the g loop

TypeError: <lambda>() takes 1 positional argument but 2 were given
```

If an additional argument is added to the constraint function (making it `lambda x, k: x[0] - 1`), I get a different error message,
```
/usr/lib64/python3.8/site-packages/scipy/optimize/slsqp.py in <listcomp>(.0)
    312     meq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    313               for c in cons['eq']]))
--> 314     mieq = sum(map(len, [atleast_1d(c['fun'](x, *c['args']))
    315                for c in cons['ineq']]))
    316     # m = The total number of constraints

TypeError: <lambda>() missing 1 required positional argument: 'k'
```

The issue can be worked around by having the constraint function take a variable number of arguments.

optimize.shgo: handle objective functions that return the gradient alongside the objective value
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

The shgo implementation seems to be unable to handle the gradient calculation option where the gradient is directly returned alongside the objective value by the objective function (as indicated by using `'jac': True`). 
I think the minimal working example below demonstrates the issue quite well.

It works fine when passing a separate Callable that calculates the gradient for `'jac'`

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```
from scipy.optimize import shgo

def func(x):
    return np.sum(np.power(x, 2)), 2 * x

shgo(
    func, 
    bounds=[[-1,1], [1,2], [-1,1], [1,2], [0,3]],  
    n=100, iters=5, 
    sampling_method="sobol", 
    minimizer_kwargs={'method': 'SLSQP', 'jac': True}
)
```

#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/home/david/anaconda3/envs/bachelor-project/lib/python3.8/site-packages/scipy/optimize/_shgo.py", line 421, in shgo
    shc.construct_complex()
  File "/home/david/anaconda3/envs/bachelor-project/lib/python3.8/site-packages/scipy/optimize/_shgo.py", line 726, in construct_complex
    self.iterate()
  File "/home/david/anaconda3/envs/bachelor-project/lib/python3.8/site-packages/scipy/optimize/_shgo.py", line 869, in iterate
    self.iterate_complex()
  File "/home/david/anaconda3/envs/bachelor-project/lib/python3.8/site-packages/scipy/optimize/_shgo.py", line 906, in iterate_delaunay
    self.sampled_surface(infty_cons_sampl=self.infty_cons_sampl)
  File "/home/david/anaconda3/envs/bachelor-project/lib/python3.8/site-packages/scipy/optimize/_shgo.py", line 1236, in sampled_surface
    self.fun_ref()
  File "/home/david/anaconda3/envs/bachelor-project/lib/python3.8/site-packages/scipy/optimize/_shgo.py", line 1447, in fun_ref
    self.F[i] = self.func(self.C[i, :], *self.args)
ValueError: setting an array element with a sequence.
```

#### Scipy/Numpy/Python version information:
```
scipy: 1.6.0 
numpy: 1.19.2 
sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)
```
-->

BUG: optimize.shgo: not using options
hi, 

I am using optimize.shgo to find the minimum of def f(x): x1,x2=x return M1(x1,x2)**2+M2(x1,x2)**2
where M1 and M2 are two complicated function which take a lot of time (9 secs) to evaluate. 
I notice that it takes a long time for optimize.shgo to finish (about 1hr) so I wan to use options dictionary to make it stop sooner then use the solution as initial guess for other minimizer but it seems that optimize.shgo doesn't care what options I give it. In fact even writing iters=0 it still runs!

can anyone here point me in the right direction?

Hints from issue discussion:
Thanks for the report @ahenriksson 

@Stefan-Endres would you mind having a look at this?
Hi @ahenriksson, thank you for the detailed report. I will write a quick fix a unittest case for this.


@Stefan-Endres 
Thank you for reporting this @cherrywoods and @pmla for the ping. This should be relatively easy to fix I just need to split the `func` argument output in `SHGO.__init__` when a bool instead of a vector function is supplied. Also need to update the documentation to match the new documentation of [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize).
Just to add that the error is also present when not using `sobol`.
That sounds good. Thanks for working on it!
@Stefan-Endres just checking in. Do you plan a PR to address some issues like these? Thanks.
Thank you for the reminder, I will try to finish the PR this weekend.
> Thank you for the reminder, I will try to finish the PR this weekend.

Thanks for the update 😃  
Hi. Thank you for reporting. Could you share us your code example reproducing the issue?.
Hello. Here is the input. M1 and M2 are defined as integrals of other functions so are integrated numerically.

def solve_moments(x):
    vw=x[0]
    Lh=x[1]
    return M1(vw,L)**2+M2(vw,L)**2
L=5/100
bnds = ((0.01, 0.99), (0.001*L, 10*L))
#opt_dict={"maxfev":2,"f_min":0,"f_tol":10e-1,"maxiter":1,"maxtime":1}
minimum=optimize.shgo(func=solve_moments,bounds=bnds,n=1, iters=0,options=None)


Here is the output: 

fun: 0.054116716357765826
    funl: array([0.05411672])
 message: 'Optimization terminated successfully.'
    nfev: 324
     nit: 2
   nlfev: 319
   nlhev: 0
   nljev: 55
 success: True
       x: array([0.76707782, 0.02817877])
      xl: array([[0.76707782, 0.02817877]])
me_bubble.guess

> I am using optimize.shgo to find the minimum of def f(x): x1,x2=x return M1(x1,x2)**2+M2(x1,x2)**2
> where M1 and M2 are two complicated function which take a lot of time (9 secs) to evaluate.
> I notice that it takes a long time for optimize.shgo to finish (about 1hr) so I wan to use options dictionary to make it stop sooner then use the solution as initial guess for other minimizer but it seems that optimize.shgo doesn't care what options I give it. In fact even writing iters=0 it still runs!

`iters` has a known bug. There are fixes for this on the develop version (for `sampling_method="sobol"`) but we still need to work on it. The way `n` works if a bit different with `sampling_method="simplicial"` than `sampling_method="sobol"` so the number do not match exactly.

@Stefan-Endres can you help here? I am not sure if this issue is not redundant with the ones we already have.

Hello, yes the `iters` bug has been fixed in the [upstream](https://github.com/Stefan-Endres/shgo) there is also a new stopping criteria called `maxtime` that can be passed to options which will stop the routine after a certain amount of time. I nearly done cleaning the code and documenting the cache library I aim to open the PR to scipy before the end of the week.
> Hello, yes the `iters` bug has been fixed in the [upstream](https://github.com/Stefan-Endres/shgo) there is also a new stopping criteria called `maxtime` that can be passed to options which will stop the routine after a certain amount of time. I nearly done cleaning the code and documenting the cache library I aim to open the PR to scipy before the end of the week.

Thanks for the prompt reply @Stefan-Endres. We will then link all these issues with your PR then, great! Let me know if I can help you with anything.
Hello again,

From the comments above it seems that I still cannot use the options {"maxfev":2,"f_min":0,"f_tol":10e-1,"maxiter":1,"maxtime":1}, right?
also the iter and n options only have an effec IF one uses sampling_method="sobol".

best regards,
Marco

> Hello again,
> 
> 
> 
> From the comments above it seems that I still cannot use the options {"maxfev":2,"f_min":0,"f_tol":10e-1,"maxiter":1,"maxtime":1}, right?
> 
> also the iter and n options only have an effec IF one uses sampling_method="sobol".
> 
> 
> 
> best regards,
> 
> Marco
> 
> 

I would have to check the code, but there are some checks in the tests regarding these options. As we said, the interface is still being worked. We will for sure see to check that with the upcoming changes. But I cannot guarantee that this would be ready for the next release as we should fix the release in a few weeks now.

Commit message: Update scipy/optimize/_shgo.py

Co-authored-by: Daniel Schmitz <40656107+dschmitz89@users.noreply.github.com>

Patch to review:
diff --git a/scipy/optimize/__init__.py b/scipy/optimize/__init__.py
index 22efc70e7020..6107f71d7a17 100644
--- a/scipy/optimize/__init__.py
+++ b/scipy/optimize/__init__.py
@@ -422,7 +422,7 @@
                            LinearConstraint,
                            Bounds)
 from ._hessian_update_strategy import HessianUpdateStrategy, BFGS, SR1
-from ._shgo import shgo
+from ._shgo import shgo, SHGO
 from ._dual_annealing import dual_annealing
 from ._qap import quadratic_assignment
 from ._direct_py import direct
diff --git a/scipy/optimize/_shgo.py b/scipy/optimize/_shgo.py
index 6556bad619d9..6ea4f7576978 100644
--- a/scipy/optimize/_shgo.py
+++ b/scipy/optimize/_shgo.py
@@ -1,24 +1,28 @@
-"""
-shgo: The simplicial homology global optimisation algorithm
-"""
-
-import numpy as np
+"""shgo: The simplicial homology global optimisation algorithm."""
+from collections import namedtuple
 import time
 import logging
 import warnings
+import sys
+
+import numpy as np
+
 from scipy import spatial
 from scipy.optimize import OptimizeResult, minimize, Bounds
+from scipy.optimize._optimize import MemoizeJac
 from scipy.optimize._constraints import new_bounds_to_old
-from ._optimize import _wrap_scalar_function
-from scipy.optimize._shgo_lib.triangulation import Complex
+from scipy._lib._util import _FunctionWrapper
 
+from scipy.optimize._shgo_lib._complex import Complex
 
 __all__ = ['shgo']
 
 
-def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
-         callback=None,
-         minimizer_kwargs=None, options=None, sampling_method='simplicial'):
+def shgo(
+    func, bounds, args=(), constraints=None, n=100, iters=1, callback=None,
+    minimizer_kwargs=None, options=None, sampling_method='simplicial', *,
+    workers=1
+):
     """
     Finds the global minimum of a function using SHG optimization.
 
@@ -77,12 +81,11 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
 
     n : int, optional
         Number of sampling points used in the construction of the simplicial
-        complex. Note that this argument is only used for ``sobol`` and other
-        arbitrary `sampling_methods`. In case of ``sobol``, it must be a
-        power of 2: ``n=2**m``, and the argument will automatically be
-        converted to the next higher power of 2. Default is 100 for
-        ``sampling_method='simplicial'`` and 128 for
-        ``sampling_method='sobol'``.
+        complex. For the default ``simplicial`` sampling method 2**dim + 1
+        sampling points are generated instead of the default `n=100`. For all
+        other specified values `n` sampling points are generated. For
+        ``sobol``, ``halton`` and other arbitrary `sampling_methods` `n=100` or
+        another speciefied number of sampling points are generated.
     iters : int, optional
         Number of iterations used in the construction of the simplicial
         complex. Default is 1.
@@ -94,7 +97,9 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
         ``scipy.optimize.minimize`` Some important options could be:
 
             * method : str
-                The minimization method, the default is ``SLSQP``.
+                The minimization method. If not given, chosen to be one of
+                BFGS, L-BFGS-B, SLSQP, depending on whether or not the
+                problem has constraints or bounds.
             * args : tuple
                 Extra arguments passed to the objective function (``func``) and
                 its derivatives (Jacobian, Hessian).
@@ -109,8 +114,8 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
         "(L)".
 
         Stopping criteria, the algorithm will terminate if any of the specified
-        criteria are met. However, the default algorithm does not require any to
-        be specified:
+        criteria are met. However, the default algorithm does not require any
+        to be specified:
 
         * maxfev : int (L)
             Maximum number of function evaluations in the feasible domain.
@@ -143,18 +148,36 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
 
         Objective function knowledge:
 
-        * symmetry : bool
-            Specify True if the objective function contains symmetric variables.
-            The search space (and therefore performance) is decreased by O(n!).
+        * symmetry : list or bool
+            Specify if the objective function contains symmetric variables.
+            The search space (and therefore performance) is decreased by up to
+            O(n!) times in the fully symmetric case. If `True` is specified
+            then all variables will be set symmetric to the first variable.
+            Default
+            is set to False.
+
+            E.g.  f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2
+
+            In this equation x_2 and x_3 are symmetric to x_1, while x_5 and
+            x_6 are symmetric to x_4, this can be specified to the solver as:
+
+            symmetry = [0,  # Variable 1
+                        0,  # symmetric to variable 1
+                        0,  # symmetric to variable 1
+                        3,  # Variable 4
+                        3,  # symmetric to variable 4
+                        3,  # symmetric to variable 4
+                        ]
 
         * jac : bool or callable, optional
             Jacobian (gradient) of objective function. Only for CG, BFGS,
             Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg. If ``jac`` is a
-            boolean and is True, ``fun`` is assumed to return the gradient along
-            with the objective function. If False, the gradient will be
+            boolean and is True, ``fun`` is assumed to return the gradient
+            along with the objective function. If False, the gradient will be
             estimated numerically. ``jac`` can also be a callable returning the
             gradient of the objective. In this case, it must accept the same
-            arguments as ``fun``. (Passed to `scipy.optimize.minmize` automatically)
+            arguments as ``fun``. (Passed to `scipy.optimize.minmize`
+            automatically)
 
         * hess, hessp : callable, optional
             Hessian (matrix of second-order derivatives) of objective function
@@ -171,13 +194,13 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
 
         * minimize_every_iter : bool
             If True then promising global sampling points will be passed to a
-            local minimization routine every iteration. If False then only the
-            final minimizer pool will be run. Defaults to False.
+            local minimization routine every iteration. If True then only the
+            final minimizer pool will be run. Defaults to True.
         * local_iter : int
             Only evaluate a few of the best minimizer pool candidates every
             iteration. If False all potential points are passed to the local
             minimization routine.
-        * infty_constraints: bool
+        * infty_constraints : bool
             If True then any sampling points generated which are outside will
             the feasible domain will be saved and given an objective function
             value of ``inf``. If False then these points will be discarded.
@@ -194,15 +217,19 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     sampling_method : str or function, optional
         Current built in sampling method options are ``halton``, ``sobol`` and
         ``simplicial``. The default ``simplicial`` provides
-        the theoretical guarantee of convergence to the global minimum in finite
-        time. ``halton`` and ``sobol`` method are faster in terms of sampling
-        point generation at the cost of the loss of
+        the theoretical guarantee of convergence to the global minimum in
+        finite time. ``halton`` and ``sobol`` method are faster in terms of
+        sampling point generation at the cost of the loss of
         guaranteed convergence. It is more appropriate for most "easier"
         problems where the convergence is relatively fast.
         User defined sampling functions must accept two arguments of ``n``
         sampling points of dimension ``dim`` per call and output an array of
         sampling points with shape `n x dim`.
 
+    workers : int  optional
+        Uses `multiprocessing.Pool <multiprocessing>`) to sample and run the
+        local serial minimizatons in parrallel.
+
     Returns
     -------
     res : OptimizeResult
@@ -259,11 +286,12 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     References
     ----------
     .. [1] Endres, SC, Sandrock, C, Focke, WW (2018) "A simplicial homology
-           algorithm for lipschitz optimisation", Journal of Global Optimization.
+           algorithm for lipschitz optimisation", Journal of Global
+           Optimization.
     .. [2] Joe, SW and Kuo, FY (2008) "Constructing Sobol' sequences with
            better  two-dimensional projections", SIAM J. Sci. Comput. 30,
            2635-2654.
-    .. [3] Hoch, W and Schittkowski, K (1981) "Test examples for nonlinear
+    .. [3] Hock, W and Schittkowski, K (1981) "Test examples for nonlinear
            programming codes", Lecture Notes in Economics and Mathematical
            Systems, 187. Springer-Verlag, New York.
            http://www.ai7.uni-bayreuth.de/test_problem_coll.pdf
@@ -275,7 +303,6 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     --------
     First consider the problem of minimizing the Rosenbrock function, `rosen`:
 
-    >>> import numpy as np
     >>> from scipy.optimize import rosen, shgo
     >>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
     >>> result = shgo(rosen, bounds)
@@ -297,6 +324,7 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     the capabilities of `shgo`.
     (https://en.wikipedia.org/wiki/Test_functions_for_optimization)
 
+    >>> import numpy as np
     >>> def eggholder(x):
     ...     return (-(x[1] + 47.0)
     ...             * np.sin(np.sqrt(abs(x[0]/2.0 + (x[1] + 47.0))))
@@ -345,9 +373,10 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     default of 1 to 3. Using ``simplicial`` this would have given us
     64 x 3 = 192 initial sampling points.
 
-    >>> result_2 = shgo(eggholder, bounds, n=64, iters=3, sampling_method='sobol')
+    >>> result_2 = shgo(eggholder,
+    ...                 bounds, n=64, iters=3, sampling_method='sobol')
     >>> len(result.xl), len(result_2.xl)
-    (12, 20)
+    (12, 23)
 
     Note the difference between, e.g., ``n=192, iters=1`` and ``n=64,
     iters=3``.
@@ -356,17 +385,18 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     sampling points for a total of 3 times.
 
     To demonstrate solving problems with non-linear constraints consider the
-    following example from Hock and Schittkowski problem 73 (cattle-feed) [3]_::
+    following example from Hock and Schittkowski problem 73 (cattle-feed)
+    [3]_::
 
         minimize: f = 24.55 * x_1 + 26.75 * x_2 + 39 * x_3 + 40.50 * x_4
 
-        subject to: 2.3 * x_1 + 5.6 * x_2 + 11.1 * x_3 + 1.3 * x_4 - 5     >= 0,
+        subject to: 2.3 * x_1 + 5.6 * x_2 + 11.1 * x_3 + 1.3 * x_4 - 5    >= 0,
 
                     12 * x_1 + 11.9 * x_2 + 41.8 * x_3 + 52.1 * x_4 - 21
                         -1.645 * sqrt(0.28 * x_1**2 + 0.19 * x_2**2 +
-                                      20.5 * x_3**2 + 0.62 * x_4**2)       >= 0,
+                                      20.5 * x_3**2 + 0.62 * x_4**2)      >= 0,
 
-                    x_1 + x_2 + x_3 + x_4 - 1                              == 0,
+                    x_1 + x_2 + x_3 + x_4 - 1                             == 0,
 
                     1 >= x_i >= 0 for all i
 
@@ -393,7 +423,7 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     ...         {'type': 'ineq', 'fun': g2},
     ...         {'type': 'eq', 'fun': h1})
     >>> bounds = [(0, 1.0),]*4
-    >>> res = shgo(f, bounds, iters=3, constraints=cons)
+    >>> res = shgo(f, bounds, n=150, constraints=cons)
     >>> res
      message: Optimization terminated successfully.
      success: True
@@ -401,8 +431,8 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
         funl: [ 2.989e+01]
            x: [ 6.355e-01  1.137e-13  3.127e-01  5.178e-02]
           xl: [[ 6.355e-01  1.137e-13  3.127e-01  5.178e-02]]
-         nit: 3
-        nfev: 114
+         nit: 1
+        nfev: 142
        nlfev: 35
        nljev: 5
        nlhev: 0
@@ -411,7 +441,6 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     (-5.062616992290714e-14, -2.9594104944408173e-12, 0.0)
 
     """
-
     # if necessary, convert bounds class to old bounds
     if isinstance(bounds, Bounds):
         bounds = new_bounds_to_old(bounds.lb, bounds.ub, len(bounds.lb))
@@ -420,14 +449,15 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
     shc = SHGO(func, bounds, args=args, constraints=constraints, n=n,
                iters=iters, callback=callback,
                minimizer_kwargs=minimizer_kwargs,
-               options=options, sampling_method=sampling_method)
+               options=options, sampling_method=sampling_method,
+               workers=workers)
 
     # Run the algorithm, process results and test success
-    shc.construct_complex()
+    shc.iterate_all()
 
     if not shc.break_routine:
         if shc.disp:
-            print("Successfully completed construction of complex.")
+            logging.info("Successfully completed construction of complex.")
 
     # Test post iterations success
     if len(shc.LMC.xl_maps) == 0:
@@ -440,6 +470,10 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
         shc.res.fun = shc.f_lowest
         shc.res.x = shc.x_lowest
         shc.res.nfev = shc.fn
+        shc.res.tnev = shc.n_sampled
+    else:
+        # Test that the optimal solutions do not violate any constraints
+        pass  # TODO
 
     # Confirm the routine ran successfully
     if not shc.break_routine:
@@ -453,19 +487,29 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
 class SHGO:
     def __init__(self, func, bounds, args=(), constraints=None, n=None,
                  iters=None, callback=None, minimizer_kwargs=None,
-                 options=None, sampling_method='sobol'):
-
+                 options=None, sampling_method='simplicial', workers=1):
         from scipy.stats import qmc
-
         # Input checks
         methods = ['halton', 'sobol', 'simplicial']
         if isinstance(sampling_method, str) and sampling_method not in methods:
             raise ValueError(("Unknown sampling_method specified."
                               " Valid methods: {}").format(', '.join(methods)))
 
+        # Split obj func if given with Jac
+        try:
+            if ((minimizer_kwargs['jac'] is True) and
+                    (not callable(minimizer_kwargs['jac']))):
+                self.func = MemoizeJac(func)
+                jac = self.func.derivative
+                minimizer_kwargs['jac'] = jac
+                func = self.func  # .fun
+            else:
+                self.func = func  # Normal definition of objective function
+        except (TypeError, KeyError):
+            self.func = func  # Normal definition of objective function
+
         # Initiate class
-        self._raw_func = func  # some methods pass args in (e.g. Complex)
-        _, self.func = _wrap_scalar_function(func, args)
+        self.func = _FunctionWrapper(func, args)
         self.bounds = bounds
         self.args = args
         self.callback = callback
@@ -489,6 +533,7 @@ def __init__(self, func, bounds, args=(), constraints=None, n=None,
 
         # Constraints
         # Process constraint dict sequence:
+        self.constraints = constraints
         if constraints is not None:
             self.min_cons = constraints
             self.g_cons = []
@@ -536,7 +581,7 @@ def __init__(self, func, bounds, args=(), constraints=None, n=None,
             self.init_options(options)
         else:  # Default settings:
             self.f_min_true = None
-            self.minimize_every_iter = False
+            self.minimize_every_iter = True
 
             # Algorithm limits
             self.maxiter = None
@@ -547,11 +592,11 @@ def __init__(self, func, bounds, args=(), constraints=None, n=None,
             self.minhgrd = None
 
             # Objective function knowledge
-            self.symmetry = False
+            self.symmetry = None
 
             # Algorithm functionality
-            self.local_iter = False
             self.infty_cons_sampl = True
+            self.local_iter = False
 
             # Feedback
             self.disp = False
@@ -597,53 +642,71 @@ def _restrict_to_keys(dictionary, goodkeys):
         self.stop_global = False  # Used in the stopping_criteria method
         self.break_routine = False  # Break the algorithm globally
         self.iters = iters  # Iterations to be ran
-        self.iters_done = 0  # Iterations to be ran
+        self.iters_done = 0  # Iterations completed
         self.n = n  # Sampling points per iteration
-        self.nc = n  # Sampling points to sample in current iteration
+        self.nc = 0  # n  # Sampling points to sample in current iteration
         self.n_prc = 0  # Processed points (used to track Delaunay iters)
-        self.n_sampled = 0  # To track number of sampling points already generated
+        self.n_sampled = 0  # To track no. of sampling points already generated
         self.fn = 0  # Number of feasible sampling points evaluations performed
         self.hgr = 0  # Homology group rank
+        # Initially attempt to build the triangulation incrementally:
+        self.qhull_incremental = True
 
         # Default settings if no sampling criteria.
+        if (self.n is None) and (self.iters is None) \
+                and (sampling_method == 'simplicial'):
+            self.n = 2 ** self.dim + 1
+            self.nc = 0  # self.n
         if self.iters is None:
             self.iters = 1
-        if self.n is None:
-            self.n = 100
-            if sampling_method == 'sobol':
-                self.n = 128
-            self.nc = self.n
+        if (self.n is None) and not (sampling_method == 'simplicial'):
+            self.n = self.n = 100
+            self.nc = 0  # self.n
+        if (self.n == 100) and (sampling_method == 'simplicial'):
+            self.n = 2 ** self.dim + 1
 
         if not ((self.maxiter is None) and (self.maxfev is None) and (
-                    self.maxev is None)
+                self.maxev is None)
                 and (self.minhgrd is None) and (self.f_min_true is None)):
             self.iters = None
 
         # Set complex construction mode based on a provided stopping criteria:
+        # Initialise sampling Complex and function cache
+        # Note that sfield_args=() since args are already wrapped in self.func
+        # using the_FunctionWrapper class.
+        self.HC = Complex(dim=self.dim, domain=self.bounds,
+                          sfield=self.func, sfield_args=(),
+                          symmetry=self.symmetry,
+                          constraints=self.constraints,
+                          workers=workers)
+
         # Choose complex constructor
         if sampling_method == 'simplicial':
             self.iterate_complex = self.iterate_hypercube
-            self.minimizers = self.simplex_minimizers
             self.sampling_method = sampling_method
 
         elif sampling_method in ['halton', 'sobol'] or \
                 not isinstance(sampling_method, str):
             self.iterate_complex = self.iterate_delaunay
-            self.minimizers = self.delaunay_complex_minimisers
             # Sampling method used
             if sampling_method in ['halton', 'sobol']:
                 if sampling_method == 'sobol':
                     self.n = int(2 ** np.ceil(np.log2(self.n)))
-                    self.nc = self.n
+                    # self.n #TODO: Should always be self.n, this is
+                    # unacceptable for shgo, check that nfev behaves as
+                    # expected.
+                    self.nc = 0
                     self.sampling_method = 'sobol'
                     self.qmc_engine = qmc.Sobol(d=self.dim, scramble=False,
-                                                seed=np.random.RandomState())
+                                                seed=0)
                 else:
                     self.sampling_method = 'halton'
                     self.qmc_engine = qmc.Halton(d=self.dim, scramble=True,
-                                                 seed=np.random.RandomState())
+                                                 seed=0)
+
                 def sampling_method(n, d):
                     return self.qmc_engine.random(n)
+
             else:
                 # A user defined sampling method:
                 self.sampling_method = 'custom'
@@ -696,7 +759,7 @@ def init_options(self, options):
                     self.minimizer_kwargs['options'].pop(opt))
 
         # Default settings:
-        self.minimize_every_iter = options.get('minimize_every_iter', False)
+        self.minimize_every_iter = options.get('minimize_every_iter', True)
 
         # Algorithm limits
         # Maximum number of iterations to perform.
@@ -719,20 +782,25 @@ def init_options(self, options):
         self.minhgrd = options.get('minhgrd', None)
 
         # Objective function knowledge
-        self.symmetry = 'symmetry' in options
-
+        self.symmetry = options.get('symmetry', False)
+        if self.symmetry:
+            self.symmetry = [0, ]*len(self.bounds)
+        else:
+            self.symmetry = None
         # Algorithm functionality
         # Only evaluate a few of the best candiates
         self.local_iter = options.get('local_iter', False)
-
         self.infty_cons_sampl = options.get('infty_constraints', True)
 
         # Feedback
         self.disp = options.get('disp', False)
 
+    def __exit__(self, *args):
+        return self.HC.V._mapwrapper.__exit__(*args)
+
     # Iteration properties
     # Main construction loop:
-    def construct_complex(self):
+    def iterate_all(self):
         """
         Construct for `iters` iterations.
 
@@ -743,7 +811,7 @@ def construct_complex(self):
 
         """
         if self.disp:
-            print('Splitting first generation')
+            logging.info('Splitting first generation')
 
         while not self.stop_global:
             if self.break_routine:
@@ -753,19 +821,25 @@ def construct_complex(self):
             self.stopping_criteria()
 
         # Build minimiser pool
-        # Final iteration only needed if pools weren't minimised every iteration
+        # Final iteration only needed if pools weren't minimised every
+        # iteration
         if not self.minimize_every_iter:
             if not self.break_routine:
                 self.find_minima()
 
-        self.res.nit = self.iters_done + 1
+        self.res.nit = self.iters_done  # + 1
+        self.fn = self.HC.V.nfev
 
     def find_minima(self):
         """
         Construct the minimizer pool, map the minimizers to local minima
         and sort the results into a global return object.
         """
+        if self.disp:
+            logging.info('Searching for minimizer pool...')
+
         self.minimizers()
+
         if len(self.X_min) != 0:
             # Minimize the pool of minimizers with local minimization methods
             # Note that if Options['local_iter'] is an `int` instead of default
@@ -780,50 +854,59 @@ def find_minima(self):
         else:
             self.find_lowest_vertex()
 
+        if self.disp:
+            logging.info(f"Minimiser pool = SHGO.X_min = {self.X_min}")
+
     def find_lowest_vertex(self):
         # Find the lowest objective function value on one of
         # the vertices of the simplicial complex
-        if self.sampling_method == 'simplicial':
-            self.f_lowest = np.inf
-            for x in self.HC.V.cache:
-                if self.HC.V[x].f < self.f_lowest:
-                    self.f_lowest = self.HC.V[x].f
-                    self.x_lowest = self.HC.V[x].x_a
-            if self.f_lowest == np.inf:  # no feasible point
-                self.f_lowest = None
-                self.x_lowest = None
-        else:
-            if self.fn == 0:
-                self.f_lowest = None
-                self.x_lowest = None
-            else:
-                self.f_I = np.argsort(self.F, axis=-1)
-                self.f_lowest = self.F[self.f_I[0]]
-                self.x_lowest = self.C[self.f_I[0]]
+        self.f_lowest = np.inf
+        for x in self.HC.V.cache:
+            if self.HC.V[x].f < self.f_lowest:
+                if self.disp:
+                    logging.info(f'self.HC.V[x].f = {self.HC.V[x].f}')
+                self.f_lowest = self.HC.V[x].f
+                self.x_lowest = self.HC.V[x].x_a
+        for lmc in self.LMC.cache:
+            if self.LMC[lmc].f_min < self.f_lowest:
+                self.f_lowest = self.LMC[lmc].f_min
+                self.x_lowest = self.LMC[lmc].x_l
+
+        if self.f_lowest == np.inf:  # no feasible point
+            self.f_lowest = None
+            self.x_lowest = None
 
     # Stopping criteria functions:
     def finite_iterations(self):
+        mi = min(x for x in [self.iters, self.maxiter] if x is not None)
+        logging.info(f'Iterations done = {self.iters_done} / {mi}')
+        logging.info(f'Iterations done = {self.iters_done} / {mi}')
         if self.iters is not None:
-            if self.iters_done >= (self.iters - 1):
+            if self.iters_done >= (self.iters):
                 self.stop_global = True
 
         if self.maxiter is not None:  # Stop for infeasible sampling
-            if self.iters_done >= (self.maxiter - 1):
+            if self.iters_done >= (self.maxiter):
                 self.stop_global = True
         return self.stop_global
 
     def finite_fev(self):
         # Finite function evals in the feasible domain
+        logging.info(f'Function evaluations done = {self.fn} / {self.maxfev}')
         if self.fn >= self.maxfev:
             self.stop_global = True
         return self.stop_global
 
     def finite_ev(self):
         # Finite evaluations including infeasible sampling points
+        logging.info(f'Sampling evaluations done = {self.n_sampled} '
+                     f'/ {self.maxev}')
         if self.n_sampled >= self.maxev:
             self.stop_global = True
 
     def finite_time(self):
+        logging.info(f'Time elapsed = {time.time() - self.init} '
+                     f'/ {self.maxtime}')
         if (time.time() - self.init) >= self.maxtime:
             self.stop_global = True
 
@@ -835,14 +918,17 @@ def finite_precision(self):
         and the tolerance with ``f_tol = options['f_tol']``
         """
         # If no minimizer has been found use the lowest sampling value
-        if len(self.LMC.xl_maps) == 0:
-            self.find_lowest_vertex()
+        self.find_lowest_vertex()
+        logging.info(f'Lowest function evaluation = {self.f_lowest}')
+        logging.info(f'Specified minimum = {self.f_min_true}')
+        # If no feasible point was return from test
+        if self.f_lowest is None:
+            return self.stop_global
 
         # Function to stop algorithm at specified percentage error:
-        if self.f_lowest == 0.0:
-            if self.f_min_true == 0.0:
-                if self.f_lowest <= self.f_tol:
-                    self.stop_global = True
+        if self.f_min_true == 0.0:
+            if self.f_lowest <= self.f_tol:
+                self.stop_global = True
         else:
             pe = (self.f_lowest - self.f_min_true) / abs(self.f_min_true)
             if self.f_lowest <= self.f_min_true:
@@ -859,6 +945,9 @@ def finite_precision(self):
         return self.stop_global
 
     def finite_homology_growth(self):
+        """
+        Stop the algorithm if homology group rank did not grow in iteration.
+        """
         if self.LMC.size == 0:
             return  # pass on no reason to stop yet.
         self.hgrd = self.LMC.size - self.hgr
@@ -866,6 +955,9 @@ def finite_homology_growth(self):
         self.hgr = self.LMC.size
         if self.hgrd <= self.minhgrd:
             self.stop_global = True
+
+        logging.info(f'Current homology growth = {self.hgrd} '
+                     f' (minimum growth = {self.minhgrd})')
         return self.stop_global
 
     def stopping_criteria(self):
@@ -890,6 +982,7 @@ def stopping_criteria(self):
             self.finite_precision()
         if self.minhgrd is not None:
             self.finite_homology_growth()
+        return self.stop_global
 
     def iterate(self):
         self.iterate_complex()
@@ -909,20 +1002,44 @@ def iterate_hypercube(self):
         Note: called with ``self.iterate_complex()`` after class initiation
         """
         # Iterate the complex
-        if self.n_sampled == 0:
-            # Initial triangulation of the hyper-rectangle. Note that
-            # we use `self.raw_func` as `self.func` is a *wrapped* function
-            # that already takes the original function arguments into
-            # account.
-            self.HC = Complex(self.dim, self._raw_func, self.args,
-                              self.symmetry, self.bounds, self.g_cons,
-                              self.g_args)
+        if self.disp:
+            logging.info('Constructing and refining simplicial complex graph '
+                         'structure')
+        if self.n is None:
+            self.HC.refine_all()
+            self.n_sampled = self.HC.V.size()  # nevs counted
         else:
-            self.HC.split_generation()
+            self.HC.refine(self.n)
+            self.n_sampled += self.n
+
+        if self.disp:
+            logging.info('Triangulation completed, evaluating all contraints '
+                         'and objective function values.')
+
+        # Readd minimisers to complex
+        if len(self.LMC.xl_maps) > 0:
+            for xl in self.LMC.cache:
+                v = self.HC.V[xl]
+                v_near = v.star()
+                for v in v.nn:
+                    v_near = v_near.union(v.nn)
+                # Reconnect vertices to complex
+                # if self.HC.connect_vertex_non_symm(tuple(self.LMC[xl].x_l),
+                #                                   near=v_near):
+                #    continue
+                # else:
+                    # If failure to find in v_near, then search all vertices
+                    # (very expensive operation:
+                #    self.HC.connect_vertex_non_symm(tuple(self.LMC[xl].x_l)
+                #                                    )
+
+        # Evaluate all constraints and functions
+        self.HC.V.process_pools()
+        if self.disp:
+            logging.info('Evaluations completed.')
 
         # feasible sampling points counted by the triangulation.py routines
         self.fn = self.HC.V.nfev
-        self.n_sampled = self.HC.V.size  # nevs counted in triangulation.py
         return
 
     def iterate_delaunay(self):
@@ -931,23 +1048,76 @@ def iterate_delaunay(self):
 
         Note: called with ``self.iterate_complex()`` after class initiation
         """
-        self.sampled_surface(infty_cons_sampl=self.infty_cons_sampl)
         self.nc += self.n
-        self.n_sampled = self.nc
+        self.sampled_surface(infty_cons_sampl=self.infty_cons_sampl)
+
+        # Add sampled points to a triangulation, construct self.Tri
+        if self.disp:
+            logging.info(f'self.n = {self.n}')
+            logging.info(f'self.nc = {self.nc}')
+            logging.info('Constructing and refining simplicial complex graph '
+                         'structure from sampling points.')
+
+        if self.dim < 2:
+            self.Ind_sorted = np.argsort(self.C, axis=0)
+            self.Ind_sorted = self.Ind_sorted.flatten()
+            tris = []
+            for ind, ind_s in enumerate(self.Ind_sorted):
+                if ind > 0:
+                    tris.append(self.Ind_sorted[ind - 1:ind + 1])
+
+            tris = np.array(tris)
+            # Store 1D triangulation:
+            self.Tri = namedtuple('Tri', ['points', 'simplices'])(self.C, tris)
+            self.points = {}
+        else:
+            if self.C.shape[0] > self.dim + 1:  # Ensure a simplex can be built
+                self.delaunay_triangulation(n_prc=self.n_prc)
+            self.n_prc = self.C.shape[0]
+
+        if self.disp:
+            logging.info('Triangulation completed, evaluating all '
+                         'contraints and objective function values.')
+
+        if hasattr(self, 'Tri'):
+            self.HC.vf_to_vv(self.Tri.points, self.Tri.simplices)
+
+        # Process all pools
+        # Evaluate all constraints and functions
+        if self.disp:
+            logging.info('Triangulation completed, evaluating all contraints '
+                         'and objective function values.')
+
+        # Evaluate all constraints and functions
+        self.HC.V.process_pools()
+        if self.disp:
+            logging.info('Evaluations completed.')
+
+        # feasible sampling points counted by the triangulation.py routines
+        self.fn = self.HC.V.nfev
+        self.n_sampled = self.nc  # nevs counted in triangulation
+        return
 
     # Hypercube minimizers
-    def simplex_minimizers(self):
+    def minimizers(self):
         """
         Returns the indexes of all minimizers
         """
         self.minimizer_pool = []
         # Note: Can implement parallelization here
         for x in self.HC.V.cache:
+            in_LMC = False
+            if len(self.LMC.xl_maps) > 0:
+                for xlmi in self.LMC.xl_maps:
+                    if np.all(np.array(x) == np.array(xlmi)):
+                        in_LMC = True
+            if in_LMC:
+                continue
+
             if self.HC.V[x].minimiser():
                 if self.disp:
                     logging.info('=' * 60)
-                    logging.info(
-                        f'v.x = {self.HC.V[x].x_a} is minimizer')
+                    logging.info(f'v.x = {self.HC.V[x].x_a} is minimizer')
                     logging.info(f'v.f = {self.HC.V[x].f} is minimizer')
                     logging.info('=' * 30)
 
@@ -961,7 +1131,6 @@ def simplex_minimizers(self):
                         logging.info(f'x = {vn.x} || f = {vn.f}')
 
                     logging.info('=' * 60)
-
         self.minimizer_pool_F = []
         self.X_min = []
         # normalized tuple in the Vertex cache
@@ -980,12 +1149,12 @@ def simplex_minimizers(self):
 
         return self.X_min
 
-    # Local minimization
-    # Minimizer pool processing
+    # Local minimisation
+    # Minimiser pool processing
     def minimise_pool(self, force_iter=False):
         """
         This processing method can optionally minimise only the best candidate
-        solutions in the minimizer pool
+        solutions in the minimiser pool
 
         Parameters
         ----------
@@ -1002,25 +1171,14 @@ def minimise_pool(self, force_iter=False):
         # Trim minimized point from current minimizer set
         self.trim_min_pool(0)
 
-        # Force processing to only
-        if force_iter:
-            self.local_iter = force_iter
-
         while not self.stop_l_iter:
             # Global stopping criteria:
-            if self.f_min_true is not None:
-                if (lres_f_min.fun - self.f_min_true) / abs(
-                        self.f_min_true) <= self.f_tol:
-                    self.stop_l_iter = True
-                    break
+            self.stopping_criteria()
+
             # Note first iteration is outside loop:
-            if self.local_iter is not None:
-                if self.disp:
-                    logging.info(
-                        'SHGO.iters in function minimise_pool = {}'.format(
-                            self.local_iter))
-                self.local_iter -= 1
-                if self.local_iter == 0:
+            if force_iter:
+                force_iter -= 1
+                if force_iter == 0:
                     self.stop_l_iter = True
                     break
 
@@ -1090,7 +1248,8 @@ def construct_lcb_simplicial(self, v_min):
         Returns
         -------
         cbounds : list of lists
-            List of size dimension with length-2 list of bounds for each dimension
+            List of size dimension with length-2 list of bounds for each
+            dimension.
 
         """
         cbounds = [[x_b_i[0], x_b_i[1]] for x_b_i in self.bounds]
@@ -1123,7 +1282,8 @@ def construct_lcb_delaunay(self, v_min, ind=None):
         Returns
         -------
         cbounds : list of lists
-            List of size dimension with length-2 list of bounds for each dimension
+            List of size dimension with length-2 list of bounds for each
+            dimension.
         """
         cbounds = [[x_b_i[0], x_b_i[1]] for x_b_i in self.bounds]
 
@@ -1151,43 +1311,43 @@ def minimize(self, x_min, ind=None):
             logging.info(f'Vertex minimiser maps = {self.LMC.v_maps}')
 
         if self.LMC[x_min].lres is not None:
+            logging.info(f'Found self.LMC[x_min].lres = '
+                         f'{self.LMC[x_min].lres}')
             return self.LMC[x_min].lres
 
-        # TODO: Check discarded bound rules
-
         if self.callback is not None:
-            print('Callback for '
+            logging.info('Callback for '
                   'minimizer starting at {}:'.format(x_min))
 
         if self.disp:
-            print('Starting '
+            logging.info('Starting '
                   'minimization at {}...'.format(x_min))
 
         if self.sampling_method == 'simplicial':
             x_min_t = tuple(x_min)
             # Find the normalized tuple in the Vertex cache:
             x_min_t_norm = self.X_min_cache[tuple(x_min_t)]
-
             x_min_t_norm = tuple(x_min_t_norm)
-
             g_bounds = self.construct_lcb_simplicial(self.HC.V[x_min_t_norm])
             if 'bounds' in self.min_solver_args:
                 self.minimizer_kwargs['bounds'] = g_bounds
+                logging.info(self.minimizer_kwargs['bounds'])
 
         else:
             g_bounds = self.construct_lcb_delaunay(x_min, ind=ind)
             if 'bounds' in self.min_solver_args:
                 self.minimizer_kwargs['bounds'] = g_bounds
+                logging.info(self.minimizer_kwargs['bounds'])
 
         if self.disp and 'bounds' in self.minimizer_kwargs:
-            print('bounds in kwarg:')
-            print(self.minimizer_kwargs['bounds'])
+            logging.info('bounds in kwarg:')
+            logging.info(self.minimizer_kwargs['bounds'])
 
         # Local minimization using scipy.optimize.minimize:
         lres = minimize(self.func, x_min, **self.minimizer_kwargs)
 
         if self.disp:
-            print(f'lres = {lres}')
+            logging.info(f'lres = {lres}')
 
         # Local function evals for all minimizers
         self.res.nlfev += lres.nfev
@@ -1238,21 +1398,21 @@ def sampled_surface(self, infty_cons_sampl=False):
         There are 2 modes, if ``infty_cons_sampl`` is True then the sampled
         points that are generated outside the feasible domain will be
         assigned an ``inf`` value in accordance with SHGO rules.
-        This guarantees convergence and usually requires less objective function
-        evaluations at the computational costs of more Delaunay triangulation
-        points.
-
-        If ``infty_cons_sampl`` is False, then the infeasible points are discarded
-        and only a subspace of the sampled points are used. This comes at the
-        cost of the loss of guaranteed convergence and usually requires more
-        objective function evaluations.
+        This guarantees convergence and usually requires less objective
+        function evaluations at the computational costs of more Delaunay
+        triangulation points.
+
+        If ``infty_cons_sampl`` is False, then the infeasible points are
+        discarded and only a subspace of the sampled points are used. This
+        comes at the cost of the loss of guaranteed convergence and usually
+        requires more objective function evaluations.
         """
         # Generate sampling points
         if self.disp:
-            print('Generating sampling points')
+            logging.info('Generating sampling points')
         self.sampling(self.nc, self.dim)
-        self.n = self.nc
-
+        if len(self.LMC.xl_maps) > 0:
+            self.C = np.vstack((self.C, np.array(self.LMC.xl_maps)))
         if not infty_cons_sampl:
             # Find subspace of feasible points
             if self.g_cons is not None:
@@ -1262,52 +1422,8 @@ def sampled_surface(self, infty_cons_sampl=False):
         self.sorted_samples()
 
         # Find objective function references
-        self.fun_ref()
-
         self.n_sampled = self.nc
 
-    def delaunay_complex_minimisers(self):
-        # Construct complex minimizers on the current sampling set.
-        # if self.fn >= (self.dim + 1):
-        if self.fn >= (self.dim + 2):
-            # TODO: Check on strange Qhull error where the number of vertices
-            # required for an initial simplex is higher than n + 1?
-            if self.dim < 2:  # Scalar objective functions
-                if self.disp:
-                    print('Constructing 1-D minimizer pool')
-
-                self.ax_subspace()
-                self.surface_topo_ref()
-                self.minimizers_1D()
-
-            else:  # Multivariate functions.
-                if self.disp:
-                    print('Constructing Gabrial graph and minimizer pool')
-
-                if self.iters == 1:
-                    self.delaunay_triangulation(grow=False)
-                else:
-                    self.delaunay_triangulation(grow=True, n_prc=self.n_prc)
-                    self.n_prc = self.C.shape[0]
-
-                if self.disp:
-                    print('Triangulation completed, building minimizer pool')
-
-                self.delaunay_minimizers()
-
-            if self.disp:
-                logging.info(
-                    f"Minimizer pool = SHGO.X_min = {self.X_min}")
-        else:
-            if self.disp:
-                print(
-                    'Not enough sampling points found in the feasible domain.')
-            self.minimizer_pool = [None]
-            try:
-                self.X_min
-            except AttributeError:
-                self.X_min = []
-
     def sampling_custom(self, n, dim):
         """
         Generates uniform sampling points in a hypercube and scales the points
@@ -1315,7 +1431,10 @@ def sampling_custom(self, n, dim):
         """
         # Generate sampling points.
         # Generate uniform sample points in [0, 1]^m \subset R^m
-        self.C = self.sampling_function(n, dim)
+        if self.n_sampled == 0:
+            self.C = self.sampling_function(n, dim)
+        else:
+            self.C = self.sampling_function(n, dim)
         # Distribute over bounds
         for i in range(len(self.bounds)):
             self.C[:, i] = (self.C[:, i] *
@@ -1334,7 +1453,7 @@ def sampling_subspace(self):
                                     + 'size.')
                 # sampling correctly for both 1-D and >1-D cases
                 if self.disp:
-                    print(self.res.message)
+                    logging.info(self.res.message)
 
     def sorted_samples(self):  # Validated
         """Find indexes of the sorted sampling points"""
@@ -1342,193 +1461,34 @@ def sorted_samples(self):  # Validated
         self.Xs = self.C[self.Ind_sorted]
         return self.Ind_sorted, self.Xs
 
-    def ax_subspace(self):  # Validated
-        """
-        Finds the subspace vectors along each component axis.
-        """
-        self.Ci = []
-        self.Xs_i = []
-        self.Ii = []
-        for i in range(self.dim):
-            self.Ci.append(self.C[:, i])
-            self.Ii.append(self.Ind_sorted[:, i])
-            self.Xs_i.append(self.Xs[:, i])
-
-    def fun_ref(self):
-        """
-        Find the objective function output reference table
-        """
-        # TODO: Replace with cached wrapper
-
-        # Note: This process can be pooled easily
-        # Obj. function returns to be used as reference table.:
-        f_cache_bool = False
-        if self.fn > 0:  # Store old function evaluations
-            Ftemp = self.F
-            fn_old = self.fn
-            f_cache_bool = True
-
-        self.F = np.zeros(np.shape(self.C)[0])
-        # NOTE: It might be easier to replace this with a cached
-        #      objective function
-        for i in range(self.fn, np.shape(self.C)[0]):
-            eval_f = True
-            if self.g_cons is not None:
-                for g in self.g_cons:
-                    if g(self.C[i, :], *self.args) < 0.0:
-                        eval_f = False
-                        break  # Breaks the g loop
-
-            if eval_f:
-                self.F[i] = self.func(self.C[i, :])
-                self.fn += 1
-            elif self.infty_cons_sampl:
-                self.F[i] = np.inf
-                self.fn += 1
-        if f_cache_bool:
-            if fn_old > 0:  # Restore saved function evaluations
-                self.F[0:fn_old] = Ftemp
-
-        return self.F
-
-    def surface_topo_ref(self):  # Validated
-        """
-        Find the BD and FD finite differences along each component vector.
-        """
-        # Replace numpy inf, -inf and nan objects with floating point numbers
-        # nan --> float
-        self.F[np.isnan(self.F)] = np.inf
-        # inf, -inf  --> floats
-        self.F = np.nan_to_num(self.F)
-
-        self.Ft = self.F[self.Ind_sorted]
-        self.Ftp = np.diff(self.Ft, axis=0)  # FD
-        self.Ftm = np.diff(self.Ft[::-1], axis=0)[::-1]  # BD
-
-    def sample_topo(self, ind):
-        # Find the position of the sample in the component axial directions
-        self.Xi_ind_pos = []
-        self.Xi_ind_topo_i = []
-
-        for i in range(self.dim):
-            for x, I_ind in zip(self.Ii[i], range(len(self.Ii[i]))):
-                if x == ind:
-                    self.Xi_ind_pos.append(I_ind)
-
-            # Use the topo reference tables to find if point is a minimizer on
-            # the current axis
-
-            # First check if index is on the boundary of the sampling points:
-            if self.Xi_ind_pos[i] == 0:
-                # if boundary is in basin
-                self.Xi_ind_topo_i.append(self.Ftp[:, i][0] > 0)
-
-            elif self.Xi_ind_pos[i] == self.fn - 1:
-                # Largest value at sample size
-                self.Xi_ind_topo_i.append(self.Ftp[:, i][self.fn - 2] < 0)
-
-            # Find axial reference for other points
-            else:
-                Xi_ind_top_p = self.Ftp[:, i][self.Xi_ind_pos[i]] > 0
-                Xi_ind_top_m = self.Ftm[:, i][self.Xi_ind_pos[i] - 1] > 0
-                self.Xi_ind_topo_i.append(Xi_ind_top_p and Xi_ind_top_m)
-
-        if np.array(self.Xi_ind_topo_i).all():
-            self.Xi_ind_topo = True
+    def delaunay_triangulation(self, n_prc=0):
+        if hasattr(self, 'Tri') and self.qhull_incremental:
+            # TODO: Uncertain if n_prc needs to add len(self.LMC.xl_maps)
+            # in self.sampled_surface
+            self.Tri.add_points(self.C[n_prc:, :])
         else:
-            self.Xi_ind_topo = False
-        self.Xi_ind_topo = np.array(self.Xi_ind_topo_i).all()
-
-        return self.Xi_ind_topo
-
-    def minimizers_1D(self):
-        """
-        Returns the indices of all minimizers
-        """
-        self.minimizer_pool = []
-        # Note: Can implement parallelization here
-        for ind in range(self.fn):
-            min_bool = self.sample_topo(ind)
-            if min_bool:
-                self.minimizer_pool.append(ind)
-
-        self.minimizer_pool_F = self.F[self.minimizer_pool]
-
-        # Sort to find minimum func value in min_pool
-        self.sort_min_pool()
-        if not len(self.minimizer_pool) == 0:
-            self.X_min = self.C[self.minimizer_pool]
-            # If function is called again and pool is found unbreak:
-        else:
-            self.X_min = []
-
-        return self.X_min
-
-    def delaunay_triangulation(self, grow=False, n_prc=0):
-        if not grow:
-            self.Tri = spatial.Delaunay(self.C)
-        else:
-            if hasattr(self, 'Tri'):
-                self.Tri.add_points(self.C[n_prc:, :])
-            else:
-                self.Tri = spatial.Delaunay(self.C, incremental=True)
+            try:
+                self.Tri = spatial.Delaunay(self.C,
+                                            incremental=self.qhull_incremental,
+                                            )
+            except spatial.QhullError:
+                if str(sys.exc_info()[1])[:6] == 'QH6239':
+                    logging.warning('QH6239 Qhull precision error detected, '
+                                    'this usually occurs when no bounds are '
+                                    'specified, Qhull can only run with '
+                                    'handling cocircular/cospherical points'
+                                    ' and in this case incremental mode is '
+                                    'switched off. The performance of shgo '
+                                    'will be reduced in this mode.')
+                    self.qhull_incremental = False
+                    self.Tri = spatial.Delaunay(self.C,
+                                                incremental=
+                                                self.qhull_incremental)
+                else:
+                    raise
 
         return self.Tri
 
-    @staticmethod
-    def find_neighbors_delaunay(pindex, triang):
-        """
-        Returns the indices of points connected to ``pindex`` on the Gabriel
-        chain subgraph of the Delaunay triangulation.
-        """
-        return triang.vertex_neighbor_vertices[1][
-               triang.vertex_neighbor_vertices[0][pindex]:
-               triang.vertex_neighbor_vertices[0][pindex + 1]]
-
-    def sample_delaunay_topo(self, ind):
-        self.Xi_ind_topo_i = []
-
-        # Find the position of the sample in the component Gabrial chain
-        G_ind = self.find_neighbors_delaunay(ind, self.Tri)
-
-        # Find finite deference between each point
-        for g_i in G_ind:
-            rel_topo_bool = self.F[ind] < self.F[g_i]
-            self.Xi_ind_topo_i.append(rel_topo_bool)
-
-        # Check if minimizer
-        self.Xi_ind_topo = np.array(self.Xi_ind_topo_i).all()
-
-        return self.Xi_ind_topo
-
-    def delaunay_minimizers(self):
-        """
-        Returns the indices of all minimizers
-        """
-        self.minimizer_pool = []
-        # Note: Can easily be parralized
-        if self.disp:
-            logging.info(f'self.fn = {self.fn}')
-            logging.info(f'self.nc = {self.nc}')
-            logging.info('np.shape(self.C)'
-                         ' = {}'.format(np.shape(self.C)))
-        for ind in range(self.fn):
-            min_bool = self.sample_delaunay_topo(ind)
-            if min_bool:
-                self.minimizer_pool.append(ind)
-
-        self.minimizer_pool_F = self.F[self.minimizer_pool]
-
-        # Sort to find minimum func value in min_pool
-        self.sort_min_pool()
-        if self.disp:
-            logging.info(f'self.minimizer_pool = {self.minimizer_pool}')
-        if not len(self.minimizer_pool) == 0:
-            self.X_min = self.C[self.minimizer_pool]
-        else:
-            self.X_min = []  # Empty pool breaks main routine
-        return self.X_min
-
 
 class LMap:
     def __init__(self, v):
@@ -1546,12 +1506,16 @@ def __init__(self):
         # Lists for search queries
         self.v_maps = []
         self.xl_maps = []
+        self.xl_maps_set = set()
         self.f_maps = []
         self.lbound_maps = []
         self.size = 0
 
     def __getitem__(self, v):
-        v = np.ndarray.tolist(v)
+        try:
+            v = np.ndarray.tolist(v)
+        except TypeError:
+            pass
         v = tuple(v)
         try:
             return self.cache[v]
@@ -1575,6 +1539,7 @@ def add_res(self, v, lres, bounds=None):
         # Cache lists for search queries
         self.v_maps.append(v)
         self.xl_maps.append(lres.x)
+        self.xl_maps_set.add(tuple(lres.x))
         self.f_maps.append(lres.fun)
         self.lbound_maps.append(bounds)
 
diff --git a/scipy/optimize/_shgo_lib/_complex.py b/scipy/optimize/_shgo_lib/_complex.py
new file mode 100644
index 000000000000..2c8344c2d3c2
--- /dev/null
+++ b/scipy/optimize/_shgo_lib/_complex.py
@@ -0,0 +1,1226 @@
+"""Base classes for low memory simplicial complex structures."""
+import copy
+import logging
+import itertools
+import decimal
+from functools import lru_cache
+
+import numpy
+
+from ._vertex import (VertexCacheField, VertexCacheIndex)
+
+
+class Complex:
+    """
+    Base class for a simplicial complex described as a cache of vertices
+    together with their connections.
+
+    Important methods:
+        Domain triangulation:
+                Complex.triangulate, Complex.split_generation
+        Triangulating arbitrary points (must be traingulable,
+            may exist outside domain):
+                Complex.triangulate(sample_set)
+        Converting another simplicial complex structure data type to the
+            structure used in Complex (ex. OBJ wavefront)
+                Complex.convert(datatype, data)
+
+    Important objects:
+        HC.V: The cache of vertices and their connection
+        HC.H: Storage structure of all vertex groups
+
+    Parameters
+    ----------
+    dim : int
+        Spatial dimensionality of the complex R^dim
+    domain : list of tuples, optional
+        The bounds [x_l, x_u]^dim of the hyperrectangle space
+        ex. The default domain is the hyperrectangle [0, 1]^dim
+        Note: The domain must be convex, non-convex spaces can be cut
+              away from this domain using the non-linear
+              g_cons functions to define any arbitrary domain
+              (these domains may also be disconnected from each other)
+    sfield :
+        A scalar function defined in the associated domain f: R^dim --> R
+    sfield_args : tuple
+        Additional arguments to be passed to `sfield`
+    vfield :
+        A scalar function defined in the associated domain
+                       f: R^dim --> R^m
+                   (for example a gradient function of the scalar field)
+    vfield_args : tuple
+        Additional arguments to be passed to vfield
+    symmetry : None or list
+            Specify if the objective function contains symmetric variables.
+            The search space (and therefore performance) is decreased by up to
+            O(n!) times in the fully symmetric case.
+
+            E.g.  f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2
+
+            In this equation x_2 and x_3 are symmetric to x_1, while x_5 and
+             x_6 are symmetric to x_4, this can be specified to the solver as:
+
+            symmetry = [0,  # Variable 1
+                        0,  # symmetric to variable 1
+                        0,  # symmetric to variable 1
+                        3,  # Variable 4
+                        3,  # symmetric to variable 4
+                        3,  # symmetric to variable 4
+                        ]
+
+    constraints : dict or sequence of dict, optional
+        Constraints definition.
+        Function(s) ``R**n`` in the form::
+
+            g(x) <= 0 applied as g : R^n -> R^m
+            h(x) == 0 applied as h : R^n -> R^p
+
+        Each constraint is defined in a dictionary with fields:
+
+            type : str
+                Constraint type: 'eq' for equality, 'ineq' for inequality.
+            fun : callable
+                The function defining the constraint.
+            jac : callable, optional
+                The Jacobian of `fun` (only for SLSQP).
+            args : sequence, optional
+                Extra arguments to be passed to the function and Jacobian.
+
+        Equality constraint means that the constraint function result is to
+        be zero whereas inequality means that it is to be
+        non-negative.constraints : dict or sequence of dict, optional
+        Constraints definition.
+        Function(s) ``R**n`` in the form::
+
+            g(x) <= 0 applied as g : R^n -> R^m
+            h(x) == 0 applied as h : R^n -> R^p
+
+        Each constraint is defined in a dictionary with fields:
+
+            type : str
+                Constraint type: 'eq' for equality, 'ineq' for inequality.
+            fun : callable
+                The function defining the constraint.
+            jac : callable, optional
+                The Jacobian of `fun` (unused).
+            args : sequence, optional
+                Extra arguments to be passed to the function and Jacobian.
+
+        Equality constraint means that the constraint function result is to
+        be zero whereas inequality means that it is to be non-negative.
+
+    workers : int  optional
+        Uses `multiprocessing.Pool <multiprocessing>`) to compute the field
+         functions in parrallel.
+    """
+    def __init__(self, dim, domain=None, sfield=None, sfield_args=(),
+                 symmetry=None, constraints=None, workers=1):
+        self.dim = dim
+
+        # Domains
+        self.domain = domain
+        if domain is None:
+            self.bounds = [(float(0), float(1.0)), ] * dim
+        else:
+            self.bounds = domain
+        self.symmetry = symmetry
+        #      here in init to avoid if checks
+
+        # Field functions
+        self.sfield = sfield
+        self.sfield_args = sfield_args
+
+        # Process constraints
+        # Constraints
+        # Process constraint dict sequence:
+        if constraints is not None:
+            self.min_cons = constraints
+            self.g_cons = []
+            self.g_args = []
+            if (type(constraints) is not tuple) and (type(constraints)
+                                                     is not list):
+                constraints = (constraints,)
+
+            for cons in constraints:
+                if cons['type'] == 'ineq':
+                    self.g_cons.append(cons['fun'])
+                    try:
+                        self.g_args.append(cons['args'])
+                    except KeyError:
+                        self.g_args.append(())
+            self.g_cons = tuple(self.g_cons)
+            self.g_args = tuple(self.g_args)
+        else:
+            self.g_cons = None
+            self.g_args = None
+
+        # Homology properties
+        self.gen = 0
+        self.perm_cycle = 0
+
+        # Every cell is stored in a list of its generation,
+        # ex. the initial cell is stored in self.H[0]
+        # 1st get new cells are stored in self.H[1] etc.
+        # When a cell is sub-generated it is removed from this list
+
+        self.H = []  # Storage structure of vertex groups
+
+        # Cache of all vertices
+        if (sfield is not None) or (self.g_cons is not None):
+            # Initiate a vertex cache and an associated field cache, note that
+            # the field case is always initiated inside the vertex cache if an
+            # associated field scalar field is defined:
+            if sfield is not None:
+                self.V = VertexCacheField(field=sfield, field_args=sfield_args,
+                                          g_cons=self.g_cons,
+                                          g_cons_args=self.g_args,
+                                          workers=workers)
+            elif self.g_cons is not None:
+                self.V = VertexCacheField(field=sfield, field_args=sfield_args,
+                                          g_cons=self.g_cons,
+                                          g_cons_args=self.g_args,
+                                          workers=workers)
+        else:
+            self.V = VertexCacheIndex()
+
+        self.V_non_symm = []  # List of non-symmetric vertices
+
+    def __call__(self):
+        return self.H
+
+    # %% Triangulation methods
+    def cyclic_product(self, bounds, origin, supremum, centroid=True):
+        """Generate initial triangulation using cyclic product"""
+        # Define current hyperrectangle
+        vot = tuple(origin)
+        vut = tuple(supremum)  # Hyperrectangle supremum
+        self.V[vot]
+        vo = self.V[vot]
+        yield vo.x
+        self.V[vut].connect(self.V[vot])
+        yield vut
+        # Cyclic group approach with second x_l --- x_u operation.
+
+        # These containers store the "lower" and "upper" vertices
+        # corresponding to the origin or supremum of every C2 group.
+        # It has the structure of `dim` times embedded lists each containing
+        # these vertices as the entire complex grows. Bounds[0] has to be done
+        # outside the loops before we have symmetric containers.
+        # NOTE: This means that bounds[0][1] must always exist
+        C0x = [[self.V[vot]]]
+        a_vo = copy.copy(list(origin))
+        a_vo[0] = vut[0]  # Update aN Origin
+        a_vo = self.V[tuple(a_vo)]
+        # self.V[vot].connect(self.V[tuple(a_vo)])
+        self.V[vot].connect(a_vo)
+        yield a_vo.x
+        C1x = [[a_vo]]
+        # C1x = [[self.V[tuple(a_vo)]]]
+        ab_C = []  # Container for a + b operations
+
+        # Loop over remaining bounds
+        for i, x in enumerate(bounds[1:]):
+            # Update lower and upper containers
+            C0x.append([])
+            C1x.append([])
+            # try to access a second bound (if not, C1 is symmetric)
+            try:
+                # Early try so that we don't have to copy the cache before
+                # moving on to next C1/C2: Try to add the operation of a new
+                # C2 product by accessing the upper bound
+                x[1]
+                # Copy lists for iteration
+                cC0x = [x[:] for x in C0x[:i + 1]]
+                cC1x = [x[:] for x in C1x[:i + 1]]
+                for j, (VL, VU) in enumerate(zip(cC0x, cC1x)):
+                    for k, (vl, vu) in enumerate(zip(VL, VU)):
+                        # Build aN vertices for each lower-upper pair in N:
+                        a_vl = list(vl.x)
+                        a_vu = list(vu.x)
+                        a_vl[i + 1] = vut[i + 1]
+                        a_vu[i + 1] = vut[i + 1]
+                        a_vl = self.V[tuple(a_vl)]
+
+                        # Connect vertices in N to corresponding vertices
+                        # in aN:
+                        vl.connect(a_vl)
+
+                        yield a_vl.x
+
+                        a_vu = self.V[tuple(a_vu)]
+                        # Connect vertices in N to corresponding vertices
+                        # in aN:
+                        vu.connect(a_vu)
+
+                        # Connect new vertex pair in aN:
+                        a_vl.connect(a_vu)
+
+                        # Connect lower pair to upper (triangulation
+                        # operation of a + b (two arbitrary operations):
+                        vl.connect(a_vu)
+                        ab_C.append((vl, a_vu))
+
+                        # Update the containers
+                        C0x[i + 1].append(vl)
+                        C0x[i + 1].append(vu)
+                        C1x[i + 1].append(a_vl)
+                        C1x[i + 1].append(a_vu)
+
+                        # Update old containers
+                        C0x[j].append(a_vl)
+                        C1x[j].append(a_vu)
+
+                        # Yield new points
+                        yield a_vu.x
+
+                # Try to connect aN lower source of previous a + b
+                # operation with a aN vertex
+                ab_Cc = copy.copy(ab_C)
+
+                for vp in ab_Cc:
+                    b_v = list(vp[0].x)
+                    ab_v = list(vp[1].x)
+                    b_v[i + 1] = vut[i + 1]
+                    ab_v[i + 1] = vut[i + 1]
+                    b_v = self.V[tuple(b_v)]  # b + vl
+                    ab_v = self.V[tuple(ab_v)]  # b + a_vl
+                    # Note o---o is already connected
+                    vp[0].connect(ab_v)  # o-s
+                    b_v.connect(ab_v)  # s-s
+
+                    # Add new list of cross pairs
+                    ab_C.append((vp[0], ab_v))
+                    ab_C.append((b_v, ab_v))
+
+            except IndexError:
+                cC0x = C0x[i]
+                cC1x = C1x[i]
+                VL, VU = cC0x, cC1x
+                for k, (vl, vu) in enumerate(zip(VL, VU)):
+                    # Build aN vertices for each lower-upper pair in N:
+                    a_vu = list(vu.x)
+                    a_vu[i + 1] = vut[i + 1]
+                    # Connect vertices in N to corresponding vertices
+                    # in aN:
+                    a_vu = self.V[tuple(a_vu)]
+                    # Connect vertices in N to corresponding vertices
+                    # in aN:
+                    vu.connect(a_vu)
+                    # Connect new vertex pair in aN:
+                    # a_vl.connect(a_vu)
+                    # Connect lower pair to upper (triangulation
+                    # operation of a + b (two arbitrary operations):
+                    vl.connect(a_vu)
+                    ab_C.append((vl, a_vu))
+                    C0x[i + 1].append(vu)
+                    C1x[i + 1].append(a_vu)
+                    # Yield new points
+                    a_vu.connect(self.V[vut])
+                    yield a_vu.x
+                    ab_Cc = copy.copy(ab_C)
+                    for vp in ab_Cc:
+                        if vp[1].x[i] == vut[i]:
+                            ab_v = list(vp[1].x)
+                            ab_v[i + 1] = vut[i + 1]
+                            ab_v = self.V[tuple(ab_v)]  # b + a_vl
+                            # Note o---o is already connected
+                            vp[0].connect(ab_v)  # o-s
+
+                            # Add new list of cross pairs
+                            ab_C.append((vp[0], ab_v))
+
+        # Clean class trash
+        try:
+            del C0x
+            del cC0x
+            del C1x
+            del cC1x
+            del ab_C
+            del ab_Cc
+        except UnboundLocalError:
+            pass
+
+        # Extra yield to ensure that the triangulation is completed
+        if centroid:
+            vo = self.V[vot]
+            vs = self.V[vut]
+            # Disconnect the origin and supremum
+            vo.disconnect(vs)
+            # Build centroid
+            vc = self.split_edge(vot, vut)
+            for v in vo.nn:
+                v.connect(vc)
+            yield vc.x
+            return vc.x
+        else:
+            yield vut
+            return vut
+
+    def triangulate(self, n=None, symmetry=None, centroid=True,
+                    printout=False):
+        """
+        Triangulate the initial domain, if n is not None then a limited number
+        of points will be generated
+
+        Parameters
+        ----------
+        n : int, Number of points to be sampled.
+        symmetry :
+
+            Ex. Dictionary/hashtable
+            f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2
+
+            symmetry = symmetry[0]: 0,  # Variable 1
+                       symmetry[1]: 0,  # symmetric to variable 1
+                       symmetry[2]: 0,  # symmetric to variable 1
+                       symmetry[3]: 3,  # Variable 4
+                       symmetry[4]: 3,  # symmetric to variable 4
+                       symmetry[5]: 3,  # symmetric to variable 4
+                        }
+        centroid : bool, if True add a central point to the hypercube
+        printout : bool, if True print out results
+
+        NOTES:
+        ------
+        Rather than using the combinatorial algorithm to connect vertices we
+        make the following observation:
+
+        The bound pairs are similar a C2 cyclic group and the structure is
+        formed using the cartesian product:
+
+        H = C2 x C2 x C2 ... x C2 (dim times)
+
+        So construct any normal subgroup N and consider H/N first, we connect
+        all vertices within N (ex. N is C2 (the first dimension), then we move
+        to a left coset aN (an operation moving around the defined H/N group by
+        for example moving from the lower bound in C2 (dimension 2) to the
+        higher bound in C2. During this operation connection all the vertices.
+        Now repeat the N connections. Note that these elements can be connected
+        in parrallel.
+        """
+        # Inherit class arguments
+        if symmetry is None:
+            symmetry = self.symmetry
+        # Build origin and supremum vectors
+        origin = [i[0] for i in self.bounds]
+        self.origin = origin
+        supremum = [i[1] for i in self.bounds]
+
+        self.supremum = supremum
+
+        if symmetry is None:
+            cbounds = self.bounds
+        else:
+            cbounds = copy.copy(self.bounds)
+            for i, j in enumerate(symmetry):
+                if i is not j:
+                    # pop second entry on second symmetry vars
+                    cbounds[i] = [self.bounds[symmetry[i]][0]]
+                    # Sole (first) entry is the sup value and there is no
+                    # origin:
+                    cbounds[i] = [self.bounds[symmetry[i]][1]]
+                    if (self.bounds[symmetry[i]] is not
+                            self.bounds[symmetry[j]]):
+                        logging.warning(f"Variable {i} was specified as "
+                                        f"symmetetric to variable {j}, however"
+                                        f", the bounds {i} ="
+                                        f" {self.bounds[symmetry[i]]} and {j}"
+                                        f" ="
+                                        f" {self.bounds[symmetry[j]]} do not "
+                                        f"match, the mismatch was ignored in "
+                                        f"the initial triangulation.")
+                        cbounds[i] = self.bounds[symmetry[j]]
+
+        if n is None:
+            # Build generator
+            self.cp = self.cyclic_product(cbounds, origin, supremum, centroid)
+            for i in self.cp:
+                i
+
+            try:
+                self.triangulated_vectors.append((tuple(self.origin),
+                                                  tuple(self.supremum)))
+            except (AttributeError, KeyError):
+                self.triangulated_vectors = [(tuple(self.origin),
+                                              tuple(self.supremum))]
+
+        else:
+            # Check if generator already exists
+            try:
+                self.cp
+            except (AttributeError, KeyError):
+                self.cp = self.cyclic_product(cbounds, origin, supremum,
+                                              centroid)
+
+            try:
+                while len(self.V.cache) < n:
+                    next(self.cp)
+            except StopIteration:
+                try:
+                    self.triangulated_vectors.append((tuple(self.origin),
+                                                      tuple(self.supremum)))
+                except (AttributeError, KeyError):
+                    self.triangulated_vectors = [(tuple(self.origin),
+                                                  tuple(self.supremum))]
+
+        if printout:
+            # for v in self.C0():
+            #   v.print_out()
+            for v in self.V.cache:
+                self.V[v].print_out()
+
+        return
+
+    def refine(self, n=1):
+        if n is None:
+            try:
+                self.triangulated_vectors
+                self.refine_all()
+                return
+            except AttributeError as ae:
+                if str(ae) == "'Complex' object has no attribute " \
+                              "'triangulated_vectors'":
+                    self.triangulate(symmetry=self.symmetry)
+                    return
+                else:
+                    raise
+
+        nt = len(self.V.cache) + n  # Target number of total vertices
+        # In the outer while loop we iterate until we have added an extra `n`
+        # vertices to the complex:
+        while len(self.V.cache) < nt:  # while loop 1
+            try:  # try 1
+                # Try to access triangulated_vectors, this should only be
+                # defined if an initial triangulation has already been
+                # performed:
+                self.triangulated_vectors
+                # Try a usual iteration of the current generator, if it
+                # does not exist or is exhausted then produce a new generator
+                try:  # try 2
+                    next(self.rls)
+                except (AttributeError, StopIteration, KeyError):
+                    vp = self.triangulated_vectors[0]
+                    self.rls = self.refine_local_space(*vp, bounds=self.bounds)
+                    next(self.rls)
+
+            except (AttributeError, KeyError):
+                # If an initial triangulation has not been completed, then
+                # we start/continue the initial triangulation targeting `nt`
+                # vertices, if nt is greater than the initial number of
+                # vertices then the `refine` routine will move back to try 1.
+                self.triangulate(nt, self.symmetry)
+        return
+
+    def refine_all(self, centroids=True):
+        """Refine the entire domain of the current complex."""
+        try:
+            self.triangulated_vectors
+            tvs = copy.copy(self.triangulated_vectors)
+            for i, vp in enumerate(tvs):
+                self.rls = self.refine_local_space(*vp, bounds=self.bounds)
+                for i in self.rls:
+                    i
+        except AttributeError as ae:
+            if str(ae) == "'Complex' object has no attribute " \
+                          "'triangulated_vectors'":
+                self.triangulate(symmetry=self.symmetry, centroid=centroids)
+            else:
+                raise
+
+        # This adds a centroid to every new sub-domain generated and defined
+        # by self.triangulated_vectors, in addition the vertices ! to complete
+        # the triangulation
+        return
+
+    def refine_local_space(self, origin, supremum, bounds, centroid=1):
+        # Copy for later removal
+        origin_c = copy.copy(origin)
+        supremum_c = copy.copy(supremum)
+
+        # Initiate local variables redefined in later inner `for` loop:
+        vl, vu, a_vu = None, None, None
+
+        # Change the vector orientation so that it is only increasing
+        s_ov = list(origin)
+        s_origin = list(origin)
+        s_sv = list(supremum)
+        s_supremum = list(supremum)
+        for i, vi in enumerate(s_origin):
+            if s_ov[i] > s_sv[i]:
+                s_origin[i] = s_sv[i]
+                s_supremum[i] = s_ov[i]
+
+        vot = tuple(s_origin)
+        vut = tuple(s_supremum)  # Hyperrectangle supremum
+
+        vo = self.V[vot]  # initiate if doesn't exist yet
+        vs = self.V[vut]
+        # Start by finding the old centroid of the new space:
+        vco = self.split_edge(vo.x, vs.x)  # Split in case not centroid arg
+
+        # Find set of extreme vertices in current local space
+        sup_set = copy.copy(vco.nn)
+        # Cyclic group approach with second x_l --- x_u operation.
+
+        # These containers store the "lower" and "upper" vertices
+        # corresponding to the origin or supremum of every C2 group.
+        # It has the structure of `dim` times embedded lists each containing
+        # these vertices as the entire complex grows. Bounds[0] has to be done
+        # outside the loops before we have symmetric containers.
+        # NOTE: This means that bounds[0][1] must always exist
+
+        a_vl = copy.copy(list(vot))
+        a_vl[0] = vut[0]  # Update aN Origin
+        if tuple(a_vl) not in self.V.cache:
+            vo = self.V[vot]  # initiate if doesn't exist yet
+            vs = self.V[vut]
+            # Start by finding the old centroid of the new space:
+            vco = self.split_edge(vo.x, vs.x)  # Split in case not centroid arg
+
+            # Find set of extreme vertices in current local space
+            sup_set = copy.copy(vco.nn)
+            a_vl = copy.copy(list(vot))
+            a_vl[0] = vut[0]  # Update aN Origin
+            a_vl = self.V[tuple(a_vl)]
+        else:
+            a_vl = self.V[tuple(a_vl)]
+
+        c_v = self.split_edge(vo.x, a_vl.x)
+        c_v.connect(vco)
+        yield c_v.x
+        Cox = [[vo]]
+        Ccx = [[c_v]]
+        Cux = [[a_vl]]
+        ab_C = []  # Container for a + b operations
+        s_ab_C = []  # Container for symmetric a + b operations
+
+        # Loop over remaining bounds
+        for i, x in enumerate(bounds[1:]):
+            # Update lower and upper containers
+            Cox.append([])
+            Ccx.append([])
+            Cux.append([])
+            # try to access a second bound (if not, C1 is symmetric)
+            try:
+                t_a_vl = list(vot)
+                t_a_vl[i + 1] = vut[i + 1]
+
+                # New: lists are used anyway, so copy all
+                # %%
+                # Copy lists for iteration
+                cCox = [x[:] for x in Cox[:i + 1]]
+                cCcx = [x[:] for x in Ccx[:i + 1]]
+                cCux = [x[:] for x in Cux[:i + 1]]
+                # Try to connect aN lower source of previous a + b
+                # operation with a aN vertex
+                ab_Cc = copy.copy(ab_C)  # NOTE: We append ab_C in the
+                # (VL, VC, VU) for-loop, but we use the copy of the list in the
+                # ab_Cc for-loop.
+                s_ab_Cc = copy.copy(s_ab_C)
+
+                # Early try so that we don't have to copy the cache before
+                # moving on to next C1/C2: Try to add the operation of a new
+                # C2 product by accessing the upper bound
+                if tuple(t_a_vl) not in self.V.cache:
+                    # Raise error to continue symmetric refine
+                    raise IndexError
+                t_a_vu = list(vut)
+                t_a_vu[i + 1] = vut[i + 1]
+                if tuple(t_a_vu) not in self.V.cache:
+                    # Raise error to continue symmetric refine:
+                    raise IndexError
+
+                for vectors in s_ab_Cc:
+                    # s_ab_C.append([c_vc, vl, vu, a_vu])
+                    bc_vc = list(vectors[0].x)
+                    b_vl = list(vectors[1].x)
+                    b_vu = list(vectors[2].x)
+                    ba_vu = list(vectors[3].x)
+
+                    bc_vc[i + 1] = vut[i + 1]
+                    b_vl[i + 1] = vut[i + 1]
+                    b_vu[i + 1] = vut[i + 1]
+                    ba_vu[i + 1] = vut[i + 1]
+
+                    bc_vc = self.V[tuple(bc_vc)]
+                    bc_vc.connect(vco)  # NOTE: Unneeded?
+                    yield bc_vc
+
+                    # Split to centre, call this centre group "d = 0.5*a"
+                    d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
+                    d_bc_vc.connect(bc_vc)
+                    d_bc_vc.connect(vectors[1])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[2])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[3])  # Connect all to centroid
+                    yield d_bc_vc.x
+                    b_vl = self.V[tuple(b_vl)]
+                    bc_vc.connect(b_vl)  # Connect aN cross pairs
+                    d_bc_vc.connect(b_vl)  # Connect all to centroid
+
+                    yield b_vl
+                    b_vu = self.V[tuple(b_vu)]
+                    bc_vc.connect(b_vu)  # Connect aN cross pairs
+                    d_bc_vc.connect(b_vu)  # Connect all to centroid
+
+                    b_vl_c = self.split_edge(b_vu.x, b_vl.x)
+                    bc_vc.connect(b_vl_c)
+
+                    yield b_vu
+                    ba_vu = self.V[tuple(ba_vu)]
+                    bc_vc.connect(ba_vu)  # Connect aN cross pairs
+                    d_bc_vc.connect(ba_vu)  # Connect all to centroid
+
+                    # Split the a + b edge of the initial triangulation:
+                    os_v = self.split_edge(vectors[1].x, ba_vu.x)  # o-s
+                    ss_v = self.split_edge(b_vl.x, ba_vu.x)  # s-s
+                    b_vu_c = self.split_edge(b_vu.x, ba_vu.x)
+                    bc_vc.connect(b_vu_c)
+                    yield os_v.x  # often equal to vco, but not always
+                    yield ss_v.x  # often equal to bc_vu, but not always
+                    yield ba_vu
+                    # Split remaining to centre, call this centre group
+                    # "d = 0.5*a"
+                    d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    yield d_bc_vc.x
+                    d_b_vl = self.split_edge(vectors[1].x, b_vl.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_b_vl)  # Connect dN cross pairs
+                    yield d_b_vl.x
+                    d_b_vu = self.split_edge(vectors[2].x, b_vu.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_b_vu)  # Connect dN cross pairs
+                    yield d_b_vu.x
+                    d_ba_vu = self.split_edge(vectors[3].x, ba_vu.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_ba_vu)  # Connect dN cross pairs
+                    yield d_ba_vu
+
+                    # comb = [c_vc, vl, vu, a_vl, a_vu,
+                    #       bc_vc, b_vl, b_vu, ba_vl, ba_vu]
+                    comb = [vl, vu, a_vu,
+                            b_vl, b_vu, ba_vu]
+                    comb_iter = itertools.combinations(comb, 2)
+                    for vecs in comb_iter:
+                        self.split_edge(vecs[0].x, vecs[1].x)
+                    # Add new list of cross pairs
+                    ab_C.append((d_bc_vc, vectors[1], b_vl, a_vu, ba_vu))
+                    ab_C.append((d_bc_vc, vl, b_vl, a_vu, ba_vu))  # = prev
+
+                for vectors in ab_Cc:
+                    bc_vc = list(vectors[0].x)
+                    b_vl = list(vectors[1].x)
+                    b_vu = list(vectors[2].x)
+                    ba_vl = list(vectors[3].x)
+                    ba_vu = list(vectors[4].x)
+                    bc_vc[i + 1] = vut[i + 1]
+                    b_vl[i + 1] = vut[i + 1]
+                    b_vu[i + 1] = vut[i + 1]
+                    ba_vl[i + 1] = vut[i + 1]
+                    ba_vu[i + 1] = vut[i + 1]
+                    bc_vc = self.V[tuple(bc_vc)]
+                    bc_vc.connect(vco)  # NOTE: Unneeded?
+                    yield bc_vc
+
+                    # Split to centre, call this centre group "d = 0.5*a"
+                    d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
+                    d_bc_vc.connect(bc_vc)
+                    d_bc_vc.connect(vectors[1])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[2])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[3])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[4])  # Connect all to centroid
+                    yield d_bc_vc.x
+                    b_vl = self.V[tuple(b_vl)]
+                    bc_vc.connect(b_vl)  # Connect aN cross pairs
+                    d_bc_vc.connect(b_vl)  # Connect all to centroid
+                    yield b_vl
+                    b_vu = self.V[tuple(b_vu)]
+                    bc_vc.connect(b_vu)  # Connect aN cross pairs
+                    d_bc_vc.connect(b_vu)  # Connect all to centroid
+                    yield b_vu
+                    ba_vl = self.V[tuple(ba_vl)]
+                    bc_vc.connect(ba_vl)  # Connect aN cross pairs
+                    d_bc_vc.connect(ba_vl)  # Connect all to centroid
+                    self.split_edge(b_vu.x, ba_vl.x)
+                    yield ba_vl
+                    ba_vu = self.V[tuple(ba_vu)]
+                    bc_vc.connect(ba_vu)  # Connect aN cross pairs
+                    d_bc_vc.connect(ba_vu)  # Connect all to centroid
+                    # Split the a + b edge of the initial triangulation:
+                    os_v = self.split_edge(vectors[1].x, ba_vu.x)  # o-s
+                    ss_v = self.split_edge(b_vl.x, ba_vu.x)  # s-s
+                    yield os_v.x  # often equal to vco, but not always
+                    yield ss_v.x  # often equal to bc_vu, but not always
+                    yield ba_vu
+                    # Split remaining to centre, call this centre group
+                    # "d = 0.5*a"
+                    d_bc_vc = self.split_edge(vectors[0].x, bc_vc.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    yield d_bc_vc.x
+                    d_b_vl = self.split_edge(vectors[1].x, b_vl.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_b_vl)  # Connect dN cross pairs
+                    yield d_b_vl.x
+                    d_b_vu = self.split_edge(vectors[2].x, b_vu.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_b_vu)  # Connect dN cross pairs
+                    yield d_b_vu.x
+                    d_ba_vl = self.split_edge(vectors[3].x, ba_vl.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_ba_vl)  # Connect dN cross pairs
+                    yield d_ba_vl
+                    d_ba_vu = self.split_edge(vectors[4].x, ba_vu.x)
+                    d_bc_vc.connect(vco)  # NOTE: Unneeded?
+                    d_bc_vc.connect(d_ba_vu)  # Connect dN cross pairs
+                    yield d_ba_vu
+                    c_vc, vl, vu, a_vl, a_vu = vectors
+
+                    comb = [vl, vu, a_vl, a_vu,
+                            b_vl, b_vu, ba_vl, ba_vu]
+                    comb_iter = itertools.combinations(comb, 2)
+                    for vecs in comb_iter:
+                        self.split_edge(vecs[0].x, vecs[1].x)
+
+                    # Add new list of cross pairs
+                    ab_C.append((bc_vc, b_vl, b_vu, ba_vl, ba_vu))
+                    ab_C.append((d_bc_vc, d_b_vl, d_b_vu, d_ba_vl, d_ba_vu))
+                    ab_C.append((d_bc_vc, vectors[1], b_vl, a_vu, ba_vu))
+                    ab_C.append((d_bc_vc, vu, b_vu, a_vl, ba_vl))
+
+                for j, (VL, VC, VU) in enumerate(zip(cCox, cCcx, cCux)):
+                    for k, (vl, vc, vu) in enumerate(zip(VL, VC, VU)):
+                        # Build aN vertices for each lower-upper C3 group in N:
+                        a_vl = list(vl.x)
+                        a_vu = list(vu.x)
+                        a_vl[i + 1] = vut[i + 1]
+                        a_vu[i + 1] = vut[i + 1]
+                        a_vl = self.V[tuple(a_vl)]
+                        a_vu = self.V[tuple(a_vu)]
+                        # Note, build (a + vc) later for consistent yields
+                        # Split the a + b edge of the initial triangulation:
+                        c_vc = self.split_edge(vl.x, a_vu.x)
+                        self.split_edge(vl.x, vu.x)  # Equal to vc
+                        # Build cN vertices for each lower-upper C3 group in N:
+                        c_vc.connect(vco)
+                        c_vc.connect(vc)
+                        c_vc.connect(vl)  # Connect c + ac operations
+                        c_vc.connect(vu)  # Connect c + ac operations
+                        c_vc.connect(a_vl)  # Connect c + ac operations
+                        c_vc.connect(a_vu)  # Connect c + ac operations
+                        yield c_vc.x
+                        c_vl = self.split_edge(vl.x, a_vl.x)
+                        c_vl.connect(vco)
+                        c_vc.connect(c_vl)  # Connect cN group vertices
+                        yield c_vl.x
+                        # yield at end of loop:
+                        c_vu = self.split_edge(vu.x, a_vu.x)
+                        c_vu.connect(vco)
+                        # Connect remaining cN group vertices
+                        c_vc.connect(c_vu)  # Connect cN group vertices
+                        yield c_vu.x
+
+                        a_vc = self.split_edge(a_vl.x, a_vu.x)  # is (a + vc) ?
+                        a_vc.connect(vco)
+                        a_vc.connect(c_vc)
+
+                        # Storage for connecting c + ac operations:
+                        ab_C.append((c_vc, vl, vu, a_vl, a_vu))
+
+                        # Update the containers
+                        Cox[i + 1].append(vl)
+                        Cox[i + 1].append(vc)
+                        Cox[i + 1].append(vu)
+                        Ccx[i + 1].append(c_vl)
+                        Ccx[i + 1].append(c_vc)
+                        Ccx[i + 1].append(c_vu)
+                        Cux[i + 1].append(a_vl)
+                        Cux[i + 1].append(a_vc)
+                        Cux[i + 1].append(a_vu)
+
+                        # Update old containers
+                        Cox[j].append(c_vl)  # !
+                        Cox[j].append(a_vl)
+                        Ccx[j].append(c_vc)  # !
+                        Ccx[j].append(a_vc)  # !
+                        Cux[j].append(c_vu)  # !
+                        Cux[j].append(a_vu)
+
+                        # Yield new points
+                        yield a_vc.x
+
+            except IndexError:
+                for vectors in ab_Cc:
+                    ba_vl = list(vectors[3].x)
+                    ba_vu = list(vectors[4].x)
+                    ba_vl[i + 1] = vut[i + 1]
+                    ba_vu[i + 1] = vut[i + 1]
+                    ba_vu = self.V[tuple(ba_vu)]
+                    yield ba_vu
+                    d_bc_vc = self.split_edge(vectors[1].x, ba_vu.x)  # o-s
+                    yield ba_vu
+                    d_bc_vc.connect(vectors[1])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[2])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[3])  # Connect all to centroid
+                    d_bc_vc.connect(vectors[4])  # Connect all to centroid
+                    yield d_bc_vc.x
+                    ba_vl = self.V[tuple(ba_vl)]
+                    yield ba_vl
+                    d_ba_vl = self.split_edge(vectors[3].x, ba_vl.x)
+                    d_ba_vu = self.split_edge(vectors[4].x, ba_vu.x)
+                    d_ba_vc = self.split_edge(d_ba_vl.x, d_ba_vu.x)
+                    yield d_ba_vl
+                    yield d_ba_vu
+                    yield d_ba_vc
+                    c_vc, vl, vu, a_vl, a_vu = vectors
+                    comb = [vl, vu, a_vl, a_vu,
+                            ba_vl,
+                            ba_vu]
+                    comb_iter = itertools.combinations(comb, 2)
+                    for vecs in comb_iter:
+                        self.split_edge(vecs[0].x, vecs[1].x)
+
+                # Copy lists for iteration
+                cCox = Cox[i]
+                cCcx = Ccx[i]
+                cCux = Cux[i]
+                VL, VC, VU = cCox, cCcx, cCux
+                for k, (vl, vc, vu) in enumerate(zip(VL, VC, VU)):
+                    # Build aN vertices for each lower-upper pair in N:
+                    a_vu = list(vu.x)
+                    a_vu[i + 1] = vut[i + 1]
+
+                    # Connect vertices in N to corresponding vertices
+                    # in aN:
+                    a_vu = self.V[tuple(a_vu)]
+                    yield a_vl.x
+                    # Split the a + b edge of the initial triangulation:
+                    c_vc = self.split_edge(vl.x, a_vu.x)
+                    self.split_edge(vl.x, vu.x)  # Equal to vc
+                    c_vc.connect(vco)
+                    c_vc.connect(vc)
+                    c_vc.connect(vl)  # Connect c + ac operations
+                    c_vc.connect(vu)  # Connect c + ac operations
+                    c_vc.connect(a_vu)  # Connect c + ac operations
+                    yield (c_vc.x)
+                    c_vu = self.split_edge(vu.x,
+                                           a_vu.x)  # yield at end of loop
+                    c_vu.connect(vco)
+                    # Connect remaining cN group vertices
+                    c_vc.connect(c_vu)  # Connect cN group vertices
+                    yield (c_vu.x)
+
+                    # Update the containers
+                    Cox[i + 1].append(vu)
+                    Ccx[i + 1].append(c_vu)
+                    Cux[i + 1].append(a_vu)
+
+                    # Update old containers
+                    s_ab_C.append([c_vc, vl, vu, a_vu])
+
+                    yield a_vu.x
+
+        # Clean class trash
+        try:
+            del Cox
+            del Ccx
+            del Cux
+            del ab_C
+            del ab_Cc
+        except UnboundLocalError:
+            pass
+
+        try:
+            self.triangulated_vectors.remove((tuple(origin_c),
+                                              tuple(supremum_c)))
+        except ValueError:
+            # Turn this into a logging warning?
+            pass
+        # Add newly triangulated vectors:
+        for vs in sup_set:
+            self.triangulated_vectors.append((tuple(vco.x), tuple(vs.x)))
+
+        # Extra yield to ensure that the triangulation is completed
+        if centroid:
+            vcn_set = set()
+            c_nn_lists = []
+            for vs in sup_set:
+                # Build centroid
+                c_nn = self.vpool(vco.x, vs.x)
+                try:
+                    c_nn.remove(vcn_set)
+                except KeyError:
+                    pass
+                c_nn_lists.append(c_nn)
+
+            for c_nn in c_nn_lists:
+                try:
+                    c_nn.remove(vcn_set)
+                except KeyError:
+                    pass
+
+            for vs, c_nn in zip(sup_set, c_nn_lists):
+                # Build centroid
+                vcn = self.split_edge(vco.x, vs.x)
+                vcn_set.add(vcn)
+                try:  # Shouldn't be needed?
+                    c_nn.remove(vcn_set)
+                except KeyError:
+                    pass
+                for vnn in c_nn:
+                    vcn.connect(vnn)
+                yield vcn.x
+        else:
+            pass
+
+        yield vut
+        return
+
+    def refine_star(self, v):
+        """Refine the star domain of a vertex `v`."""
+        # Copy lists before iteration
+        vnn = copy.copy(v.nn)
+        v1nn = []
+        d_v0v1_set = set()
+        for v1 in vnn:
+            v1nn.append(copy.copy(v1.nn))
+
+        for v1, v1nn in zip(vnn, v1nn):
+            vnnu = v1nn.intersection(vnn)
+
+            d_v0v1 = self.split_edge(v.x, v1.x)
+            for o_d_v0v1 in d_v0v1_set:
+                d_v0v1.connect(o_d_v0v1)
+            d_v0v1_set.add(d_v0v1)
+            for v2 in vnnu:
+                d_v1v2 = self.split_edge(v1.x, v2.x)
+                d_v0v1.connect(d_v1v2)
+        return
+
+    @lru_cache(maxsize=None)
+    def split_edge(self, v1, v2):
+        v1 = self.V[v1]
+        v2 = self.V[v2]
+        # Destroy original edge, if it exists:
+        v1.disconnect(v2)
+        # Compute vertex on centre of edge:
+        try:
+            vct = (v2.x_a - v1.x_a) / 2.0 + v1.x_a
+        except TypeError:  # Allow for decimal operations
+            vct = (v2.x_a - v1.x_a) / decimal.Decimal(2.0) + v1.x_a
+
+        vc = self.V[tuple(vct)]
+        # Connect to original 2 vertices to the new centre vertex
+        vc.connect(v1)
+        vc.connect(v2)
+        return vc
+
+    def vpool(self, origin, supremum):
+        vot = tuple(origin)
+        vst = tuple(supremum)
+        # Initiate vertices in case they don't exist
+        vo = self.V[vot]
+        vs = self.V[vst]
+
+        # Remove origin - supremum disconnect
+
+        # Find the lower/upper bounds of the refinement hyperrectangle
+        bl = list(vot)
+        bu = list(vst)
+        for i, (voi, vsi) in enumerate(zip(vot, vst)):
+            if bl[i] > vsi:
+                bl[i] = vsi
+            if bu[i] < voi:
+                bu[i] = voi
+
+        #      NOTE: This is mostly done with sets/lists because we aren't sure
+        #            how well the numpy arrays will scale to thousands of
+        #             dimensions.
+        vn_pool = set()
+        vn_pool.update(vo.nn)
+        vn_pool.update(vs.nn)
+        cvn_pool = copy.copy(vn_pool)
+        for vn in cvn_pool:
+            for i, xi in enumerate(vn.x):
+                if bl[i] <= xi <= bu[i]:
+                    pass
+                else:
+                    try:
+                        vn_pool.remove(vn)
+                    except KeyError:
+                        pass  # NOTE: Not all neigbouds are in initial pool
+        return vn_pool
+
+    def vf_to_vv(self, vertices, simplices):
+        """
+        Convert a vertex-face mesh to a vertex-vertex mesh used by this class
+
+        Parameters
+        ----------
+        vertices : list
+            Vertices
+        simplices : list
+            Simplices
+        """
+        if self.dim > 1:
+            for s in simplices:
+                edges = itertools.combinations(s, self.dim)
+                for e in edges:
+                    self.V[tuple(vertices[e[0]])].connect(
+                        self.V[tuple(vertices[e[1]])])
+        else:
+            for e in simplices:
+                self.V[tuple(vertices[e[0]])].connect(
+                    self.V[tuple(vertices[e[1]])])
+        return
+
+    def connect_vertex_non_symm(self, v_x, near=None):
+        """
+        Adds a vertex at coords v_x to the complex that is not symmetric to the
+        initial triangulation and sub-triangulation.
+
+        If near is specified (for example; a star domain or collections of
+        cells known to contain v) then only those simplices containd in near
+        will be searched, this greatly speeds up the process.
+
+        If near is not specified this method will search the entire simplicial
+        complex structure.
+
+        Parameters
+        ----------
+        v_x : tuple
+            Coordinates of non-symmetric vertex
+        near : set or list
+            List of vertices, these are points near v to check for
+        """
+        if near is None:
+            star = self.V
+        else:
+            star = near
+        # Create the vertex origin
+        if tuple(v_x) in self.V.cache:
+            if self.V[v_x] in self.V_non_symm:
+                pass
+            else:
+                return
+
+        self.V[v_x]
+        found_nn = False
+        S_rows = []
+        for v in star:
+            S_rows.append(v.x)
+
+        S_rows = numpy.array(S_rows)
+        A = numpy.array(S_rows) - numpy.array(v_x)
+        # Iterate through all the possible simplices of S_rows
+        for s_i in itertools.combinations(range(S_rows.shape[0]),
+                                          r=self.dim + 1):
+            # Check if connected, else s_i is not a simplex
+            valid_simplex = True
+            for i in itertools.combinations(s_i, r=2):
+                # Every combination of vertices must be connected, we check of
+                # the current iteration of all combinations of s_i are
+                # connected we break the loop if it is not.
+                if ((self.V[tuple(S_rows[i[1]])] not in
+                        self.V[tuple(S_rows[i[0]])].nn)
+                    and (self.V[tuple(S_rows[i[0]])] not in
+                         self.V[tuple(S_rows[i[1]])].nn)):
+                    valid_simplex = False
+                    break
+
+            S = S_rows[tuple([s_i])]
+            if valid_simplex:
+                if self.deg_simplex(S, proj=None):
+                    valid_simplex = False
+
+            # If s_i is a valid simplex we can test if v_x is inside si
+            if valid_simplex:
+                # Find the A_j0 value from the precalculated values
+                A_j0 = A[tuple([s_i])]
+                if self.in_simplex(S, v_x, A_j0):
+                    found_nn = True
+                    # breaks the main for loop, s_i is the target simplex:
+                    break
+
+        # Connect the simplex to point
+        if found_nn:
+            for i in s_i:
+                self.V[v_x].connect(self.V[tuple(S_rows[i])])
+        # Attached the simplex to storage for all non-symmetric vertices
+        self.V_non_symm.append(self.V[v_x])
+        # this bool value indicates a successful connection if True:
+        return found_nn
+
+    def in_simplex(self, S, v_x, A_j0=None):
+        """Check if a vector v_x is in simplex `S`.
+
+        Parameters
+        ----------
+        S : array_like
+            Array containing simplex entries of vertices as rows
+        v_x :
+            A candidate vertex
+        A_j0 : array, optional,
+            Allows for A_j0 to be pre-calculated
+
+        Returns
+        -------
+        res : boolean
+            True if `v_x` is in `S`
+        """
+        A_11 = numpy.delete(S, 0, 0) - S[0]
+
+        sign_det_A_11 = numpy.sign(numpy.linalg.det(A_11))
+        if sign_det_A_11 == 0:
+            # NOTE: We keep the variable A_11, but we loop through A_jj
+            # ind=
+            # while sign_det_A_11 == 0:
+            #    A_11 = numpy.delete(S, ind, 0) - S[ind]
+            #    sign_det_A_11 = numpy.sign(numpy.linalg.det(A_11))
+
+            sign_det_A_11 = -1  # TODO: Choose another det of j instead?
+            # TODO: Unlikely to work in many cases
+
+        if A_j0 is None:
+            A_j0 = S - v_x
+
+        for d in range(self.dim + 1):
+            det_A_jj = (-1)**d * sign_det_A_11
+            # TODO: Note that scipy might be faster to add as an optional
+            #       dependency
+            sign_det_A_j0 = numpy.sign(numpy.linalg.det(numpy.delete(A_j0, d,
+                                                                     0)))
+            # TODO: Note if sign_det_A_j0 == then the point is coplanar to the
+            #       current simplex facet, so perhaps return True and attach?
+            if det_A_jj == sign_det_A_j0:
+                continue
+            else:
+                return False
+
+        return True
+
+    def deg_simplex(self, S, proj=None):
+        """Test a simplex S for degeneracy (linear dependence in R^dim).
+
+        Parameters
+        ----------
+        S : np.array
+            Simplex with rows as vertex vectors
+        proj : array, optional,
+            If the projection S[1:] - S[0] is already
+            computed it can be added as an optional argument.
+        """
+        # Strategy: we test all combination of faces, if any of the
+        # determinants are zero then the vectors lie on the same face and is
+        # therefore linearly dependent in the space of R^dim
+        if proj is None:
+            proj = S[1:] - S[0]
+
+        # TODO: Is checking the projection of one vertex against faces of other
+        #       vertices sufficient? Or do we need to check more vertices in
+        #       dimensions higher than 2?
+        # TODO: Literature seems to suggest using proj.T, but why is this
+        #       needed?
+        if numpy.linalg.det(proj) == 0.0:  # TODO: Repalace with tolerance?
+            return True  # Simplex is degenerate
+        else:
+            return False  # Simplex is not degenerate
diff --git a/scipy/optimize/_shgo_lib/_vertex.py b/scipy/optimize/_shgo_lib/_vertex.py
new file mode 100644
index 000000000000..c983618b6d16
--- /dev/null
+++ b/scipy/optimize/_shgo_lib/_vertex.py
@@ -0,0 +1,458 @@
+import collections
+from abc import ABC, abstractmethod
+
+import numpy as np
+
+from scipy._lib._util import MapWrapper
+
+
+class VertexBase(ABC):
+    """
+    Base class for a vertex.
+    """
+    def __init__(self, x, nn=None, index=None):
+        """
+        Initiation of a vertex object.
+
+        Parameters
+        ----------
+        x : tuple or vector
+            The geometric location (domain).
+        nn : list, optional
+            Nearest neighbour list.
+        index : int, optional
+            Index of vertex.
+        """
+        self.x = x
+        self.hash = hash(self.x)  # Save precomputed hash
+
+        if nn is not None:
+            self.nn = set(nn)  # can use .indexupdate to add a new list
+        else:
+            self.nn = set()
+
+        self.index = index
+
+    def __hash__(self):
+        return self.hash
+
+    def __getattr__(self, item):
+        if item not in ['x_a']:
+            raise AttributeError(f"{type(self)} object has no attribute "
+                                 f"'{item}'")
+        if item == 'x_a':
+            self.x_a = np.array(self.x)
+            return self.x_a
+
+    @abstractmethod
+    def connect(self, v):
+        raise NotImplementedError("This method is only implemented with an "
+                                  "associated child of the base class.")
+
+    @abstractmethod
+    def disconnect(self, v):
+        raise NotImplementedError("This method is only implemented with an "
+                                  "associated child of the base class.")
+
+    def star(self):
+        """Returns the star domain ``st(v)`` of the vertex.
+
+        Parameters
+        ----------
+        v :
+            The vertex ``v`` in ``st(v)``
+
+        Returns
+        -------
+        st : set
+            A set containing all the vertices in ``st(v)``
+        """
+        self.st = self.nn
+        self.st.add(self)
+        return self.st
+
+
+class VertexScalarField(VertexBase):
+    """
+    Add homology properties of a scalar field f: R^n --> R associated with
+    the geometry built from the VertexBase class
+    """
+
+    def __init__(self, x, field=None, nn=None, index=None, field_args=(),
+                 g_cons=None, g_cons_args=()):
+        """
+        Parameters
+        ----------
+        x : tuple,
+            vector of vertex coordinates
+        field : callable, optional
+            a scalar field f: R^n --> R associated with the geometry
+        nn : list, optional
+            list of nearest neighbours
+        index : int, optional
+            index of the vertex
+        field_args : tuple, optional
+            additional arguments to be passed to field
+        g_cons : callable, optional
+            constraints on the vertex
+        g_cons_args : tuple, optional
+            additional arguments to be passed to g_cons
+
+        """
+        super().__init__(x, nn=nn, index=index)
+
+        # Note Vertex is only initiated once for all x so only
+        # evaluated once
+        # self.feasible = None
+
+        # self.f is externally defined by the cache to allow parallel
+        # processing
+        # None type that will break arithmetic operations unless defined
+        # self.f = None
+
+        self.check_min = True
+        self.check_max = True
+
+    def connect(self, v):
+        """Connects self to another vertex object v.
+
+        Parameters
+        ----------
+        v : VertexBase or VertexScalarField object
+        """
+        if v is not self and v not in self.nn:
+            self.nn.add(v)
+            v.nn.add(self)
+
+            # Flags for checking homology properties:
+            self.check_min = True
+            self.check_max = True
+            v.check_min = True
+            v.check_max = True
+
+    def disconnect(self, v):
+        if v in self.nn:
+            self.nn.remove(v)
+            v.nn.remove(self)
+
+            # Flags for checking homology properties:
+            self.check_min = True
+            self.check_max = True
+            v.check_min = True
+            v.check_max = True
+
+    def minimiser(self):
+        """Check whether this vertex is strictly less than all its
+           neighbours"""
+        if self.check_min:
+            self._min = all(self.f < v.f for v in self.nn)
+            self.check_min = False
+
+        return self._min
+
+    def maximiser(self):
+        """
+        Check whether this vertex is strictly greater than all its
+        neighbours.
+        """
+        if self.check_max:
+            self._max = all(self.f > v.f for v in self.nn)
+            self.check_max = False
+
+        return self._max
+
+
+class VertexVectorField(VertexBase):
+    """
+    Add homology properties of a scalar field f: R^n --> R^m associated with
+    the geometry built from the VertexBase class.
+    """
+
+    def __init__(self, x, sfield=None, vfield=None, field_args=(),
+                 vfield_args=(), g_cons=None,
+                 g_cons_args=(), nn=None, index=None):
+        super().__init__(x, nn=nn, index=index)
+
+        raise NotImplementedError("This class is still a work in progress")
+
+
+class VertexCacheBase:
+    """Base class for a vertex cache for a simplicial complex."""
+    def __init__(self):
+
+        self.cache = collections.OrderedDict()
+        self.nfev = 0  # Feasible points
+        self.index = -1
+
+    def __iter__(self):
+        for v in self.cache:
+            yield self.cache[v]
+        return
+
+    def size(self):
+        """Returns the size of the vertex cache."""
+        return self.index + 1
+
+    def print_out(self):
+        headlen = len(f"Vertex cache of size: {len(self.cache)}:")
+        print('=' * headlen)
+        print(f"Vertex cache of size: {len(self.cache)}:")
+        print('=' * headlen)
+        for v in self.cache:
+            self.cache[v].print_out()
+
+
+class VertexCube(VertexBase):
+    """Vertex class to be used for a pure simplicial complex with no associated
+    differential geometry (single level domain that exists in R^n)"""
+    def __init__(self, x, nn=None, index=None):
+        super().__init__(x, nn=nn, index=index)
+
+    def connect(self, v):
+        if v is not self and v not in self.nn:
+            self.nn.add(v)
+            v.nn.add(self)
+
+    def disconnect(self, v):
+        if v in self.nn:
+            self.nn.remove(v)
+            v.nn.remove(self)
+
+
+class VertexCacheIndex(VertexCacheBase):
+    def __init__(self):
+        """
+        Class for a vertex cache for a simplicial complex without an associated
+        field. Useful only for building and visualising a domain complex.
+
+        Parameters
+        ----------
+        """
+        super().__init__()
+        self.Vertex = VertexCube
+
+    def __getitem__(self, x, nn=None):
+        try:
+            return self.cache[x]
+        except KeyError:
+            self.index += 1
+            xval = self.Vertex(x, index=self.index)
+            # logging.info("New generated vertex at x = {}".format(x))
+            # NOTE: Surprisingly high performance increase if logging
+            # is commented out
+            self.cache[x] = xval
+            return self.cache[x]
+
+
+class VertexCacheField(VertexCacheBase):
+    def __init__(self, field=None, field_args=(), g_cons=None, g_cons_args=(),
+                 workers=1):
+        """
+        Class for a vertex cache for a simplicial complex with an associated
+        field.
+
+        Parameters
+        ----------
+        field : callable
+            Scalar or vector field callable.
+        field_args : tuple, optional
+            Any additional fixed parameters needed to completely specify the
+            field function
+        g_cons : dict or sequence of dict, optional
+            Constraints definition.
+            Function(s) ``R**n`` in the form::
+        g_cons_args : tuple, optional
+            Any additional fixed parameters needed to completely specify the
+            constraint functions
+        workers : int  optional
+            Uses `multiprocessing.Pool <multiprocessing>`) to compute the field
+             functions in parrallel.
+
+        """
+        super().__init__()
+        self.index = -1
+        self.Vertex = VertexScalarField
+        self.field = field
+        self.field_args = field_args
+        self.wfield = FieldWraper(field, field_args)  # if workers is not 1
+
+        self.g_cons = g_cons
+        self.g_cons_args = g_cons_args
+        self.wgcons = ConstraintWraper(g_cons, g_cons_args)
+        self.gpool = set()  # A set of tuples to process for feasibility
+
+        # Field processing objects
+        self.fpool = set()  # A set of tuples to process for scalar function
+        self.sfc_lock = False  # True if self.fpool is non-Empty
+
+        if workers == 1:
+            self.process_gpool = self.proc_gpool
+            if g_cons is None:
+                self.process_fpool = self.proc_fpool_nog
+            else:
+                self.process_fpool = self.proc_fpool_g
+        else:
+            self.workers = workers
+            self._mapwrapper = MapWrapper(workers)
+            self.pool = self._mapwrapper.pool
+            self.process_gpool = self.pproc_gpool
+            if g_cons is None:
+                self.process_fpool = self.pproc_fpool_nog
+            else:
+                self.process_fpool = self.pproc_fpool_g
+
+    def __getitem__(self, x, nn=None):
+        try:
+            return self.cache[x]
+        except KeyError:
+            self.index += 1
+            xval = self.Vertex(x, field=self.field, nn=nn, index=self.index,
+                               field_args=self.field_args,
+                               g_cons=self.g_cons,
+                               g_cons_args=self.g_cons_args)
+
+            self.cache[x] = xval  # Define in cache
+            self.gpool.add(xval)  # Add to pool for processing feasibility
+            self.fpool.add(xval)  # Add to pool for processing field values
+            return self.cache[x]
+
+    def __getstate__(self):
+        self_dict = self.__dict__.copy()
+        del self_dict['pool']
+        return self_dict
+
+    def process_pools(self):
+        if self.g_cons is not None:
+            self.process_gpool()
+        self.process_fpool()
+        self.proc_minimisers()
+
+    def feasibility_check(self, v):
+        v.feasible = True
+        for g, args in zip(self.g_cons, self.g_cons_args):
+            if g(v.x_a, *args) < 0.0:
+                v.f = np.inf
+                v.feasible = False
+                break
+
+    def compute_sfield(self, v):
+        """Compute the scalar field values of a vertex object `v`.
+
+        Parameters
+        ----------
+        v : VertexBase or VertexScalarField object
+        """
+        try:
+            v.f = self.field(v.x_a, *self.field_args)
+            self.nfev += 1
+        except AttributeError:
+            v.f = np.inf
+            # logging.warning(f"Field function not found at x = {self.x_a}")
+        if np.isnan(v.f):
+            v.f = np.inf
+
+    def proc_gpool(self):
+        """Process all constraints."""
+        if self.g_cons is not None:
+            for v in self.gpool:
+                self.feasibility_check(v)
+        # Clean the pool
+        self.gpool = set()
+
+    def pproc_gpool(self):
+        """Process all constraints in parallel."""
+        gpool_l = []
+        for v in self.gpool:
+            gpool_l.append(v.x_a)
+
+        G = self.pool.map(self.wgcons.gcons, gpool_l)
+        for v, g in zip(self.gpool, G):
+            v.feasible = g  # set vertex object attribute v.feasible = g (bool)
+
+    def proc_fpool_g(self):
+        """Process all field functions with constraints supplied."""
+        for v in self.fpool:
+            if v.feasible:
+                self.compute_sfield(v)
+        # Clean the pool
+        self.fpool = set()
+
+    def proc_fpool_nog(self):
+        """Process all field functions with no constraints supplied."""
+        for v in self.fpool:
+            self.compute_sfield(v)
+        # Clean the pool
+        self.fpool = set()
+
+    def pproc_fpool_g(self):
+        """
+        Process all field functions with constraints supplied in parallel.
+        """
+        self.wfield.func
+        fpool_l = []
+        for v in self.fpool:
+            if v.feasible:
+                fpool_l.append(v.x_a)
+            else:
+                v.f = np.inf
+        F = self.pool.map(self.wfield.func, fpool_l)
+        for va, f in zip(fpool_l, F):
+            vt = tuple(va)
+            self[vt].f = f  # set vertex object attribute v.f = f
+            self.nfev += 1
+        # Clean the pool
+        self.fpool = set()
+
+    def pproc_fpool_nog(self):
+        """
+        Process all field functions with no constraints supplied in parallel.
+        """
+        self.wfield.func
+        fpool_l = []
+        for v in self.fpool:
+            fpool_l.append(v.x_a)
+        F = self.pool.map(self.wfield.func, fpool_l)
+        for va, f in zip(fpool_l, F):
+            vt = tuple(va)
+            self[vt].f = f  # set vertex object attribute v.f = f
+            self.nfev += 1
+        # Clean the pool
+        self.fpool = set()
+
+    def proc_minimisers(self):
+        """Check for minimisers."""
+        for v in self:
+            v.minimiser()
+            v.maximiser()
+
+
+class ConstraintWraper:
+    """Object to wrap constraints to pass to `multiprocessing.Pool`."""
+    def __init__(self, g_cons, g_cons_args):
+        self.g_cons = g_cons
+        self.g_cons_args = g_cons_args
+
+    def gcons(self, v_x_a):
+        vfeasible = True
+        for g, args in zip(self.g_cons, self.g_cons_args):
+            if g(v_x_a, *args) < 0.0:
+                vfeasible = False
+                break
+        return vfeasible
+
+
+class FieldWraper:
+    """Object to wrap field to pass to `multiprocessing.Pool`."""
+    def __init__(self, field, field_args):
+        self.field = field
+        self.field_args = field_args
+
+    def func(self, v_x_a):
+        try:
+            v_f = self.field(v_x_a, *self.field_args)
+        except Exception:
+            v_f = np.inf
+        if np.isnan(v_f):
+            v_f = np.inf
+
+        return v_f
diff --git a/scipy/optimize/_shgo_lib/meson.build b/scipy/optimize/_shgo_lib/meson.build
index 98ffbd55e095..65e3b9abe14e 100644
--- a/scipy/optimize/_shgo_lib/meson.build
+++ b/scipy/optimize/_shgo_lib/meson.build
@@ -1,6 +1,7 @@
 py3.install_sources([
     '__init__.py',
-    'triangulation.py'
+    '_complex.py',
+    '_vertex.py'
   ],
   subdir: 'scipy/optimize/_shgo_lib'
 )
diff --git a/scipy/optimize/_shgo_lib/triangulation.py b/scipy/optimize/_shgo_lib/triangulation.py
deleted file mode 100644
index ffcd9c15c15b..000000000000
--- a/scipy/optimize/_shgo_lib/triangulation.py
+++ /dev/null
@@ -1,661 +0,0 @@
-import numpy as np
-import copy
-
-
-class Complex:
-    def __init__(self, dim, func, func_args=(), symmetry=False, bounds=None,
-                 g_cons=None, g_args=()):
-        self.dim = dim
-        self.bounds = bounds
-        self.symmetry = symmetry  # TODO: Define the functions to be used
-        #      here in init to avoid if checks
-        self.gen = 0
-        self.perm_cycle = 0
-
-        # Every cell is stored in a list of its generation,
-        # e.g., the initial cell is stored in self.H[0]
-        # 1st get new cells are stored in self.H[1] etc.
-        # When a cell is subgenerated it is removed from this list
-
-        self.H = []  # Storage structure of cells
-        # Cache of all vertices
-        self.V = VertexCache(func, func_args, bounds, g_cons, g_args)
-
-        # Generate n-cube here:
-        self.n_cube(dim, symmetry=symmetry)
-
-        # TODO: Assign functions to a the complex instead
-        if symmetry:
-            self.generation_cycle = 1
-            # self.centroid = self.C0()[-1].x
-            # self.C0.centroid = self.centroid
-        else:
-            self.add_centroid()
-
-        self.H.append([])
-        self.H[0].append(self.C0)
-        self.hgr = self.C0.homology_group_rank()
-        self.hgrd = 0  # Complex group rank differential
-        # self.hgr = self.C0.hg_n
-
-        # Build initial graph
-        self.graph_map()
-
-        self.performance = []
-        self.performance.append(0)
-        self.performance.append(0)
-
-    def __call__(self):
-        return self.H
-
-    def n_cube(self, dim, symmetry=False, printout=False):
-        """
-        Generate the simplicial triangulation of the N-D hypercube
-        containing 2**n vertices
-        """
-        origin = list(np.zeros(dim, dtype=int))
-        self.origin = origin
-        supremum = list(np.ones(dim, dtype=int))
-        self.supremum = supremum
-
-        # tuple versions for indexing
-        origintuple = tuple(origin)
-        supremumtuple = tuple(supremum)
-
-        x_parents = [origintuple]
-
-        if symmetry:
-            self.C0 = Simplex(0, 0, 0, self.dim)  # Initial cell object
-            self.C0.add_vertex(self.V[origintuple])
-
-            i_s = 0
-            self.perm_symmetry(i_s, x_parents, origin)
-            self.C0.add_vertex(self.V[supremumtuple])
-        else:
-            self.C0 = Cell(0, 0, origin, supremum)  # Initial cell object
-            self.C0.add_vertex(self.V[origintuple])
-            self.C0.add_vertex(self.V[supremumtuple])
-
-            i_parents = []
-            self.perm(i_parents, x_parents, origin)
-
-        if printout:
-            print("Initial hyper cube:")
-            for v in self.C0():
-                v.print_out()
-
-    def perm(self, i_parents, x_parents, xi):
-        # TODO: Cut out of for if outside linear constraint cutting planes
-        xi_t = tuple(xi)
-
-        # Construct required iterator
-        iter_range = [x for x in range(self.dim) if x not in i_parents]
-
-        for i in iter_range:
-            i2_parents = copy.copy(i_parents)
-            i2_parents.append(i)
-            xi2 = copy.copy(xi)
-            xi2[i] = 1
-            # Make new vertex list a hashable tuple
-            xi2_t = tuple(xi2)
-            # Append to cell
-            self.C0.add_vertex(self.V[xi2_t])
-            # Connect neighbors and vice versa
-            # Parent point
-            self.V[xi2_t].connect(self.V[xi_t])
-
-            # Connect all family of simplices in parent containers
-            for x_ip in x_parents:
-                self.V[xi2_t].connect(self.V[x_ip])
-
-            x_parents2 = copy.copy(x_parents)
-            x_parents2.append(xi_t)
-
-            # Permutate
-            self.perm(i2_parents, x_parents2, xi2)
-
-    def perm_symmetry(self, i_s, x_parents, xi):
-        # TODO: Cut out of for if outside linear constraint cutting planes
-        xi_t = tuple(xi)
-        xi2 = copy.copy(xi)
-        xi2[i_s] = 1
-        # Make new vertex list a hashable tuple
-        xi2_t = tuple(xi2)
-        # Append to cell
-        self.C0.add_vertex(self.V[xi2_t])
-        # Connect neighbors and vice versa
-        # Parent point
-        self.V[xi2_t].connect(self.V[xi_t])
-
-        # Connect all family of simplices in parent containers
-        for x_ip in x_parents:
-            self.V[xi2_t].connect(self.V[x_ip])
-
-        x_parents2 = copy.copy(x_parents)
-        x_parents2.append(xi_t)
-
-        i_s += 1
-        if i_s == self.dim:
-            return
-        # Permutate
-        self.perm_symmetry(i_s, x_parents2, xi2)
-
-    def add_centroid(self):
-        """Split the central edge between the origin and supremum of
-        a cell and add the new vertex to the complex"""
-        self.centroid = list(
-            (np.array(self.origin) + np.array(self.supremum)) / 2.0)
-        self.C0.add_vertex(self.V[tuple(self.centroid)])
-        self.C0.centroid = self.centroid
-
-        # Disconnect origin and supremum
-        self.V[tuple(self.origin)].disconnect(self.V[tuple(self.supremum)])
-
-        # Connect centroid to all other vertices
-        for v in self.C0():
-            self.V[tuple(self.centroid)].connect(self.V[tuple(v.x)])
-
-        self.centroid_added = True
-        return
-
-    # Construct incidence array:
-    def incidence(self):
-        if self.centroid_added:
-            self.structure = np.zeros([2 ** self.dim + 1, 2 ** self.dim + 1],
-                                         dtype=int)
-        else:
-            self.structure = np.zeros([2 ** self.dim, 2 ** self.dim],
-                                         dtype=int)
-
-        for v in self.HC.C0():
-            for v2 in v.nn:
-                self.structure[v.index, v2.index] = 1
-
-        return
-
-    # A more sparse incidence generator:
-    def graph_map(self):
-        """ Make a list of size 2**n + 1 where an entry is a vertex
-        incidence, each list element contains a list of indexes
-        corresponding to that entries neighbors"""
-
-        self.graph = [[v2.index for v2 in v.nn] for v in self.C0()]
-
-    # Graph structure method:
-    # 0. Capture the indices of the initial cell.
-    # 1. Generate new origin and supremum scalars based on current generation
-    # 2. Generate a new set of vertices corresponding to a new
-    #    "origin" and "supremum"
-    # 3. Connected based on the indices of the previous graph structure
-    # 4. Disconnect the edges in the original cell
-
-    def sub_generate_cell(self, C_i, gen):
-        """Subgenerate a cell `C_i` of generation `gen` and
-        homology group rank `hgr`."""
-        origin_new = tuple(C_i.centroid)
-        centroid_index = len(C_i()) - 1
-
-        # If not gen append
-        try:
-            self.H[gen]
-        except IndexError:
-            self.H.append([])
-
-        # Generate subcubes using every extreme vertex in C_i as a supremum
-        # and the centroid of C_i as the origin
-        H_new = []  # list storing all the new cubes split from C_i
-        for i, v in enumerate(C_i()[:-1]):
-            supremum = tuple(v.x)
-            H_new.append(
-                self.construct_hypercube(origin_new, supremum, gen, C_i.hg_n))
-
-        for i, connections in enumerate(self.graph):
-            # Present vertex V_new[i]; connect to all connections:
-            if i == centroid_index:  # Break out of centroid
-                break
-
-            for j in connections:
-                C_i()[i].disconnect(C_i()[j])
-
-        # Destroy the old cell
-        if C_i is not self.C0:  # Garbage collector does this anyway; not needed
-            del C_i
-
-        # TODO: Recalculate all the homology group ranks of each cell
-        return H_new
-
-    def split_generation(self):
-        """
-        Run sub_generate_cell for every cell in the current complex self.gen
-        """
-        no_splits = False  # USED IN SHGO
-        try:
-            for c in self.H[self.gen]:
-                if self.symmetry:
-                    # self.sub_generate_cell_symmetry(c, self.gen + 1)
-                    self.split_simplex_symmetry(c, self.gen + 1)
-                else:
-                    self.sub_generate_cell(c, self.gen + 1)
-        except IndexError:
-            no_splits = True  # USED IN SHGO
-
-        self.gen += 1
-        return no_splits  # USED IN SHGO
-
-    def construct_hypercube(self, origin, supremum, gen, hgr,
-                            printout=False):
-        """
-        Build a hypercube with triangulations symmetric to C0.
-
-        Parameters
-        ----------
-        origin : vec
-        supremum : vec (tuple)
-        gen : generation
-        hgr : parent homology group rank
-        """
-        # Initiate new cell
-        v_o = np.array(origin)
-        v_s = np.array(supremum)
-
-        C_new = Cell(gen, hgr, origin, supremum)
-        C_new.centroid = tuple((v_o + v_s) * .5)
-
-        # Build new indexed vertex list
-        V_new = []
-
-        for i, v in enumerate(self.C0()[:-1]):
-            v_x = np.array(v.x)
-            sub_cell_t1 = v_o - v_o * v_x
-            sub_cell_t2 = v_s * v_x
-
-            vec = sub_cell_t1 + sub_cell_t2
-
-            vec = tuple(vec)
-            C_new.add_vertex(self.V[vec])
-            V_new.append(vec)
-
-        # Add new centroid
-        C_new.add_vertex(self.V[C_new.centroid])
-        V_new.append(C_new.centroid)
-
-        # Connect new vertices #TODO: Thread into other loop; no need for V_new
-        for i, connections in enumerate(self.graph):
-            # Present vertex V_new[i]; connect to all connections:
-            for j in connections:
-                self.V[V_new[i]].connect(self.V[V_new[j]])
-
-        if printout:
-            print("A sub hyper cube with:")
-            print(f"origin: {origin}")
-            print(f"supremum: {supremum}")
-            for v in C_new():
-                v.print_out()
-
-        # Append the new cell to the to complex
-        self.H[gen].append(C_new)
-
-        return C_new
-
-    def split_simplex_symmetry(self, S, gen):
-        """
-        Split a hypersimplex S into two sub simplices by building a hyperplane
-        which connects to a new vertex on an edge (the longest edge in
-        dim = {2, 3}) and every other vertex in the simplex that is not
-        connected to the edge being split.
-
-        This function utilizes the knowledge that the problem is specified
-        with symmetric constraints
-
-        The longest edge is tracked by an ordering of the
-        vertices in every simplices, the edge between first and second
-        vertex is the longest edge to be split in the next iteration.
-        """
-        # If not gen append
-        try:
-            self.H[gen]
-        except IndexError:
-            self.H.append([])
-
-        # Find new vertex.
-        # V_new_x = tuple((np.array(C()[0].x) + np.array(C()[1].x)) / 2.0)
-        s = S()
-        firstx = s[0].x
-        lastx = s[-1].x
-        V_new = self.V[tuple((np.array(firstx) + np.array(lastx)) / 2.0)]
-
-        # Disconnect old longest edge
-        self.V[firstx].disconnect(self.V[lastx])
-
-        # Connect new vertices to all other vertices
-        for v in s[:]:
-            v.connect(self.V[V_new.x])
-
-        # New "lower" simplex
-        S_new_l = Simplex(gen, S.hg_n, self.generation_cycle,
-                          self.dim)
-        S_new_l.add_vertex(s[0])
-        S_new_l.add_vertex(V_new)  # Add new vertex
-        for v in s[1:-1]:  # Add all other vertices
-            S_new_l.add_vertex(v)
-
-        # New "upper" simplex
-        S_new_u = Simplex(gen, S.hg_n, S.generation_cycle, self.dim)
-
-        # First vertex on new long edge
-        S_new_u.add_vertex(s[S_new_u.generation_cycle + 1])
-
-        for v in s[1:-1]:  # Remaining vertices
-            S_new_u.add_vertex(v)
-
-        for k, v in enumerate(s[1:-1]):  # iterate through inner vertices
-            if k == S.generation_cycle:
-                S_new_u.add_vertex(V_new)
-            else:
-                S_new_u.add_vertex(v)
-
-        S_new_u.add_vertex(s[-1])  # Second vertex on new long edge
-
-        self.H[gen].append(S_new_l)
-        self.H[gen].append(S_new_u)
-
-        return
-
-    # Plots
-    def plot_complex(self):
-        """
-             Here, C is the LIST of simplexes S in the
-             2- or 3-D complex
-
-             To plot a single simplex S in a set C, use e.g., [C[0]]
-        """
-        from matplotlib import pyplot
-        if self.dim == 2:
-            pyplot.figure()
-            for C in self.H:
-                for c in C:
-                    for v in c():
-                        if self.bounds is None:
-                            x_a = np.array(v.x, dtype=float)
-                        else:
-                            x_a = np.array(v.x, dtype=float)
-                            for i in range(len(self.bounds)):
-                                x_a[i] = (x_a[i] * (self.bounds[i][1]
-                                                    - self.bounds[i][0])
-                                          + self.bounds[i][0])
-
-                        # logging.info('v.x_a = {}'.format(x_a))
-
-                        pyplot.plot([x_a[0]], [x_a[1]], 'o')
-
-                        xlines = []
-                        ylines = []
-                        for vn in v.nn:
-                            if self.bounds is None:
-                                xn_a = np.array(vn.x, dtype=float)
-                            else:
-                                xn_a = np.array(vn.x, dtype=float)
-                                for i in range(len(self.bounds)):
-                                    xn_a[i] = (xn_a[i] * (self.bounds[i][1]
-                                                          - self.bounds[i][0])
-                                               + self.bounds[i][0])
-
-                            # logging.info('vn.x = {}'.format(vn.x))
-
-                            xlines.append(xn_a[0])
-                            ylines.append(xn_a[1])
-                            xlines.append(x_a[0])
-                            ylines.append(x_a[1])
-
-                        pyplot.plot(xlines, ylines)
-
-            if self.bounds is None:
-                pyplot.ylim([-1e-2, 1 + 1e-2])
-                pyplot.xlim([-1e-2, 1 + 1e-2])
-            else:
-                pyplot.ylim(
-                    [self.bounds[1][0] - 1e-2, self.bounds[1][1] + 1e-2])
-                pyplot.xlim(
-                    [self.bounds[0][0] - 1e-2, self.bounds[0][1] + 1e-2])
-
-            pyplot.show()
-
-        elif self.dim == 3:
-            fig = pyplot.figure()
-            ax = fig.add_subplot(111, projection='3d')
-
-            for C in self.H:
-                for c in C:
-                    for v in c():
-                        x = []
-                        y = []
-                        z = []
-                        # logging.info('v.x = {}'.format(v.x))
-                        x.append(v.x[0])
-                        y.append(v.x[1])
-                        z.append(v.x[2])
-                        for vn in v.nn:
-                            x.append(vn.x[0])
-                            y.append(vn.x[1])
-                            z.append(vn.x[2])
-                            x.append(v.x[0])
-                            y.append(v.x[1])
-                            z.append(v.x[2])
-                            # logging.info('vn.x = {}'.format(vn.x))
-
-                        ax.plot(x, y, z, label='simplex')
-
-            pyplot.show()
-        else:
-            print("dimension higher than 3 or wrong complex format")
-        return
-
-
-class VertexGroup:
-    def __init__(self, p_gen, p_hgr):
-        self.p_gen = p_gen  # parent generation
-        self.p_hgr = p_hgr  # parent homology group rank
-        self.hg_n = None
-        self.hg_d = None
-
-        # Maybe add parent homology group rank total history
-        # This is the sum off all previously split cells
-        # cumulatively throughout its entire history
-        self.C = []
-
-    def __call__(self):
-        return self.C
-
-    def add_vertex(self, V):
-        if V not in self.C:
-            self.C.append(V)
-
-    def homology_group_rank(self):
-        """
-        Returns the homology group order of the current cell
-        """
-        if self.hg_n is None:
-            self.hg_n = sum(1 for v in self.C if v.minimiser())
-
-        return self.hg_n
-
-    def homology_group_differential(self):
-        """
-        Returns the difference between the current homology group of the
-        cell and its parent group
-        """
-        if self.hg_d is None:
-            self.hgd = self.hg_n - self.p_hgr
-
-        return self.hgd
-
-    def polytopial_sperner_lemma(self):
-        """
-        Returns the number of stationary points theoretically contained in the
-        cell based information currently known about the cell
-        """
-        pass
-
-    def print_out(self):
-        """
-        Print the current cell to console
-        """
-        for v in self():
-            v.print_out()
-
-
-class Cell(VertexGroup):
-    """
-    Contains a cell that is symmetric to the initial hypercube triangulation
-    """
-
-    def __init__(self, p_gen, p_hgr, origin, supremum):
-        super().__init__(p_gen, p_hgr)
-
-        self.origin = origin
-        self.supremum = supremum
-        self.centroid = None  # (Not always used)
-        # TODO: self.bounds
-
-
-class Simplex(VertexGroup):
-    """
-    Contains a simplex that is symmetric to the initial symmetry constrained
-    hypersimplex triangulation
-    """
-
-    def __init__(self, p_gen, p_hgr, generation_cycle, dim):
-        super().__init__(p_gen, p_hgr)
-
-        self.generation_cycle = (generation_cycle + 1) % (dim - 1)
-
-
-class Vertex:
-    def __init__(self, x, bounds=None, func=None, func_args=(), g_cons=None,
-                 g_cons_args=(), nn=None, index=None):
-        self.x = x
-        self.order = sum(x)
-        x_a = np.array(x, dtype=float)
-        if bounds is not None:
-            for i, (lb, ub) in enumerate(bounds):
-                x_a[i] = x_a[i] * (ub - lb) + lb
-
-        # TODO: Make saving the array structure optional
-        self.x_a = x_a
-
-        # Note Vertex is only initiated once for all x so only
-        # evaluated once
-        if func is not None:
-            self.feasible = True
-            if g_cons is not None:
-                for g, args in zip(g_cons, g_cons_args):
-                    if g(self.x_a, *args) < 0.0:
-                        self.f = np.inf
-                        self.feasible = False
-                        break
-            if self.feasible:
-                self.f = func(x_a, *func_args)
-
-        if nn is not None:
-            self.nn = nn
-        else:
-            self.nn = set()
-
-        self.fval = None
-        self.check_min = True
-
-        # Index:
-        if index is not None:
-            self.index = index
-
-    def __hash__(self):
-        return hash(self.x)
-
-    def connect(self, v):
-        if v is not self and v not in self.nn:
-            self.nn.add(v)
-            v.nn.add(self)
-
-            if self.minimiser():
-                v._min = False
-                v.check_min = False
-
-            # TEMPORARY
-            self.check_min = True
-            v.check_min = True
-
-    def disconnect(self, v):
-        if v in self.nn:
-            self.nn.remove(v)
-            v.nn.remove(self)
-            self.check_min = True
-            v.check_min = True
-
-    def minimiser(self):
-        """Check whether this vertex is strictly less than all its neighbors"""
-        if self.check_min:
-            self._min = all(self.f < v.f for v in self.nn)
-            self.check_min = False
-
-        return self._min
-
-    def print_out(self):
-        print(f"Vertex: {self.x}")
-        constr = 'Connections: '
-        for vc in self.nn:
-            constr += f'{vc.x} '
-
-        print(constr)
-        print(f'Order = {self.order}')
-
-
-class VertexCache:
-    def __init__(self, func, func_args=(), bounds=None, g_cons=None,
-                 g_cons_args=(), indexed=True):
-
-        self.cache = {}
-        self.func = func
-        self.g_cons = g_cons
-        self.g_cons_args = g_cons_args
-        self.func_args = func_args
-        self.bounds = bounds
-        self.nfev = 0
-        self.size = 0
-
-        if indexed:
-            self.index = -1
-
-    def __getitem__(self, x, indexed=True):
-        try:
-            return self.cache[x]
-        except KeyError:
-            if indexed:
-                self.index += 1
-                xval = Vertex(x, bounds=self.bounds,
-                              func=self.func, func_args=self.func_args,
-                              g_cons=self.g_cons,
-                              g_cons_args=self.g_cons_args,
-                              index=self.index)
-            else:
-                xval = Vertex(x, bounds=self.bounds,
-                              func=self.func, func_args=self.func_args,
-                              g_cons=self.g_cons,
-                              g_cons_args=self.g_cons_args)
-
-            # logging.info("New generated vertex at x = {}".format(x))
-            # NOTE: Surprisingly high performance increase if logging is commented out
-            self.cache[x] = xval
-
-            # TODO: Check
-            if self.func is not None:
-                if self.g_cons is not None:
-                    if xval.feasible:
-                        self.nfev += 1
-                        self.size += 1
-                    else:
-                        self.size += 1
-                else:
-                    self.nfev += 1
-                    self.size += 1
-
-            return self.cache[x]
diff --git a/scipy/optimize/tests/test__shgo.py b/scipy/optimize/tests/test__shgo.py
index 5ccaa20937d3..ef3d7d4503f7 100644
--- a/scipy/optimize/tests/test__shgo.py
+++ b/scipy/optimize/tests/test__shgo.py
@@ -1,9 +1,11 @@
 import logging
 import numpy
+import time
 from numpy.testing import assert_allclose
 import pytest
 from pytest import raises as assert_raises, warns
-from scipy.optimize import shgo, Bounds, minimize
+from scipy.optimize import (shgo, Bounds, minimize_scalar, minimize, rosen,
+                            rosen_der, rosen_hess)
 from scipy.optimize._shgo import SHGO
 
 
@@ -53,7 +55,7 @@ def g(x):
 
 class StructTest2(StructTestFunction):
     """
-    Scalar function with several minima to test all minimizer retrievals
+    Scalar function with several minima to test all minimiser retrievals
     """
 
     def f(self, x):
@@ -68,7 +70,8 @@ def g(x):
 test2_1 = StructTest2(bounds=[(0, 60)],
                       expected_x=[1.53567906],
                       expected_fun=-28.44677132,
-                      # Important: test that funl return is in the correct order
+                      # Important: test that funl return is in the correct
+                      # order
                       expected_xl=numpy.array([[1.53567906],
                                                [55.01782167],
                                                [7.80894889],
@@ -191,10 +194,12 @@ def f(self, x):
 
 class StructTestLJ(StructTestFunction):
     """
-    LennardJones objective function. Used to test symmetry constraints settings.
+    LennardJones objective function. Used to test symmetry constraints
+    settings.
     """
 
     def f(self, x, *args):
+        print(f'x = {x}')
         self.N = args[0]
         k = int(self.N / 3)
         s = 0.0
@@ -222,15 +227,31 @@ def f(self, x, *args):
 
 testLJ = StructTestLJ(bounds=boundsLJ,
                       expected_fun=[-1.0],
-                      expected_x=[-2.71247337e-08,
-                                  -2.71247337e-08,
-                                  -2.50000222e+00,
-                                  -2.71247337e-08,
-                                  -2.71247337e-08,
-                                  -1.50000222e+00]
+                      expected_x=None,
+                      # expected_x=[-2.71247337e-08,
+                      #            -2.71247337e-08,
+                      #            -2.50000222e+00,
+                      #            -2.71247337e-08,
+                      #            -2.71247337e-08,
+                      #            -1.50000222e+00]
                       )
 
 
+class StructTestS(StructTestFunction):
+    def f(self, x):
+        return ((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2
+                + (x[2] - 0.5) ** 2 + (x[3] - 0.5) ** 2)
+
+    g = None
+    cons = wrap_constraints(g)
+
+
+test_s = StructTestS(bounds=[(0, 2.0), ] * 4,
+                     expected_fun=0.0,
+                     expected_x=numpy.ones(4) - 0.5
+                     )
+
+
 class StructTestTable(StructTestFunction):
     def f(self, x):
         if x[0] == 3.0 and x[1] == 3.0:
@@ -277,16 +298,17 @@ def g4(x):
                                        )
 
 
-def run_test(test, args=(), test_atol=1e-5, n=128, iters=None,
+@pytest.mark.skip("Not a test")
+def run_test(test, args=(), test_atol=1e-5, n=100, iters=None,
              callback=None, minimizer_kwargs=None, options=None,
-             sampling_method='sobol'):
+             sampling_method='sobol', workers=1):
     res = shgo(test.f, test.bounds, args=args, constraints=test.cons,
                n=n, iters=iters, callback=callback,
                minimizer_kwargs=minimizer_kwargs, options=options,
-               sampling_method=sampling_method)
-
-    logging.info(res)
+               sampling_method=sampling_method, workers=workers)
 
+    print(f'res = {res}')
+    logging.info(f'res = {res}')
     if test.expected_x is not None:
         numpy.testing.assert_allclose(res.x, test.expected_x,
                                       rtol=test_atol,
@@ -313,7 +335,7 @@ def run_test(test, args=(), test_atol=1e-5, n=128, iters=None,
 # Base test functions:
 class TestShgoSobolTestFunctions:
     """
-    Global optimization tests with Sobol sampling:
+    Global optimisation tests with Sobol sampling:
     """
 
     # Sobol algorithm
@@ -330,7 +352,8 @@ def test_f1_2_sobol(self):
     def test_f1_3_sobol(self):
         """Multivariate test function 1:
         x[0]**2 + x[1]**2 with bounds=[(None, None),(None, None)]"""
-        run_test(test1_3)
+        options = {'disp': True}
+        run_test(test1_3, options=options)
 
     def test_f2_1_sobol(self):
         """Univariate test function on
@@ -348,29 +371,29 @@ def test_f3_sobol(self):
 
     @pytest.mark.slow
     def test_f4_sobol(self):
-        """NLP: (High-dimensional) Hock and Schittkowski 11 problem (HS11)"""
-        # run_test(test4_1, n=500)
-        # run_test(test4_1, n=800)
+        """NLP: (High dimensional) Hock and Schittkowski 11 problem (HS11)"""
         options = {'infty_constraints': False}
-        run_test(test4_1, n=2048, options=options)
+        # run_test(test4_1, n=990, options=options)
+        run_test(test4_1, n=990 * 2, options=options)
 
     def test_f5_1_sobol(self):
         """NLP: Eggholder, multimodal"""
-        run_test(test5_1, n=64)
+        # run_test(test5_1, n=30)
+        run_test(test5_1, n=60)
 
     def test_f5_2_sobol(self):
         """NLP: Eggholder, multimodal"""
         # run_test(test5_1, n=60, iters=5)
-        run_test(test5_1, n=128, iters=5)
+        run_test(test5_1, n=60, iters=5)
 
         # def test_t911(self):
-        #    """1-D tabletop function"""
+        #    """1D tabletop function"""
         #    run_test(test11_1)
 
 
 class TestShgoSimplicialTestFunctions:
     """
-    Global optimization tests with Simplicial sampling:
+    Global optimisation tests with Simplicial sampling:
     """
 
     def test_f1_1_simplicial(self):
@@ -386,13 +409,13 @@ def test_f1_2_simplicial(self):
     def test_f1_3_simplicial(self):
         """Multivariate test function 1: x[0]**2 + x[1]**2
         with bounds=[(None, None),(None, None)]"""
-        run_test(test1_3, n=1, sampling_method='simplicial')
+        run_test(test1_3, n=5, sampling_method='simplicial')
 
     def test_f2_1_simplicial(self):
         """Univariate test function on
         f(x) = (x - 30) * sin(x) with bounds=[(0, 60)]"""
         options = {'minimize_every_iter': False}
-        run_test(test2_1, iters=7, options=options,
+        run_test(test2_1, n=200, iters=7, options=options,
                  sampling_method='simplicial')
 
     def test_f2_2_simplicial(self):
@@ -406,18 +429,74 @@ def test_f3_simplicial(self):
 
     @pytest.mark.slow
     def test_f4_simplicial(self):
-        """NLP: (High-dimensional) Hock and Schittkowski 11 problem (HS11)"""
+        """NLP: (High dimensional) Hock and Schittkowski 11 problem (HS11)"""
         run_test(test4_1, n=1, sampling_method='simplicial')
 
-    def test_lj_symmetry(self):
+    def test_lj_symmetry_old(self):
         """LJ: Symmetry-constrained test function"""
         options = {'symmetry': True,
                    'disp': True}
         args = (6,)  # Number of atoms
-        run_test(testLJ, args=args, n=None,
-                 options=options, iters=4,
+        run_test(testLJ, args=args, n=300,
+                 options=options, iters=1,
                  sampling_method='simplicial')
 
+    def test_f5_1_lj_symmetry(self):
+        """LJ: Symmetry constrained test function"""
+        options = {'symmetry': [0, ] * 6,
+                   'disp': True}
+        args = (6,)  # No. of atoms
+
+        run_test(testLJ, args=args, n=300,
+                 options=options, iters=1,
+                 sampling_method='simplicial')
+
+    def test_f5_2_cons_symmetry(self):
+        """Symmetry constrained test function"""
+        options = {'symmetry': [0, 0],
+                   'disp': True}
+
+        run_test(test1_1, n=200,
+                 options=options, iters=1,
+                 sampling_method='simplicial')
+
+    def test_f5_3_cons_symmetry(self):
+        """Assymmetrically constrained test function"""
+        options = {'symmetry': [0, 0, 0, 3],
+                   'disp': True}
+
+        run_test(test_s, n=10000,
+                 options=options,
+                 iters=1,
+                 sampling_method='simplicial')
+
+    @pytest.mark.skip("Not a test")
+    def test_f0_min_variance(self):
+        """Return a minimum on a perfectly symmetric problem, based on
+            gh10429"""
+        avg = 0.5  # Given average value of x
+        cons = {'type': 'eq', 'fun': lambda x: numpy.mean(x) - avg}
+
+        # Minimize the variance of x under the given constraint
+        res = shgo(numpy.var, bounds=6 * [(0, 1)], constraints=cons)
+        assert res.success
+        assert_allclose(res.fun, 0, atol=1e-15)
+        assert_allclose(res.x, 0.5)
+
+    @pytest.mark.skip("Not a test")
+    def test_f0_min_variance_1D(self):
+        """Return a minimum on a perfectly symmetric 1D problem, based on
+            gh10538"""
+
+        def fun(x):
+            return x * (x - 1.0) * (x - 0.5)
+
+        bounds = [(0, 1)]
+        res = shgo(fun, bounds=bounds)
+        ref = minimize_scalar(fun, bounds=bounds[0])
+        assert res.success
+        assert_allclose(res.fun, ref.fun)
+        assert_allclose(res.x, ref.x, rtol=1e-6)
 
 # Argument test functions
 class TestShgoArguments:
@@ -428,7 +507,7 @@ def test_1_1_simpl_iter(self):
     def test_1_2_simpl_iter(self):
         """Iterative simplicial on TestFunction 2 (univariate)"""
         options = {'minimize_every_iter': False}
-        run_test(test2_1, n=None, iters=7, options=options,
+        run_test(test2_1, n=None, iters=9, options=options,
                  sampling_method='simplicial')
 
     def test_2_1_sobol_iter(self):
@@ -445,7 +524,8 @@ def test_2_2_sobol_iter(self):
         numpy.testing.assert_allclose(res.fun, test2_1.expected_fun, atol=1e-5)
 
     def test_3_1_disp_simplicial(self):
-        """Iterative sampling on TestFunction 1 and 2  (multi- and univariate)"""
+        """Iterative sampling on TestFunction 1 and 2  (multi and univariate)
+        """
 
         def callback_func(x):
             print("Local minimization callback test")
@@ -458,7 +538,7 @@ def callback_func(x):
                  callback=callback_func, options={'disp': True})
 
     def test_3_2_disp_sobol(self):
-        """Iterative sampling on TestFunction 1 and 2 (multi- and univariate)"""
+        """Iterative sampling on TestFunction 1 and 2 (multi and univariate)"""
 
         def callback_func(x):
             print("Local minimization callback test")
@@ -471,9 +551,11 @@ def callback_func(x):
                  callback=callback_func, options={'disp': True})
 
     def test_args_gh14589(self):
-        # Using `args` used to cause `shgo` to fail; see #14589, #15986, #16506
-        res = shgo(func=lambda x, y, z: x*z + y, bounds=[(0, 3)], args=(1, 2))
-        ref = shgo(func=lambda x: 2*x + 1, bounds=[(0, 3)])
+        """Using `args` used to cause `shgo` to fail; see #14589, #15986,
+        #16506"""
+        res = shgo(func=lambda x, y, z: x * z + y, bounds=[(0, 3)], args=(1, 2)
+                   )
+        ref = shgo(func=lambda x: 2 * x + 1, bounds=[(0, 3)])
         assert_allclose(res.fun, ref.fun)
         assert_allclose(res.x, ref.x)
 
@@ -501,22 +583,8 @@ def test_4_2_known_f_min(self):
         run_test(test4_1, n=None, test_atol=1e-5, options=options,
                  sampling_method='simplicial')
 
-    @pytest.mark.slow
-    def test_4_3_known_f_min(self):
-        """Test Global mode limiting local evalutions"""
-        options = {  # Specify known function value
-            'f_min': test4_1.expected_fun,
-            'f_tol': 1e-6,
-            # Specify number of local iterations to perform+
-            'minimize_every_iter': True,
-            'local_iter': 1,
-            'infty_constraints': False}
-
-        run_test(test4_1, n=1024, test_atol=1e-5, options=options,
-                 sampling_method='sobol')
-
     def test_4_4_known_f_min(self):
-        """Test Global mode limiting local evalutions for 1-D functions"""
+        """Test Global mode limiting local evaluations for 1D funcs"""
         options = {  # Specify known function value
             'f_min': test2_1.expected_fun,
             'f_tol': 1e-6,
@@ -570,10 +638,9 @@ def test_7_1_minkwargs(self):
             # unittests which run test4_1 normally
             minimizer_kwargs = {'method': solver,
                                 'constraints': test3_1.cons}
-            print(f"Solver = {solver}")
-            print("=" * 100)
-            run_test(test3_1, n=128, test_atol=1e-3,
-                     minimizer_kwargs=minimizer_kwargs, sampling_method='sobol')
+            run_test(test3_1, n=100, test_atol=1e-3,
+                     minimizer_kwargs=minimizer_kwargs,
+                     sampling_method='sobol')
 
     def test_7_2_minkwargs(self):
         """Test the minimizer_kwargs default inits"""
@@ -598,8 +665,9 @@ def hess(x):
                                 'hess': hess}
             logging.info(f"Solver = {solver}")
             logging.info("=" * 100)
-            run_test(test1_1, n=128, test_atol=1e-3,
-                     minimizer_kwargs=minimizer_kwargs, sampling_method='sobol')
+            run_test(test1_1, n=100, test_atol=1e-3,
+                     minimizer_kwargs=minimizer_kwargs,
+                     sampling_method='sobol')
 
     def test_8_homology_group_diff(self):
         options = {'minhgrd': 1,
@@ -615,27 +683,52 @@ def test_9_cons_g(self):
     def test_10_finite_time(self):
         """Test single function constraint passing"""
         options = {'maxtime': 1e-15}
-        shgo(test1_1.f, test1_1.bounds, n=1, iters=None,
-             options=options, sampling_method='sobol')
 
-    def test_11_f_min_time(self):
+        def f(x):
+            time.sleep(1e-14)
+            return 0.0
+
+        res = shgo(f, test1_1.bounds, iters=5, options=options)
+        # Assert that only 1 rather than 5 requested iterations ran:
+        assert res.nit == 1
+
+    def test_11_f_min_0(self):
         """Test to cover the case where f_lowest == 0"""
-        options = {'maxtime': 1e-15,
-                   'f_min': 0.0}
-        shgo(test1_2.f, test1_2.bounds, n=1, iters=None,
-             options=options, sampling_method='sobol')
+        options = {'f_min': 0.0,
+                   'disp': True}
+        res = shgo(test1_2.f, test1_2.bounds, n=10, iters=None,
+                   options=options, sampling_method='sobol')
+        numpy.testing.assert_equal(0, res.x[0])
+        numpy.testing.assert_equal(0, res.x[1])
 
+    # @nottest
+    @pytest.mark.skip(reason="no way of currently testing this")
     def test_12_sobol_inf_cons(self):
         """Test to cover the case where f_lowest == 0"""
+        # TODO: This test doesn't cover anything new, it is unknown what the
+        # original test was intended for as it was never complete. Delete or
+        # replace in the future.
         options = {'maxtime': 1e-15,
                    'f_min': 0.0}
-        shgo(test1_2.f, test1_2.bounds, n=1, iters=None,
-             options=options, sampling_method='sobol')
+        res = shgo(test1_2.f, test1_2.bounds, n=1, iters=None,
+                   options=options, sampling_method='sobol')
+        numpy.testing.assert_equal(0.0, res.fun)
+
+    def test_13_high_sobol(self):
+        """Test init of high-dimensional sobol sequences"""
+
+        def f(x):
+            return 0
+
+        bounds = [(None, None), ] * 41
+        SHGOc = SHGO(f, bounds, sampling_method='sobol')
+        # SHGOc.sobol_points(2, 50)
+        SHGOc.sampling_function(2, 50)
 
     def test_14_local_iter(self):
         """Test limited local iterations for a pseudo-global mode"""
         options = {'local_iter': 4}
-        run_test(test5_1, n=64, options=options)
+        run_test(test5_1, n=60, options=options)
 
     def test_15_min_every_iter(self):
         """Test minimize every iter options and cover function cache"""
@@ -652,8 +745,9 @@ def test_16_disp_bounds_minimizer(self):
 
     def test_17_custom_sampling(self):
         """Test the functionality to add custom sampling methods to shgo"""
+
         def sample(n, d):
-            return numpy.random.uniform(size=(n,d))
+            return numpy.random.uniform(size=(n, d))
 
         run_test(test1_1, n=30, sampling_method=sample)
 
@@ -678,17 +772,154 @@ def f(x):
         numpy.testing.assert_allclose(res_new_bounds.x,
                                       res_old_bounds.x)
 
+    def test_19_parallelization(self):
+        """Test the functionality to add custom sampling methods to shgo"""
+
+        run_test(test1_1, n=30, workers=1)  # Constrained
+        run_test(test_s, n=30, workers=1)  # Unconstrained
+
+    def test_20_constrained_args(self):
+        """Test that constraints can be passed to arguments"""
+
+        def eggholder(x):
+            return (-(x[1] + 47.0)
+                    * numpy.sin(numpy.sqrt(abs(x[0] / 2.0 + (x[1] + 47.0))))
+                    - x[0] * numpy.sin(numpy.sqrt(abs(x[0] - (x[1] + 47.0))))
+                    )
+
+        def f(x):  # (cattle-feed)
+            return 24.55 * x[0] + 26.75 * x[1] + 39 * x[2] + 40.50 * x[3]
+
+        bounds = [(0, 1.0), ] * 4
+
+        def g1_modified(x, i):
+            return i * 2.3 * x[0] + i * 5.6 * x[1] + 11.1 * x[2] + 1.3 * x[
+                3] - 5  # >=0
+
+        def g2(x):
+            return (12 * x[0] + 11.9 * x[1] + 41.8 * x[2] + 52.1 * x[3] - 21
+                    - 1.645 * numpy.sqrt(0.28 * x[0] ** 2 + 0.19 * x[1] ** 2
+                                         + 20.5 * x[2] ** 2 + 0.62 * x[3] ** 2)
+                    )  # >=0
+
+        def h1(x):
+            return x[0] + x[1] + x[2] + x[3] - 1  # == 0
+
+        cons = ({'type': 'ineq', 'fun': g1_modified, "args": (0,)},
+                {'type': 'ineq', 'fun': g2},
+                {'type': 'eq', 'fun': h1})
+
+        shgo(f, bounds, n=300, iters=1, constraints=cons)
+        # using constrain with arguments AND sampling method sobol
+        shgo(f, bounds, n=300, iters=1, constraints=cons,
+             sampling_method='sobol')
+
+    def test_21_1_jac_true(self):
+        """Test that shgo can handle objective functions that return the
+        gradient alongside the objective value. Fixes gh-13547"""
+        # previous
+        def func(x):
+            return numpy.sum(numpy.power(x, 2)), 2 * x
+
+        shgo(
+            func,
+            bounds=[[-1, 1], [1, 2]],
+            n=100, iters=5,
+            sampling_method="sobol",
+            minimizer_kwargs={'method': 'SLSQP', 'jac': True}
+        )
+
+        # new
+        def func(x):
+            return numpy.sum(x ** 2), 2 * x
+
+        bounds = [[-1, 1], [1, 2], [-1, 1], [1, 2], [0, 3]]
+
+        res = shgo(func, bounds=bounds, sampling_method="sobol",
+                   minimizer_kwargs={'method': 'SLSQP', 'jac': True})
+        ref = minimize(func, x0=[1, 1, 1, 1, 1], bounds=bounds,
+                       jac=True)
+        assert res.success
+        assert_allclose(res.fun, ref.fun)
+        assert_allclose(res.x, ref.x, atol=1e-15)
+
+    @pytest.mark.parametrize('derivative', ['jac', 'hess', 'hessp'])
+    def test_21_2_derivative_options(self, derivative):
+        """shgo used to raise an error when passing `options` with 'jac'
+        # see gh-12963. check that this is resolved
+        """
+
+        def objective(x):
+            return 3 * x[0] * x[0] + 2 * x[0] + 5
+
+        def gradient(x):
+            return 6 * x[0] + 2
+
+        def hess(x):
+            return 6
+
+        def hessp(x, p):
+            return 6 * p
+
+        derivative_funcs = {'jac': gradient, 'hess': hess, 'hessp': hessp}
+        options = {derivative: derivative_funcs[derivative]}
+        minimizer_kwargs = {'method': 'trust-constr'}
+
+        bounds = [(-100, 100)]
+        res = shgo(objective, bounds, minimizer_kwargs=minimizer_kwargs,
+                   options=options)
+        ref = minimize(objective, x0=[0], bounds=bounds, **minimizer_kwargs,
+                       **options)
+
+        assert res.success
+        numpy.testing.assert_allclose(res.fun, ref.fun)
+        numpy.testing.assert_allclose(res.x, ref.x)
+
+    def test_21_3_hess_options_rosen(self):
+        """Ensure the Hessian gets passed correctly to the local minimizer
+        routine. Previous report gh-14533.
+        """
+        bounds = [(0, 1.6), (0, 1.6), (0, 1.4), (0, 1.4), (0, 1.4)]
+        options = {'jac': rosen_der, 'hess': rosen_hess}
+        minimizer_kwargs = {'method': 'Newton-CG'}
+        res = shgo(rosen, bounds, minimizer_kwargs=minimizer_kwargs,
+                   options=options)
+        ref = minimize(rosen, numpy.zeros(5), method='Newton-CG',
+                       **options)
+        assert res.success
+        assert_allclose(res.fun, ref.fun)
+        assert_allclose(res.x, ref.x, atol=1e-15)
+
+    def test_21_arg_tuple_sobol(self):
+        """shgo used to raise an error when passing `args` with Sobol sampling
+        # see gh-12114. check that this is resolved"""
+
+        def fun(x, k):
+            return x[0] ** k
+
+        constraints = ({'type': 'ineq', 'fun': lambda x: x[0] - 1})
+
+        bounds = [(0, 10)]
+        res = shgo(fun, bounds, args=(1,), constraints=constraints,
+                   sampling_method='sobol')
+        ref = minimize(fun, numpy.zeros(1), bounds=bounds, args=(1,),
+                       constraints=constraints)
+        assert res.success
+        assert_allclose(res.fun, ref.fun)
+        assert_allclose(res.x, ref.x)
+
 
 # Failure test functions
 class TestShgoFailures:
     def test_1_maxiter(self):
         """Test failure on insufficient iterations"""
         options = {'maxiter': 2}
-        res = shgo(test4_1.f, test4_1.bounds, n=4, iters=None,
+        res = shgo(test4_1.f, test4_1.bounds, n=2, iters=None,
                    options=options, sampling_method='sobol')
 
         numpy.testing.assert_equal(False, res.success)
-        numpy.testing.assert_equal(4, res.nfev)
+        # numpy.testing.assert_equal(4, res.nfev)
+        numpy.testing.assert_equal(4, res.tnev)
 
     def test_2_sampling(self):
         """Rejection of unknown sampling method"""
@@ -699,12 +930,13 @@ def test_3_1_no_min_pool_sobol(self):
         """Check that the routine stops when no minimiser is found
            after maximum specified function evaluations"""
         options = {'maxfev': 10,
+                   # 'maxev': 10,
                    'disp': True}
-        res = shgo(test_table.f, test_table.bounds, n=4, options=options,
+        res = shgo(test_table.f, test_table.bounds, n=3, options=options,
                    sampling_method='sobol')
         numpy.testing.assert_equal(False, res.success)
-
-        numpy.testing.assert_equal(16, res.nfev)
+        # numpy.testing.assert_equal(9, res.nfev)
+        numpy.testing.assert_equal(12, res.nfev)
 
     def test_3_2_no_min_pool_simplicial(self):
         """Check that the routine stops when no minimiser is found
@@ -728,11 +960,11 @@ def test_4_2_bound_err(self):
     def test_5_1_1_infeasible_sobol(self):
         """Ensures the algorithm terminates on infeasible problems
            after maxev is exceeded. Use infty constraints option"""
-        options = {'maxev': 64,
+        options = {'maxev': 100,
                    'disp': True}
 
         res = shgo(test_infeasible.f, test_infeasible.bounds,
-                   constraints=test_infeasible.cons, n=64, options=options,
+                   constraints=test_infeasible.cons, n=100, options=options,
                    sampling_method='sobol')
 
         numpy.testing.assert_equal(False, res.success)
@@ -740,12 +972,12 @@ def test_5_1_1_infeasible_sobol(self):
     def test_5_1_2_infeasible_sobol(self):
         """Ensures the algorithm terminates on infeasible problems
            after maxev is exceeded. Do not use infty constraints option"""
-        options = {'maxev': 64,
+        options = {'maxev': 100,
                    'disp': True,
                    'infty_constraints': False}
 
         res = shgo(test_infeasible.f, test_infeasible.bounds,
-                   constraints=test_infeasible.cons, n=64, options=options,
+                   constraints=test_infeasible.cons, n=100, options=options,
                    sampling_method='sobol')
 
         numpy.testing.assert_equal(False, res.success)
@@ -763,7 +995,7 @@ def test_5_2_infeasible_simplicial(self):
         numpy.testing.assert_equal(False, res.success)
 
     def test_6_1_lower_known_f_min(self):
-        """Test Global mode limiting local evalutions with f* too high"""
+        """Test Global mode limiting local evaluations with f* too high"""
         options = {  # Specify known function value
             'f_min': test2_1.expected_fun + 2.0,
             'f_tol': 1e-6,
@@ -780,33 +1012,42 @@ def test_6_1_lower_known_f_min(self):
                   }
         warns(UserWarning, shgo, *args, **kwargs)
 
-    @pytest.mark.parametrize('derivative', ['jac', 'hess', 'hessp'])
-    def test_21_2_derivative_options(self, derivative):
-        """shgo used to raise an error when passing `options` with 'jac'
-        # see gh-12829. check that this is resolved
-        """
-        def objective(x):
-            return 3 * x[0] * x[0] + 2 * x[0] + 5
+    def test(self):
+        from scipy.optimize import rosen, shgo
+        bounds = [(0, 2), (0, 2), (0, 2), (0, 2), (0, 2)]
 
-        def gradient(x):
-            return 6 * x[0] + 2
+        def fun(x):
+            fun.nfev += 1
+            return rosen(x)
 
-        def hess(x):
-            return 6
+        fun.nfev = 0
 
-        def hessp(x, p):
-            return 6 * p
+        result = shgo(fun, bounds)
+        print(result.x, result.fun, fun.nfev)  # 50
 
-        derivative_funcs = {'jac': gradient, 'hess': hess, 'hessp': hessp}
-        options = {derivative: derivative_funcs[derivative]}
-        minimizer_kwargs = {'method': 'trust-constr'}
 
-        bounds = [(-100, 100)]
-        res = shgo(objective, bounds, minimizer_kwargs=minimizer_kwargs,
-                   options=options)
-        ref = minimize(objective, x0=[0], bounds=bounds, **minimizer_kwargs,
-                       **options)
+# Returns
+class TestShgoReturns:
+    def test_1_nfev_simplicial(self):
+        bounds = [(0, 2), (0, 2), (0, 2), (0, 2), (0, 2)]
 
-        assert res.success
-        numpy.testing.assert_allclose(res.fun, ref.fun)
-        numpy.testing.assert_allclose(res.x, ref.x)
+        def fun(x):
+            fun.nfev += 1
+            return rosen(x)
+
+        fun.nfev = 0
+
+        result = shgo(fun, bounds)
+        numpy.testing.assert_equal(fun.nfev, result.nfev)
+
+    def test_1_nfev_sobol(self):
+        bounds = [(0, 2), (0, 2), (0, 2), (0, 2), (0, 2)]
+
+        def fun(x):
+            fun.nfev += 1
+            return rosen(x)
+
+        fun.nfev = 0
+
+        result = shgo(fun, bounds, sampling_method='sobol')
+        numpy.testing.assert_equal(fun.nfev, result.nfev)


Code review comments:
diff_hunk: @@ -422,7 +422,7 @@
                            LinearConstraint,
                            Bounds)
 from ._hessian_update_strategy import HessianUpdateStrategy, BFGS, SR1
-from ._shgo import shgo
+from ._shgo import shgo, SHGO
path: scipy/optimize/__init__.py
text: @user1:
```suggestion
from ._shgo import shgo
```
This is causing the final refguide failure, I think. What was the intent here? I'll remove it if I find locally that it's not being imported anywhere.

@author:
This is leftover from another branch update where `SHGO.HC.save_complex` and `SHGO.HC.load_complex` is used to save and load the progress of the optimization run to and from the hard disk drive, which is needed on clusters where user runtime is limited.

Currently there aren't many users that use this feature yet (and the code is not well tested), so I will make a separate PR in the future, it should also be possible to make a cleaner implementation for the `shgo` function to avoid a second class object in the namespace.