Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] MAINT: Update optimize.shgo #14013

Closed
wants to merge 1,236 commits into from
Closed

[WIP] MAINT: Update optimize.shgo #14013

wants to merge 1,236 commits into from

Conversation

Stefan-Endres
Copy link
Contributor

@Stefan-Endres Stefan-Endres commented May 9, 2021

Reference issues/fixes

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 ENH: optimize: Class based Optimizers #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 (resolves scipy.optimize.shgo fails on one-dimensional objective functions with particular symmetry #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)

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

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:

Proposed reviewers

@tupui @andyfaff

@tupui
Copy link
Member

tupui commented May 10, 2021

Thank you for all these improvements! I will try to have a look on it fast for this to have a chance to be in the release.

@tupui tupui added defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize maintenance Items related to regular maintenance tasks enhancement A new feature or improvement labels May 10, 2021
Copy link
Member

@tupui tupui left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a huge update, thanks again for putting this together.

I did a first pass and I would first advise to clean these in order to facilitate the review. Things like reverting American spellings would lower the diff size a lot.
For workers try to follow the same usage as in differential_evolution. There is also a helper function to handler the pool (_util.MapWrapper).
Also, there are lot of TODOs. I would try to clean these as well (put them in the issue with a checklist instead).

scipy/optimize/_shgo.py Outdated Show resolved Hide resolved
Comment on lines 27 to 28
minimizer_kwargs=None, options=None, sampling_method='simplicial',
workers=None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

New arguments should be keywords only.

Suggested change
minimizer_kwargs=None, options=None, sampling_method='simplicial',
workers=None):
minimizer_kwargs=None, options=None, sampling_method='simplicial', *,
workers=None):

Comment on lines 182 to 183
local minimisation routine every iteration. If False then only the
final minimiser pool will be run. Defaults to False.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

American spelling with z should be used as the whole documentation uses optimize. Here and in other places.

Suggested change
local minimisation routine every iteration. If False then only the
final minimiser pool will be run. Defaults to False.
local minimization routine every iteration. If False then only the
final minimizer pool will be run. Defaults to False.

Comment on lines 214 to 216
workers : int optional
Uses `multiprocessing.Pool <multiprocessing>`) to sample and run the
local serial minimizatons in parrallel.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add default and information about -1 (copy from differential_evolution).

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=None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should default to 1 here as differential_evolution and other ones.

