stefan-endres / shgo Goto Github PK
View Code? Open in Web Editor NEWSimplicial Homology Global Optimization
Home Page: https://stefan-endres.github.io/shgo
License: MIT License
Simplicial Homology Global Optimization
Home Page: https://stefan-endres.github.io/shgo
License: MIT License
I think it's fine to run only on Python 3, but currently the code fails on 3.5, which is still quite a large installed base. We can use Travis or Tox to test for different version compatibility.
Force round-off of a local minimisation of the smallest sampling value every iteration, even in the case where no minimiser homology is found. This solves the issue of the symmetric sampling method not finding solutions for symmetrically structured problems such as the issue raised in the scipy repository.
from shgo import shgo
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))))
- x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47.0))))
)
def f(x): # (cattle-feed)
return 24.55x[0] + 26.75x[1] + 39x[2] + 40.50x[3]
def g1(x):
return 2.3x[0] + 5.6x[1] + 11.1x[2] + 1.3x[3] - 5 # >=0
def g2(x):
return (12x[0] + 11.9x[1] +41.8x[2] + 52.1x[3] - 21
- 1.645 * np.sqrt(0.28x[0]**2 + 0.19x[1]**2
+ 20.5x[2]**2 + 0.62x[3]**2)
) # >=0
def h1(x):
return x[0] + x[1] + x[2] + x[3] - 1 # == 0
cons = ({'type': 'ineq', 'fun': g1},
{'type': 'ineq', 'fun': g2},
{'type': 'eq', 'fun': h1})
bounds = [(0, 1.0),]*4
res = shgo(f, bounds, iters=3, constraints=cons)
res
##################################################
#################### modified ####################
#using sampling method sobol
res_sobol = shgo(f, bounds, iters=3, constraints=cons,sampling_method='sobol')
res_sobol # still working
#using constrain with arguments
def g1_modified(x,i):
return i2.3x[0] + i5.6x[1] + 11.1x[2] + 1.3x[3] - 5 # >=0
cons = ({'type': 'ineq', 'fun': g1_modified, "args":(0,)},
{'type': 'ineq', 'fun': g2},
{'type': 'eq', 'fun': h1})
res_modConstr = shgo(f, bounds, iters=3, constraints=cons)
res_modConstr #still working
#using constrain with arguments AND sampling method sobol
res_modConstr_sobol = shgo(f, bounds, iters=3, constraints=cons,sampling_method='sobol')
res_modConstr_sobol #TypeError: g1_modified() missing 1 required positional argument: 'i'
There are functions in sobol_seq.py which are not used.
is_prime()
is only used in prime_ge
, which is not used anywhere else.i4_uniform()
is not used anywhere elseThe less code we have, the less we need to maintain.
pytest is superior to UnitTest in almost every way. SciPy is also moving their test suite to pytest in the 1.0.0 release
Instead of the pattern of
if options is not None:
if 'f_tol' in options:
self.minimizer_kwargs['options']['ftol'] = options['f_tol']
...
I think we can rather do
if options is not None:
self.minimizer_kwargs.update(options)
This will automatically overwrite keys if they exist in options. Less code to maintain.
The 1-dimensional Sobol sequences can't be implemented as Delaunay triangulation object since qhull only accepts dimensions 2 and higher.
This is causing major issues in refactoring and expansion of the code since we have to deal with two data structures in all the functions.
A common data structure for both cases similar to the simplicial sampling method would ease future changes to the code.
Currently in the go_bench_f_min.py
file uses the interruptingcow module (interruptingcow.timeout
) to stop the entire while loop in finite time. We cannot asses the time pass at the start of every loop since building on the simplicial complex (by either using delaunay triangulation or the hypercube sampling) can take extremely long on a single iteration.
We want to minimise dependencies for shgo since we want to include it in many other modules.
Therefore we need a way to stop the entire loop in a similar way that interruptingcow by wrapping the entire loop in a with statement and raising a RuntimeError when the processing time limit expires.
It seems that the equality constraints, referenced in the documentation as eq
, are not support. See https://github.com/Stefan-Endres/shgo/blob/master/shgo/_shgo.py#L495.
for cons in constraints:
if cons['type'] is '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)
If I've missed how to specify, please let me know.
The universal_iters branch is a major restructuring of the code with the following benefits:
sampling_method
argument.Unfortunately, it is not a simple refactoring since the inputs (n
and iters
arguments) actually need to be changed slightly to produce the same outputs in the overall routine, although the same behaviour should (at least theoretically) be possible to replicate given the correct usage of these parameters. The 3 unittests for ambiguous unittests also need to be removed since they are no longer ambiguous, but can used in a customized optimisation run. All other unittests are kept and working with little to no modifications.
However, just because the unittests are working does not mean the performance has been maintained. Therefore we need a way to compare performance of the two branches in addition to replicating the performance shown in publications. One possibility is adding unittests assessing the function evaluations of the current master branch. Finally some of the changes allowing for Python compatibility other than Python 3.6 which were committed to the master branch need to be rechecked with the new code.
I've noticed that {'maxfev':50} often results in >200 actual function evaluations. Is there a way to ensure this number is closer to what one desires, for example when the function to be evaluated is very expensive?
README.rst should be split into files at the top level headings
It seems plausible that an SHGO variant might be capable of finding good plateaus, which is to say relatively stable local minima where there is no particularly steep dropoff nearby. I realize that isn't a rigorous definition, but do you have thoughts on how this might be hacked to achieve the intent, approximately?
Not all tests are running. I receive the following warnings:
shgo/tests/test__shgo.py::TestFunction
cannot collect test class 'TestFunction' because it has a __init__ constructor
shgo/tests/test__shgo.py::Test1
cannot collect test class 'Test1' because it has a __init__ constructor
shgo/tests/test__shgo.py::Test2
cannot collect test class 'Test2' because it has a __init__ constructor
shgo/tests/test__shgo.py::Test3
cannot collect test class 'Test3' because it has a __init__ constructor
shgo/tests/test__shgo.py::Test4
cannot collect test class 'Test4' because it has a __init__ constructor
shgo/tests/test__shgo.py::TestLJ
cannot collect test class 'TestLJ' because it has a __init__ constructor
shgo/tests/test__shgo.py::TestTable
cannot collect test class 'TestTable' because it has a __init__ constructor
shgo/tests/test__shgo.py::TestInfeasible
cannot collect test class 'TestInfeasible' because it has a __init__ constructor
-- Docs: http://doc.pytest.org/en/latest/warnings.html
Configure codeclimate to run test coverage on the repo.
https://docs.codeclimate.com/v1.0/docs/configuring-test-coverage
We should not be checking for explicit type membership, but rather use isinstance
.
Incorporate the FunctionCache class from pointclasses.py into the main code. While the global routine will not evaluate the same value twice (and it is very unlikely that a non-deterministic local minimization would), it will allow us to reduce some complexity in the random sampling code to make it more similar to the simplicial sampling.
Is there an option to specify starting point? Something like this: res = shgo(f, bounds, x0=x0)
where the x0 is the user best guess for the optimum.
Some local minima are missing even if setting a very large n
. Is this the method's limit or the improper settings of parameters? Is it possible to find all local minima with some scale through this method?
e.g. the rastrigin function
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import shgo
def rastrigin(X):
return 20. + sum([(x**2 - 10. * np.cos(2. * np.pi * x)) for x in X])
func = rastrigin
xmin, xmax = -5.12, 5.12
x = np.linspace(xmin, xmax, 100)
y = np.linspace(xmin, xmax, 100)
xx, yy = np.meshgrid(x, y)
zz = func((xx, yy))
bounds = [(xmin, xmax)] * 2
ret = shgo(func, bounds,
n=100000, iters=1, sampling_method='sobol')
print('nfev: {:5d}'.format(ret.nfev))
print('nxl: {:5d}'.format(len(ret.xl)))
plt.figure()
plt.pcolormesh(x, y, zz)
for x1, y1 in ret.xl:
plt.plot(x1, y1, 'ro')
plt.show()
Coveralls appears to be tracking only the coverage of the test script itself, not the other files. This is why we had such great coverage to begin with and why it has dropped a lot now, since the tests are not all running (see #13).
We need to make it so that the coverage is tracked across all files in the project not just the test script.
There are many sub-routines in shgo that are low hanging fruit for parallelization, most importantly the sampling mapping of the objective function itself (in cases where it is possible to parallelize the problem). In addition we can optimize the sampling size of an iteration based on the computational resources available.
We need to be careful with both dependencies and code structure changes we include for several reasons. First our dependencies on scipy.optimize.minimize
and scipy.spatial.Delaunay
. Secondly our ambition to include shgo
in scipy.optimize
means it should ideally have the same dependencies and structure. Finally we want to minimize reliance on maintenance from other packages which can lead to such issues that we had with using multiprocessing_on_dill
in tgo
.
My suggestion is to use numba to avoid needing to change the code structure at all. We can do this using tricks such as the one used in poliastro:
https://github.com/poliastro/poliastro/blob/0.6.x/src/poliastro/jit.py
https://github.com/poliastro/poliastro/blob/master/setup.py#L40
which simply maps the decorator to the dependency if it is installed or does nothing. Obviously it is possible that we can use the same tricks for other libraries and methods by redefining range functions etc in our code.
While it is possible that SciPy will eventually include numba as a dependency, based on discussions in the scipy-dev mailing lists this will not happen in the near future:
https://mail.python.org/pipermail/scipy-dev/2018-March/022576.html
However, for now we should be able to maintain numba as an optional dependency as describe in the rest of this post. My idea is to provide two main modes of parallelization, based on both CPU and GPU parallelization. So therefore using numba for GPU parallelization would be ideal since it avoids extra dependencies. Finally numba provides us with access to LLVM that can be used in the sampling generation.
User architecture can be found or specified and then we can map it to our own decorator.
We can use @numba.cuda.jit
, I propose an early test of the objective function so that the user can be warned if this fails.
https://devblogs.nvidia.com/seven-things-numba/
https://numba.pydata.org/numba-doc/latest/cuda/index.html
We can use @hsa.jit(device=True)
https://numba.pydata.org/numba-doc/latest/hsa/overview.html
Our options for CPU parallelization
http://numba.pydata.org/numba-doc/dev/user/parallel.html
or multiprocessing_on_dill etc.
However, it appears that parallelization with numba isn't as simple as just adding @jit(parallel = True)
https://stackoverflow.com/questions/45610292/how-to-make-numba-jit-use-all-cpu-cores-parallelize-numba-jit
So we should also do a few tests on non-trivial functions to see if it is worth implementing.
There are a couple of places in the code where we use is
to check for string equality, so
if a is 'teststring':
# whatever
This is not guaranteed to work for all strings as it relies on an implementation detail in cpython (where small strings are cached. It should be taken away wherever it is used.
Fix the typos in the documentation found here and improve the readability as well as detailing the use of the arguments inside the algorithm #27.
The examples should also include usage of advanced features such as the symmetry argument of the algorithm to increase visibility.
Hi,
The function I would like to optimize has (amongst others) one parameter that only allows a boolean value (so the parameter can be zero or one).
How could I set this as a constraint for shgo?
Thanks in advance!
I'm trying to use shgo
to tune the gains of a control algorithm in a black-box fashion, by iteratively running the controller in a simulation (each simulation only takes about 20 seconds to complete) and computing a cost from the simulation results, which can then be used as the objective function f
for shgo
. However, having left shgo
running for the past few hours, I'm not seeing any progress being made. What I have noticed is the following:
shgo
seems to be stuck switching the values in the optimization parameter vector x
between the bounds x_min
and x_max
that were given to it, and doesn't seem to have done any searching within the bounds themselves.
In the first few minutes of running shgo
, it iterates quickly through each simulation run, however, as time goes on it slows down substantially, spending minutes between each iteration. What is shgo
doing in the minutes between simulation runs, is there a way to produce verbose output to investigate further?
Due to the strange results seen above, there could be two likely culprits that initially come to mind. The first is that my optimization vector x
is of dimension 15
, and after reading the shgo
documentation it looks like anything above 10
dimensions will be a struggle for shgo
to solve quickly (due to symmetries in the system dynamics, I can likely get this number below 10
though). The second potential problem is that the objective function (i.e. the cost output from the simulation) can return inf
(or sys.float_info.max
) if the controller causes the system that it is controlling to become unstable and have its states blow up, can shgo
deal with objective functions that return inf
?
Hi!
Is it correct that you overwrite the symmetry passed?
if self.symmetry:
self.symmetry = [0, ]*len(self.bounds)
The code expects to pass a skip argument to sobol_points
, but sobol_points_10k
doesn't accept such an argument.
This shows two things:
I understand why you have left the code there, but if it's not working it shouldn't ship effectively commented out this way. It's better to move this stuff to another branch where it is active and fix it there.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.