SHGO algorithm

[reneeotten]

Description

SciPy v1.2 introduced two new global optimization algorithms and these should be supported in lmfit as well (see #527). This PR adds the shgo (simplicial homology global optimization) method.

Tests were added to check make sure the implementation in lmfit gives the same result as calling scipy.optimize.shgo directly, using the example in thes SciPy function docstring. It seems to me there might be an upstream issue when using bounds at 0 or None, but I'll need to narrow that down a bit more to make sure that's indeed the case. The tests in test_covariance_matrix.py and a few others I tried myself did not work when a boundary was exactly set to zero (throwing an error from scipy/optimize/_shgo_lib/triangulation.py). Replacing the 0 with a very small value (I picked machine precision for a float64) seems to solve that problem. That required a bit of work before calling prepare_fit, but I think the logic works there.

Similarly, I think something looks weird when no bounds are specified (shgo will internally set the boundaries then to + and - 1e50), and it looks like the initial starting point for the minimization is then zero (at least when specifying finite bounds, it will take the the value in the middle between min and max), which throws the same error as when having bounds at zero. At least I cannot get it to pass the test_numdifftools_no_bounds in test_covariance_matrix.py, so for now I did not add the shgo method there.

I think this has nothing to do with the implementation in lmfit, but I would welcome a few more pairs of eyes to make sure... So, please take a look at the PR, try it on a few of your favorite minimization problems to see how it behaves and let me know.

Type of Changes

• Bug fix
• New feature: add SHGO algorithm
• Refactoring / maintenance: a few small fixes after the last PR merge
• Documentation / examples

Tested on

Python: 2.7.16, 3.6.8, and 3.7.2
lmfit: 0.9.12+77.ga05de1a, scipy: 1.2.1, numpy: 1.16.2, asteval: 0.9.13, uncertainties: 3.0.3, six: 1.12.0

Verification

Have you

• included docstrings that follow PEP 257?
• referenced existing Issue and/or provided relevant link to mailing list?
• verified that existing tests pass locally?
• squashed/minimized your commits and written descriptive commit messages?
• added or updated existing tests to cover the changes?
• updated the documentation?
• added an example?

[reneeotten]

fixed the logic to skip tests when scipy <= 1.2 ; last failing test is unrelated to this PR and is caused by minimal differences in the uncertaintes (only for PY35, so probably related to scipy 1.1). We can probably change the number of decimals for that test to three and than it will not be as flaky anymore..

[newville]

@reneeotten The behavior of SHGO for bounds or initial values near zero is pretty strange, and should probably be raised as an issue on scipy.

Other than that, I think this looks great. I like the decorator for specifying the algorithms for testing, and thanks for the better solution for suppressing the flake8 message about unused imports.

I think it's OK to leave the "tweak 0 to epsilon" code in, at least until it gets solved upstream, but maybe add a brief note about the issue...

Anyone have concerns or objections to merging this in the next couple of days?

[reneeotten]

@newville let's keep it open for a little bit; I'll take a look at the actual paper first to see if I am missing something obvious. I would like to get a small, self-contained example and submit this to the author his GitHub (or scipy) to see what they think. Something that I also don't fully understand yet is that while the solutions for the eggholder-test are the same when directly using scipy or through lmfit, the number of function evaluations in some steps are not. It's my understanding that the method should be deterministic, which means that this cannot be the case.

Either way, I would prefer people to try out this PR for a bit and see how it behaves before merging. Of course, it's hard to tell if that ever happens or if it will only get used when merged into master... What about giving me (and others) till the end of the week or so to do a bit more testing?

[newville]

@reneeotten OK, that's fine with me!

[andyfaff]

If you have a MWE that displays issues in the shgo implementation, then scipy is always interested in fixing them.

…

[reneeotten]

If you have a MWE that displays issues in the shgo implementation, then scipy is always interested in fixing them.

I know @andyfaff. But before saying that there are issues in the implementation I want to make sure it's not me who introduced these issue when adding the algorithm to lmfit ;) I will investigate a bit more this week before opening an issue with the author and/or scipy.

[newville]

@reneeotten is this merge-able? I'd like to start working on preparing 0.9.13 soon. Should this be included?

[reneeotten]

@newville no, I am sorry but it's not ready to merge yet. I did look into it a bit, but haven't finished it yet... so feel free to prepare a release without this.

I do think though that the weirdness I saw when bounds at zero is from lmfit; it appears that during the shgo algorithm it minimizes in smaller pieces in the parameter space within the bounds. So that means that in some cases it will minimize in a small area, where for example in the VoigtModel, the "best" value for sigma in that case is exactly zero. This will throw an ZeroDivisionError from asteval for obvious reasons in calculating the line shape. So we probably will need to add the tiny parameter in this (and some other line shape functions) as well to avoid this.

[newville]

@reneeotten OK, that's fine.

For the weird behaviour seen with bounds at 0, I'm not sure I've see that. But I am increasingly wary of auto-selecting initial values (including the lineshape models built into lmfit). If Max and Min are set but not an initial value, assuming that Initial Value is (Max-Min)/2 may cause serious problems. For sure, handling an initial value of 0 is a challenge. Like, did that mean ~1.e-3 or ~1.e-150?

[reneeotten]

okay, I think this should be good to go now. It turned out that the issue with "bounds at zero" was actually "our" fault and not in the shgo algorithm. It's reproducible also with brute and caused by the fact that for the lineshape function the minimum for sigma is set to 0.0. Methods that will sample points between bounds (i.e., brute, shgo) assume that this is a valid point, which is actually not the case as it will throw (correctly) a ZeroDivisionError.

I have updated the built-in models so that the default bounds are set to "machine precision" instead of zero, and also included such check in the lineshape functions. Please double-check that I didn't miss anything or messed something up... Of course it will still fail is a user explicitly sets the minimum bound for sigma to 0.0, but there isn't much we can do about that.

[codecov]

Codecov Report

Merging #545 into master will increase coverage by 1.58%.
The diff coverage is 95.94%.

@@            Coverage Diff             @@
##           master     #545      +/-   ##
==========================================
+ Coverage   83.28%   84.86%   +1.58%     
==========================================
  Files          11       11              
  Lines        3170     3198      +28     
==========================================
+ Hits         2640     2714      +74     
+ Misses        530      484      -46

Impacted Files	Coverage Δ
lmfit/lineshapes.py	87.3% <100%> (+15.54%)	⬆️
lmfit/models.py	87.2% <100%> (+5.23%)	⬆️
lmfit/printfuncs.py	86.77% <50%> (-0.11%)	⬇️
lmfit/minimizer.py	90.22% <94.28%> (+0.17%)	⬆️
lmfit/parameter.py	81.84% <0%> (+0.26%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 7d9bce8...d1524d8. Read the comment docs.

[newville]

@reneeotten Oh, that makes sense, and it's great that we can add SHGO. But a general comment, and maybe a suggestion for maybe thinking about how this could go differently.

For division-by-zero, we may want to be more careful than just saying that any sigma parameter must be > eps (~2.e-16). For one thing, eps represents the smallest detectable relative change between two floats, but not the smallest value that can be represented. The smallest non-zero more like 1e-308 (1/np.finfo(np.float).max).

In my experience, working with tiny values (say, below 1.e-12) is definitely prone to several "close to zero" errors. But, floating point is supposed to handle that. Like, 1..e150 / 1.e-150 does not actually raise ZeroDivisionError (it is close to giving Inf, which causes other problems). I think the point is that some very small values ought to be an allowed, but that we probably set an upper limit to guard against division-by-zero. Whether that is 1.e-13, 2.e-16, 1.e-50, or 1.e-308 is probably less important than applying it consistently.

But it also seems to me that the problem isn't with sigma, but with divide-by-zero. So, it might be better to guard against division-by-zero explicitly where the division happens, replacing .... / sigma with ... / max(TINY, sigma) (or something similar) where TINY could be 2.e-16 or 1.e-50 or 1/np.finfo(np.float).max.

I would guess that if lineshapes.py were fixed to explicitly guard against division-by-zero, then the rest of the code could retain min=0 for sigma and other parameters, which would probably be more readable and preserve the convenient (if not exactly correct) notion that one does not need to think about "how small of a zero do you mean?".

[reneeotten]

For division-by-zero, we may want to be more careful than just saying that any sigma parameter must be > eps (~2.e-16). For one thing, eps represents the smallest detectable relative change between two floats, but not the smallest value that can be represented. The smallest non-zero more like 1e-308 (1/np.finfo(np.float).max).

In my experience, working with tiny values (say, below 1.e-12) is definitely prone to several "close to zero" errors. But, floating point is supposed to handle that. Like, 1..e150 / 1.e-150 does not actually raise ZeroDivisionError (it is close to giving Inf, which causes other problems). I think the point is that some very small values ought to be an allowed, but that we probably set an upper limit to guard against division-by-zero. Whether that is 1.e-13, 2.e-16, 1.e-50, or 1.e-308 is probably less important than applying it consistently.

Sure, I see your point and am open to suggestions for a more informed decision about what a "small value to replace zero" should be. I can look a bit more and see if, for example, numpy or scipy guard for ZeroDivisionErrors in their code and how they do it.

But it also seems to me that the problem isn't with sigma, but with divide-by-zero. So, it might be better to guard against division-by-zero explicitly where the division happens, replacing .... / sigma with ... / max(TINY, sigma) (or something similar) where TINY could be 2.e-16 or 1.e-50 or 1/np.finfo(np.float).max.

Yes, the problem is as you said with divide-by-zero. It was my intention to change all these occurences in lineshapes.py and models.py, but might have missed something. It just happens to be mainly the variable sigma in our lineshape functions. I decided to do this replacement at the top of each lineshape function instead of in the actual calculation of the lineshape for readability - that can, of course, be changed.

I would guess that if lineshapes.py were fixed to explicitly guard against division-by-zero, then the rest of the code could retain min=0 for sigma and other parameters, which would probably be more readable and preserve the convenient (if not exactly correct) notion that one does not need to think about "how small of a zero do you mean?".

That's what I initially tried, but unfortunately that alone is not sufficient. The error occurs in the __residual function where it calls params.update_constraints().

/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/minimizer.py:598: in penalty
    r = self.__residual(fvars, apply_bounds_transformation)
/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/minimizer.py:526: in __residual
    params.update_constraints()
/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/parameter.py:216: in update_constraints
    _update_param(name)
/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/parameter.py:212: in _update_param
    self._asteval.symtable[name] = par.value
/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/parameter.py:788: in value
    return self._getval()
/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/parameter.py:770: in _getval
    check_ast_errors(self._expr_eval)
/Users/renee/Library/Python/3.7/lib/python/site-packages/lmfit-0.9.13+5.g9ff6068-py3.7.egg/lmfit/parameter.py:24: in check_ast_errors
    expr_eval.raise_exception(None)

The only way to circumvent this problem (as far as I can tell) is to actually set minimum bounds to something "slightly larger than zero" in the actual model class. Again, I am open to other suggestions.

[newville]

@reneeotten FWIW, I don't really have a very strong preference for what "TINY" to avoid divide-by-zero.

I would guess that if lineshapes.py were fixed to explicitly guard against division-by-zero, then the rest of the code could retain min=0 for sigma and other parameters, which would probably be more readable and preserve the convenient (if not exactly correct) notion that one does not need to think about "how small of a zero do you mean?".

That's what I initially tried, but unfortunately that alone is not sufficient. The error occurs in the __residual function where it calls params.update_constraints().

Oh, that seems kind of weird and seems like that might not be limited to one particular solver. Do you have a simple example for this?

[reneeotten]

Oh, that seems kind of weird and seems like that might not be limited to one particular solver. Do you have a simple example for this?

Yes, as I said earlier the same issue happens with brute. I think though that everything works as intended... at least I would expect the error if you try to evaluate gamma=0.0 in the example below (just put it in the tests directory to run). The issue will only show up for methods that do not actually perform a minimization: brute just evaluates the function at specified parameter values and the error in shgo occurs when setting up the sampling points, where a similar evaluation happens.

import os

import numpy as np

from lmfit.models import VoigtModel


# load data to be fitted, assuming this file is in the "tests" directory
data = np.loadtxt(os.path.join(os.path.dirname(__file__), '..', 'examples',
                               'test_peak.dat'))
x = data[:, 0]
y = data[:, 1]

# define the model and initialize parameters
mod = VoigtModel()
params = mod.guess(y, x=x)
params['amplitude'].set(min=25, max=70)
params['sigma'].set(min=0, max=1)
params['center'].set(min=5, max=15)

result = mod.fit(y, params, x=x, method='brute')

[newville]

@reneeotten Thanks. That suggests to me that when sigma=0 the expression for the height constraint fails. So, maybe change that definition to something like

amplitude/max(1.e-99, sigma*sqrt(2*pi)))*wofz((1j*gamma)/max(1.e-99, sigma*sqrt(2))).real

or whatever "smallest positive number" should be..... 2e-16 would be fine too.

[reneeotten]

okay, I have made changes to the lineshapes and calculation of heigth and fwhm in built-in models so that they'll result in finite output values and not throw a ZeroDivisionError anymore. I hope that's what your intention was as well @newville? Please do double-check and make sure that I din't make a mistake somewhere....

[newville]

@reneeotten yes, that's great! That should solve the problem for all solvers or anyone evaluating these models.

Ready to merge?

[reneeotten]

Ready to merge?
Yes, I think it's good to go!

[newville]

@reneeotten OK, merging -- this is really great!!