self.symmetry, self.bounds, self.g_cons,
self.g_args)
if self.disp:
print('Constructing and refining simplicial complex graph '
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Print or logging? Same bellow.

logging.info('=' * 30)

if self.HC.V[x] not in self.minimizer_pool:
self.minimizer_pool.append(self.HC.V[x])

if self.disp:
logging.info('Neighbors:')
logging.info('Neighbours:')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are not needed either (while in this case it does not really matter, it increases the diff which makes it harder to review). In general we use American spelling.

results['x'] = self.xl_maps[ind_sorted[0]] # Save global minima
results['fun'] = self.f_maps[ind_sorted[0]] # Save global fun value

self.xl_maps = np.ndarray.tolist(self.xl_maps)
self.f_maps = np.ndarray.tolist(self.f_maps)
return results

# TODO: In scipy version delete this
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also for caching there is functools.lru_cache

import pytest
from pytest import raises as assert_raises, warns
from nose.tools import nottest
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't use nose. If you want to skip a test, use @pytest.mark.skip(reason="This is not run because ...")

self.process_fpool = self.proc_fpool_g
else:
self.workers = workers
self.pool = mp.Pool(processes=workers) #TODO: Move this pool to
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe rely on _util.MapWrapper instead

@Stefan-Endres
Copy link
Contributor Author

Thank you very much @tupui . Apologies about the doctrings I think I accidentally pulled the wrong direction in the last merge. I will update and clean this as soon as possible to get it up to review standard.

@tupui
Copy link
Member

tupui commented May 10, 2021

Thank you very much @tupui . Apologies about the doctrings I think I accidentally pulled the wrong direction in the last merge. I will update and clean this as soon as possible to get it up to review standard.

No problem at all 😃

@tupui
Copy link
Member

tupui commented May 14, 2021

BTW, we recently added bounds to Nelder-Mead, so the solver_args needs to be updated. In general, maybe we should have this list generally accessible some place else. Otherwise we will always forget to update it here.

@tupui
Copy link
Member

tupui commented Jul 20, 2021

Hi @Stefan-Endres just checking in 😃 Any update?

@Stefan-Endres
Copy link
Contributor Author

Hi, apologies for letting the time slip on this, I will have the weekend free to work on an hopefully finish this.

@tupui
Copy link
Member

tupui commented Jul 22, 2021

Hi, apologies for letting the time slip on this, I will have the weekend free to work on an hopefully finish this.

No problem 😃 Thanks!

Copy link
Member

@tupui tupui left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did another superficial pass.

There are unchecked item in your list. Are you still working on these or is the PR ready for more review?

Did you check that the resource usage is constant? (CPU, RAM). With more workers in the loop and the rewrite, things can go bad fast.

@@ -241,7 +260,7 @@ def shgo(func, bounds, args=(), constraints=None, n=None, iters=1,
specified using the `bounds` argument.

While most of the theoretical advantages of SHGO are only proven for when
``f(x)`` is a Lipschitz smooth function, the algorithm is also proven to
``f(x)`` is a Lipschitz smooth function. the algorithm is also proven to
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
``f(x)`` is a Lipschitz smooth function. the algorithm is also proven to
``f(x)`` is a Lipschitz smooth function, the algorithm is also proven to

sfield=self.func, sfield_args=self.args,
symmetry=self.symmetry,
constraints=self.constraints,
# constraints=self.g_cons,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# constraints=self.g_cons,

self.fn = 0 # Number of feasible sampling points evaluations performed
self.hgr = 0 # Homology group rank

# Default settings if no sampling criteria.
if (self.n is None) and (self.iters is None):
self.n = 128
self.nc = 0 # self.n
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
self.nc = 0 # self.n
self.nc = 0

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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
self.nc = 0 # n # Sampling points to sample in current iteration
self.nc = 0 # Sampling points to sample in current iteration

Comment on lines +529 to +530
if (not isinstance(constraints, tuple)) and (not
isinstance(constraints, list)):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (not isinstance(constraints, tuple)) and (not
isinstance(constraints, list)):
if (not isinstance(constraints, tuple)) and \
(not isinstance(constraints, list)):

# Generate sampling points.
# Generate uniform sample points in [0, 1]^m \subset R^m
if self.n_sampled == 0:
self.C = self.sobol_points(n, dim)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where is sobol_points defined?

import itertools
import json
import decimal
from functools import lru_cache # For Python 3 only
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
from functools import lru_cache # For Python 3 only
from functools import lru_cache

HC.V: The cache of vertices and their connection
HC.H: Storage structure of all vertex groups

:param dim: int, Spatial dimensionality of the complex R^dim
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The doc is in the wrong format here, should be numpydoc.


# Printing
if printout:
print("=" * 19)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some print left

"""
Cache objects
"""
class VertexCacheBase(object):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not needed anymore

Suggested change
class VertexCacheBase(object):
class VertexCacheBase:

@Stefan-Endres
Copy link
Contributor Author

Hi @tupui apologies I should've commented that I need a bit more time before another review pass. Aside from the TODOs there is also a few more things to implement from your earlier comments in this PR.

I compared CPU times on the GO benchmarks, but not RAM usage yet I didn't think about that. I suppose the approriate measure would be to benchmark memory usage on the for workers=1 vs workers >=2 on the scipy GO testsuite, I will run these and post the results here.

@tupui
Copy link
Member

tupui commented Jul 28, 2021

Hi @tupui apologies I should've commented that I need a bit more time before another review pass. Aside from the TODOs there is also a few more things to implement from your earlier comments in this PR.

No worries, let me know when you're ready or need help with something 😃

I compared CPU times on the GO benchmarks, but not RAM usage yet I didn't think about that. I suppose the approriate measure would be to benchmark memory usage on the for workers=1 vs workers >=2 on the scipy GO testsuite, I will run these and post the results here.

Thanks! Would be great so that we have full confidence with all the changes.

@dschmitz89
Copy link
Contributor

@Stefan-Endres : I was thinking of adding support for the Bounds class for shgo. Is it worthwhile at the moment before this PR is merged?

@Stefan-Endres
Copy link
Contributor Author

@Stefan-Endres : I was thinking of adding support for the Bounds class for shgo. Is it worthwhile at the moment before this PR is merged?

I think it's worth doing now as I would need to rebase this PR anyway, I am not certain when I will have full free weekend again to finish this.

@tupui
Copy link
Member

tupui commented Jan 24, 2022

FYI we've merged the improvement @Stefan-Endres.

@wilrop
Copy link

wilrop commented Mar 10, 2022

Hi, I'm just checking in to see if this update is already in the latest version of scipy because I'm still having problems of this kind. If it is already included, I can open a new issue so that maybe someone can look at it.

@tupui
Copy link
Member

tupui commented Mar 10, 2022

Hi, I'm just checking in to see if this update is already in the latest version of scipy because I'm still having problems of this kind. If it is already included, I can open a new issue so that maybe someone can look at it.

Hi @wilrop, thanks for checking in. Please don't open any new issue on this. This PR is still not merged and would address the issues. @Stefan-Endres is still working on it.

@tupui tupui added this to the 1.9.0 milestone Mar 14, 2022
@tupui
Copy link
Member

tupui commented Mar 14, 2022

@Stefan-Endres I am adding a release milestone. We have until around June to finish this.

@tupui
Copy link
Member

tupui commented May 19, 2022

@Stefan-Endres FYI the branching for 1.9 is next week (26 of May). Let me know if I can help. Otherwise this will be pushed to 1.10.

@Stefan-Endres
Copy link
Contributor Author

Stefan-Endres commented May 20, 2022

@tupui I will work on finishing this tonight. I think the only thing we really need is to test for memory leaks? The remaining issue #14533 could be done in a second PR.

@tupui
Copy link
Member

tupui commented May 20, 2022

Great thank you! I will be on the lookout to review. The memory leak is just (I hope) a sanity check, so yes just running it twice with a low count and with a high count of worker is fine.

@Stefan-Endres
Copy link
Contributor Author

I ran into compilation rebasing this PR last night, currently running the test suite.

Some preliminary results indicate memory usage of multiprocessing is only about 0.025% bigger (expected due to bigger class initiated to handle the cache processing).

Running mem_tests.py:

import sys

from memory_profiler import profile
from pympler import asizeof
from scipy.optimize._shgo import SHGO
from scipy.optimize import rosen

# Define SHGO class
bounds = [(0,2), (0, 2), (0, 2), (0, 2)]
shc1 = SHGO(rosen, bounds, iters=4)  # With only 1 worker
shc2 = SHGO(rosen, bounds, iters=4, workers=2)  # With 2 workers using multiprocessing
shc1.iterate_all(), shc2.iterate_all()
print('Results with 1 cpu:')
print(f'asizeof.asizeof(shc1) (SHGO class) = {asizeof.asizeof(shc1)}')
print(f'asizeof.asizeof(shc1.HC) (Cache container) = {asizeof.asizeof(shc2.HC)}')
print('-')
print('Results with 2 cpus:')
print(f'asizeof.asizeof(shc2) (SHGO class) = {asizeof.asizeof(shc2)}')
print(f'asizeof.asizeof(shc2.HC) (Cache container) = {asizeof.asizeof(shc2.HC)}')
print('-')
print('Difference in memory usage:')
print(f'SHGO class: asizeof.asizeof(shc1) - asizeof.asizeof(shc2)  = {asizeof.asizeof(shc1) - asizeof.asizeof(shc2)}')
print(f'Cache containers: asizeof.asizeof(shc1.HC) - asizeof.asizeof(shc2.HC) = {asizeof.asizeof(shc1.HC) - asizeof.asizeof(shc2.HC)}')

produces

$ python mem_tests.py
Results with 1 cpu:
asizeof.asizeof(shc1) (SHGO class) = 29499688
asizeof.asizeof(shc1.HC) (Cache container) = 68495504
-
Results with 2 cpus:
asizeof.asizeof(shc2) (SHGO class) = 29517392
asizeof.asizeof(shc2.HC) (Cache container) = 68513208
-
Difference in memory usage:
SHGO class: asizeof.asizeof(shc1) - asizeof.asizeof(shc2)  = -17704
Cache containers: asizeof.asizeof(shc1.HC) - asizeof.asizeof(shc2.HC) = -17704

However, using the shgo library it is only about 0.003% bigger:

from shgo._shgo import SHGO  # Use shgo instead of scipy version
from scipy.optimize import rosen

Produces

$ python -i mem_tests.py 
Results with 1 cpu:
asizeof.asizeof(shc1) (SHGO class) = 29514720
asizeof.asizeof(shc1.HC) (Cache container) = 68510536
-
Results with 2 cpus:
asizeof.asizeof(shc2) (SHGO class) = 29516480
asizeof.asizeof(shc2.HC) (Cache container) = 68512264
-
Difference in memory usage:
SHGO class: asizeof.asizeof(shc1) - asizeof.asizeof(shc2)  = -1760
Cache containers: asizeof.asizeof(shc1.HC) - asizeof.asizeof(shc2.HC) = -1728

I don't fully understand why the scipy code which is idenitical has a bigger overhead for multiprocessing.

Next step is to finish rebasing and rerun the same tests on all the scipy GO benchmarks, please let me know if anyone prefers using different methods of measurement for testing the memory usage.

tylerjereddy and others added 6 commits May 22, 2022 22:44
…ns that occur for large parameter values (#14036)

Avoid the numerical issues in `roots_jacobi`, `roots_gegenbauer`, and
`_gen_roots_and_weights` that occur for large values of a and b in
`roots_jacobi(n,a,b)`. Also helps with large values of n (i.e. n > 200).

[ci skip]
@Stefan-Endres
Copy link
Contributor Author

Stefan-Endres commented May 22, 2022

Apologies no idea why this happend, @tupui perhaps it would be better to close this PR and start a new one? Despite rebasing I still get this error when trying to build and run the tests:

RuntimeError: Running scipy/stats/_generate_pyx.py failed
Cythonizing sources
Traceback (most recent call last):
  File "/home/stefan_endres/projects/scipy_fucked/setup.py", line 534, in <module>
    setup_package()
  File "/home/stefan_endres/projects/scipy_fucked/setup.py", line 518, in setup_package
    generate_cython()
  File "/home/stefan_endres/projects/scipy_fucked/setup.py", line 267, in generate_cython
    raise RuntimeError("Running cythonize failed!")
RuntimeError: Running cythonize failed!

@mdhaber mdhaber reopened this May 22, 2022
@mdhaber mdhaber changed the base branch from main to maintenance/0.5.2.x May 22, 2022 22:14
@mdhaber mdhaber changed the base branch from maintenance/0.5.2.x to main May 22, 2022 22:14
@mdhaber mdhaber closed this May 22, 2022
@tupui
Copy link
Member

tupui commented May 23, 2022

Seems like you rebased agains a maintenance branch and not main.

@Stefan-Endres feel free to open a new PR if that's easier for you. FYI, since you opened this PR the infrastructure evolved a bit and we have a new convenient developer interface. If you have issues, I would suggest you looking into creating a new conda environment with the meson yaml and using python do.py.

@Stefan-Endres
Copy link
Contributor Author

@Stefan-Endres feel free to open a new PR if that's easier for you.

I'm not sure how to undo the rebase and GitHub automatically adding a bunch of reviewer requests etc. I think it will be easier if you don't mind.

FYI, since you opened this PR the infrastructure evolved a bit and we have a new convenient developer interface. If you have issues, I would suggest you looking into creating a new conda environment with the meson yaml and using python do.py.

Thank you, I will try this after work tonight instead of the old build instructions.

@tupui
Copy link
Member

tupui commented May 23, 2022

Let me know if you need help. We can make a call as well 😃

@wilrop
Copy link

wilrop commented May 23, 2022

Hi, I'm probably super late to the party but I think there is a (different?) bug in how constraints are interpreted. If I add a single dictionary as a constraint everything is fine, but if I give a list of dictionaries as constraints it completely ignores all of them. Is this something that is handled with the current improvements or completely new? I've developed a minimal example, so I can create a new issue if that helps.

@tupui
Copy link
Member

tupui commented May 23, 2022

@wilrop if this bug is different from existing ones you are free to create an issue. Another possibility is that you wait a bit for an update here so we can try to reproduce your bug with latest changes (better IMO).

@wilrop
Copy link

wilrop commented May 23, 2022

Okay no problem, I'll wait and open an issue if the problem persists.

@Stefan-Endres
Copy link
Contributor Author

@wilrop please open an issue with your code, there are unittests that have tuples of active constraints so this is possibly a new bug.

@wilrop
Copy link

wilrop commented May 24, 2022

I've just discovered that the bug is not specific to the list of constraints but rather with indexing in the x variable for constraints. I'll open a new issue anyway and tag @Stefan-Endres and @tupui as you have already been involved here.

@Stefan-Endres Stefan-Endres mentioned this pull request May 29, 2022
5 tasks
o-alexandre-felipe added a commit to o-alexandre-felipe/scipy that referenced this pull request Jun 30, 2022
Solution suggested by jjramsey in
scipy#14589 (comment)

The answer at the time was that an important refactor was in progress
scipy#14013

This pull request is now closed and the bug persists in version 1.8.1 accordingly to this stack overflow question
https://stackoverflow.com/q/72794609/12750353
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected enhancement A new feature or improvement maintenance Items related to regular maintenance tasks scipy.optimize
Projects
None yet