Replacing inaccurate Voigt1D algorithm with Humlicek approximation or scipy.special.wofz

[dhomeier]

Description

This pull request is to address inaccuracies and errors in the McLean et al. algorithm used for Voigt1D as reported in https://dx.doi.org/10.1016/j.jqsrt.2018.03.019

Hopes to fix #7256

This PR replaces the default algorithm with the complex error function computed from one of the 16-coefficient algorithms suggested by Franz Schreier in #7256 (comment) and available at https://atmos.eoc.dlr.de/tools/lbl4IR/ .
As an alternative implementation I have added the compiled Faddeeva function from scipy.special.wofz, which is still more accurate at almost the same speed, but of course depends on scipy – algorithms can be selected with the new option method.
In addition I have added some clarification on the normalisation of the function discussed in #7256 (comment) and on Slack #general.

I could still need some advice on correctly implementing the new option. Originally I was hoping to simply set the method in __init__ so that the Voigt1D instance would call the appropriate function as cls._wofz. But settings in __init__ are apparently not passed through to the classmethods evaluate and fit_deriv, so I needed to set the whole stuff up in __new__. This still looks messy and unnecessary convoluted, so I'd be grateful for suggestions for a cleaner solution.
Moreover I was expecting this would at least create class instances with the given method set, e.g. for

voi_w = models.Voigt1D(amplitude_L=2.0/np.pi, fwhm_L=1.0, fwhm_G=0.9, method='wofz')
voi_h = models.Voigt1D(amplitude_L=2.0/np.pi, fwhm_L=1.0, fwhm_G=0.9, method='humlicek2')

voi_w(x) would always compute profiles with the scipy function and voi_h(x) with the Humlicek algorithm; however after the second instantiation the method for voi_w is set to humlicek2 as well!
Otherwise it seems to work, but I am afraid the behaviour above could lead to confusion and unintended results.

Update: instantiation now works as desired; the respective algorithms can be queried with voi_h.method, voi_w.method.

[pllim]

Thanks, @dhomeier ! I am tentatively milestoning this for next release with no backport, given looks like a behavorial change that is not backward compatible (at a glance). Can you please elaborate the background for this PR, if the conversation has happened on Slack etc, for future reference?

[dhomeier]

Sorry for pushing the Open button prematurely when I was trying to pull down the Draft option!
The topic had been discussed in the linked issue for some time, although the normalisation (integral is neither 1 nor sqrt(pi)) only came up in the recent Slack thread.
In principle it constitutes a bug fix, since the function call will still work the same without any additional arguments (except producing less inaccurate results).

[pllim]

Thanks for the clarification! I updated the milestone.

[dhomeier]

Some numbers on the performance:

In [5]: lam = np.linspace(-10, 10, 200001)

In [6]: for f_l in 10**np.arange(-6, 6.1, 0.5):
   ...:     print(f_l)
   ...:     v = Voigt1D(x_0=5, amplitude_L=2/(np.pi*f_l), fwhm_L=f_l, fwhm_G=1, method='humlicek2')
   ...:     %timeit v(lam + 1.e-9 * np.random.rand())
   ...:     v = Voigt1D(x_0=5, amplitude_L=2/(np.pi*f_l), fwhm_L=f_l, fwhm_G=1, method='faddeeva')
   ...:     %timeit v(lam + 1.e-9 * np.random.rand())
   ...: 
