Improve handling of spectral units

[ajjackson]

It looks like spectral intensities are sometimes handled in energy units (such as when imported from CASTEP DOS). If these are converted to/from their proper (reciprocal energy) units, this can cause divide-by-zero errors as pint performs a wavenumber <-> wavelength conversion and takes the reciprocal of the spectral intensities, including empty bins.

Here we introduce another mechanism for dealing with spectral units: a few extra unit conversions in a special "reciprocal_spectroscopy" context (which comes an "rs" alias for convenience). This stops Pint from choking on the 1/(1/cm) <-> 1/meV conversions.

The broadening system is made a bit more robust; previously the variable-width broadening assumed that things can always be converted to energy units.

Perhaps the most serious issue addressed here is the way the Spectrum data "getters" (hidden behind a @property interface) work.

The current spectrum getters create a unit conversion factor and then
multiply the raw array by it. There are two problems with this:

• Multiplying any iterator by a Quantity or Unit leads to Pint
  checking through all the data for consistency. This can be quite
  expensive, and is unnecessary for these trusted array objects. When
  iterating over spectra yielded from a Spectrum1DCollection, the cost
  can become quite noticeable.

• Some unit conversions in the spectroscopy context involve taking
  a reciprocal, e.g. conversion from cm -> 1/cm is allowed. But if
  this is done to a factor the data itself will not be inverted. If
  we start with a Spectrum1D with x_data in wavenumber (1/cm) then

  spectrum.x_data.to("cm")

  and

  spectrum.x_data_unit = "cm"
  spectrum.x_data

  will not give the same result; the second method leaves the actual
  values untouched (as cm / (1/cm) = 1).

To further improve these conversions I included some technically-redundant energy-wavenumber conversions in the "reciprocal_spectroscopy" context. With these available, .to(..., "rs") conversions can avoid needlessly taking reciprocals of the data and raising divide-by-zero errors. That should be more efficient, too, but it's not clear if that works in practice.

[github-actions]

Test Results

    22 files  ±0      22 suites  ±0   29m 15s ⏱️ - 1m 15s
 1 048 tests ±0   1 042 ✅ ±0   6 💤 ±0  0 ❌ ±0 
10 420 runs  ±0  10 357 ✅ ±0  63 💤 ±0  0 ❌ ±0 

Results for commit 517e8d0. ± Comparison against base commit 75fe42e.

♻️ This comment has been updated with latest results.

[ajjackson]

I performed a little speed test of conversions in different contexts. With this definition file:

@context reciprocal_spectroscopy = rs
    1 / [frequency] <-> 1 / [length]: 1 / value / speed_of_light
    1 / [energy] -> 1 / [frequency]: planck_constant * value
    1 / [frequency] -> 1 / [energy]: value / planck_constant
    [length] -> 1 / [energy]: value / planck_constant / speed_of_light
    1 / [energy] -> [length]: value * planck_constant * speed_of_light
    1 / [length] -> [energy]: value * planck_constant * speed_of_light
    [energy] -> 1 / [length]: value / planck_constant / speed_of_light
@end

@context bracketed_spectroscopy = bs
    1 / [frequency] <-> 1 / [length]: 1 / (value * speed_of_light)
    1 / [energy] -> 1 / [frequency]: planck_constant * value
    1 / [frequency] -> 1 / [energy]: value / planck_constant
    [length] -> 1 / [energy]: value / (planck_constant * speed_of_light)
    1 / [energy] -> [length]: value * (planck_constant * speed_of_light)
    1 / [length] -> [energy]: value * (planck_constant * speed_of_light)
    [energy] -> 1 / [length]: value / (planck_constant * speed_of_light)
@end

and this speed test code:

import timeit

import numpy as np
from pint import UnitRegistry


ureg = UnitRegistry()
ureg.load_definitions("./definitions.txt")

sparse_data = np.random.random((10000, 200))
sparse_data[sparse_data > 0.5] = 0

sparse_spectrum = ureg.Quantity(sparse_data, "1/cm")

print("Default spectroscopy context:")
print(timeit.timeit("sparse_spectrum.to('meV', 'spectroscopy')", globals=globals(), number=1000))

print("Reciprocal spectroscopy context:")
print(timeit.timeit("sparse_spectrum.to('meV', 'reciprocal_spectroscopy')", globals=globals(), number=1000))

print("Reciprocal spectroscopy context with parenthesis groups:")
print(timeit.timeit("sparse_spectrum.to('meV', 'bracketed_spectroscopy')", globals=globals(), number=1000))

the output is

Default spectroscopy context:
/Users/adam.jackson/micromamba/envs/euphonic-np-test/lib/python3.10/site-packages/pint/facets/plain/quantity.py:1054: RuntimeWarning: divide by zero encountered in divide
  return self.__class__(other_magnitude / self._magnitude, 1 / self._units)
4.869501542067155
Reciprocal spectroscopy context:
3.7165777920745313
Reciprocal spectroscopy context with parenthesis groups:
2.6495610419660807

So we see that

• .to() is quite expensive
• there is almost a performance factor x2 gained by preparing appropriate definitions

Usually a factor x2 is not a big thing to scream about, but this is so slow that it is a credible bottleneck and worth sorting out.

That being said, things narrow a lot when replacing .to() with .ito() to convert in-place (which has a 100x impact in this example!), so we should look for ways to benefit from that inside Euphonic/Abins:

Default spectroscopy context:
/Users/adam.jackson/micromamba/envs/euphonic-np-test/lib/python3.10/site-packages/pint/facets/plain/quantity.py:1054: RuntimeWarning: divide by zero encountered in divide
  return self.__class__(other_magnitude / self._magnitude, 1 / self._units)
0.04182825004681945
Reciprocal spectroscopy context:
0.03280979208648205
Reciprocal spectroscopy context with parenthesis groups:
0.03438079194165766

[ajjackson]

I think it could be worth refactoring how spectra work; currently they store a numpy array and unit information, and dynamically construct a Quantity when properties are accessed. This is quite expensive, and it would be simpler to store a full Quantity. The current behaviour of "lazily" dealing with unit changes could be retained by deferring an ito() to the property getter.

Presumably some features work directly on the private array, and these would have to deal with an extra call to magnitude. It seems unlikely that this performance benefit outweight the downsides, but we can do some benchmarking if unsure.

Anyway, this can all be a separate PR!

[ajjackson]

@mducle @oerc0122 This could use another set of eyes before merging. In principle it is straightforward stuff, but messing around with unit conversions always scares me a bit. Are there any additional sanity checks or stress tests we can try in order to make sure we didn't accidentally give permission to do anything wrong?

[oerc0122]

All looks good if it works, just a couple of stylistic questions.