1e-06
Voigt1D: using algorithm hum2zpf16c
25.5 ms ± 481 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
28.8 ms ± 137 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
3.162277660168379e-06
Voigt1D: using algorithm hum2zpf16c
25.2 ms ± 133 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
28.8 ms ± 127 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
1e-05
Voigt1D: using algorithm hum2zpf16c
25.2 ms ± 229 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
29.2 ms ± 467 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
3.1622776601683795e-05
Voigt1D: using algorithm hum2zpf16c
25.8 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
28.8 ms ± 119 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.0001
Voigt1D: using algorithm hum2zpf16c
25.4 ms ± 534 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
29.3 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.00031622776601683794
Voigt1D: using algorithm hum2zpf16c
25.2 ms ± 129 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
28.9 ms ± 102 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.001
Voigt1D: using algorithm hum2zpf16c
25.3 ms ± 165 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
29 ms ± 562 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.0031622776601683794
Voigt1D: using algorithm hum2zpf16c
25.4 ms ± 348 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
28.9 ms ± 229 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.01
Voigt1D: using algorithm hum2zpf16c
25.5 ms ± 713 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
28.8 ms ± 87.2 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.03162277660168379
Voigt1D: using algorithm hum2zpf16c
32.3 ms ± 8.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
33 ms ± 5.65 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.1
Voigt1D: using algorithm hum2zpf16c
24.8 ms ± 370 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
31.6 ms ± 3.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
0.31622776601683794
Voigt1D: using algorithm hum2zpf16c
45.4 ms ± 16.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
37.9 ms ± 14.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
1.0
Voigt1D: using algorithm hum2zpf16c
40.5 ms ± 22.2 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
35.8 ms ± 7.33 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
3.1622776601683795
Voigt1D: using algorithm hum2zpf16c
23 ms ± 1.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Voigt1D: using algorithm wofz
26.3 ms ± 84.1 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
10.0
Voigt1D: using algorithm hum2zpf16c
16.4 ms ± 141 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
18 ms ± 130 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
31.622776601683793
Voigt1D: using algorithm hum2zpf16c
13.3 ms ± 547 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
14.4 ms ± 46.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
100.0
Voigt1D: using algorithm hum2zpf16c
13.4 ms ± 699 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
16.3 ms ± 1.44 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
316.22776601683796
Voigt1D: using algorithm hum2zpf16c
13 ms ± 139 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
12.3 ms ± 49.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
1000.0
Voigt1D: using algorithm hum2zpf16c
13.4 ms ± 552 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
11.5 ms ± 102 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
3162.2776601683795
Voigt1D: using algorithm hum2zpf16c
12.8 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
11.3 ms ± 105 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
10000.0
Voigt1D: using algorithm hum2zpf16c
12.4 ms ± 211 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
9.44 ms ± 44 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
31622.776601683792
Voigt1D: using algorithm hum2zpf16c
24.2 ms ± 4.82 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
17.2 ms ± 431 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
100000.0
Voigt1D: using algorithm hum2zpf16c
16.2 ms ± 862 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
15.5 ms ± 1.05 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)
316227.7660168379
Voigt1D: using algorithm hum2zpf16c
14 ms ± 514 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
14.1 ms ± 333 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
1000000.0
Voigt1D: using algorithm hum2zpf16c
14.4 ms ± 427 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Voigt1D: using algorithm wofz
14.3 ms ± 631 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So for small fwhm_G/fwhm_L ratios Humlicek is slightly faster, not really a big gain, but the accuracy is certainly good enough for all practical purposes (see Fig. 3 of https://dx.doi.org/10.1093/mnras/sty1680), and it is a pure Numpy implementation.

[keflavich]

lgtm. I agree it would be better to avoid the duplication in __new__ and __init__, but I don't know how to do that.

I've left some minor UI-related comments.

I haven't tested this. It might be nice to add a Voigt1D fitting example to the docs somewhere to showcase this?

[saimn]

You cannot use a class attribute to store the method, which should be specific to the instance. Otherwise you could get weird results:

>>> voi_w = models.Voigt1D(amplitude_L=2.0/np.pi, fwhm_L=1.0, fwhm_G=doppler, method='wofz')
>>> voi_w.__class__._wofz is wofz
True
>>> voi_h = models.Voigt1D(amplitude_L=2.0/np.pi, fwhm_L=1.0, fwhm_G=doppler, method='humlicek2')
>>> voi_w.__class__._wofz is wofz
False

Using voi_w after instantiating voi_h will not give you what you expect.

[saimn]

evaluate and fit_deriv being class methods, the only solution would be to have two different classes.

[dhomeier]

The first docs build warning

/home/docs/checkouts/readthedocs.org/user_builds/astropy/envs/11177/lib/python3.8/site-packages/astropy/modeling/fitting.py:1166: AstropyUserWarning: The fit may be unsuccessful; check fit_info['message'] for more information.
warnings.warn("The fit may be unsuccessful; check "

seems to show up in other builds as well; not sure what to make of

/home/docs/checkouts/readthedocs.org/user_builds/astropy/envs/11177/lib/python3.8/site-packages/astropy/modeling/functional_models.py:docstring of astropy.modeling.functional_models.Voigt1D.hum2zpf16c:4: WARNING: Undefined substitution referenced: "x".

is it complaining that the docstring mentions |x| and y which do not actually appear in the function (being z.real and z.imag, really)?

[dhomeier]

You cannot use a class attribute to store the method, which should be specific to the instance. Otherwise you could get weird results:

>>> voi_w = models.Voigt1D(amplitude_L=2.0/np.pi, fwhm_L=1.0, fwhm_G=doppler, method='wofz')
>>> voi_w.__class__._wofz is wofz
True
>>> voi_h = models.Voigt1D(amplitude_L=2.0/np.pi, fwhm_L=1.0, fwhm_G=doppler, method='humlicek2')
>>> voi_w.__class__._wofz is wofz
False

Using voi_w after instantiating voi_h will not give you what you expect.

Yes, exactly what I was worrying about. It was not quite clear why the functions needed to be classmethods; but not sure if it's worth creating separate classes, if that is the only way to solve this.

[dhomeier]

I would have been content to just handle method as an additional argument, but as it's not listed with the Parameters I could not even get it to be accepted as kwarg without the __new__ and __init__ voodoo...

[perrygreenfield]

@dhomeier I think the fact that evaluate was previously a class method is just a historical accident and isn't required to be (I haven't recently looked at the code, but I can't think of a reason it needs to be). There are other examples of models that define evaluate also as static methods or instance methods. Perhaps it is as simple as just trying it as an instance method. If that works then setting which algorithm becomes much simpler.

[dhomeier]

Thanks for the clarification @perrygreenfield – using instance methods and setting the implementation method on __init__ makes this much easier!

[dhomeier]

@plim my understanding was that the changelog entry should go into section 4.3 here, since this will go directly into the master branch, and would be moved accordingly when preparing the backport to 4.0. Or should I move it into 4.0.5 here?

The previous commit had 3 unrelated test failures in coordinates and other sub-packages, everything else passed.

I'll look into an example for the docs as suggested by @keflavich ; functionally this should be ready now.

[pllim]

@dhomeier , change log should be in the section that is defined by milestone. Thanks!

[perrygreenfield]

This is a great improvement! Just a couple of minor comments or questions.

[perrygreenfield]

Somehow I missed that the last requested change had been made. Sorry for the late review.

[nden]

@dhomeier Could you move the change log to the same section as the milestone.
Also there are some failures in the doc build which need to be fixed before merging.

[dhomeier]

@nden I have apparently sufficiently clarified the definition of z to get rid of that warning.
On the remaining

/home/docs/checkouts/readthedocs.org/user_builds/astropy/envs/11177/lib/python3.8/site-packages/scipy/optimize/optimize.py:282: RuntimeWarning: Values in x were outside bounds during a minimize step, clipping to bounds
  warnings.warn("Values in x were outside bounds during a "
/home/docs/checkouts/readthedocs.org/user_builds/astropy/envs/11177/lib/python3.8/site-packages/astropy/modeling/fitting.py:1166: AstropyUserWarning: The fit may be unsuccessful; check fit_info['message'] for more information.
  warnings.warn("The fit may be unsuccessful; check "

The out-of-bounds RuntimeWarning appears when running the example from compound-models.rst with 4.2 as well, and seems to have been newly introduced in Scipy 1.16.
The second warning is raised by the "Tied Constraints" example in example-fitting-contraints.rst even with Scipy 1.5.4, because the fit does not complete in 100 iterations. No idea why it has not been triggered before, but both seem clearly outside the scope of this PR.

[dhomeier]

FTR
fitted_model = fitter(hbeta_combo, wave, flux, maxiter=200)
does converge in 111 iterations for me.

[dhomeier]

@kevlavich I have thought about your suggestion of adding a new example for Voigt1D to the docs, but haven't come up with a good showcase highlighting the specific properties of the Voigt function so far.
I think the docs in general could use some clarification of the various meanings of amplitude (peak or whatever else) in the pre-defined models, and their relation to the normalisation properties, especially since they are labelled under "Profiles". When treated as a profile (i.e. for convolution), one should really use a form integrating to unity.
But that goes a bit beyond the case of the Voigt function, so should go into a separate PR.

[dhomeier]

Or should I rebase in case there is a fix in master that I have missed somehow?

[nden]

Thanks @dhomeier!

[dhomeier]

Thanks for troubleshooting, @nden!