Only strip the numbers from a unit if it's dimensionless

[henrikjacobsenfys]

(I've complained about this at standups, but it's bothering me too much now and I have to get it fixed)

General

scipp and its underlying unit handling library stores units weirdly. They mostly get stored how you set them, e.g. if you make a variable with the unit meV, it'll keep that unit. But powers of units are another matter. They get converted to SI units, but so that .value is unchanged. So 117 1/meV^2 get stored as 117 4.558e+45 1/J^2. .value is correct, the unit is also correct, except it looks funny.

It's silly when the units are compatible. For example, 1 cm/1 m can be stored as 1 0.01 . In our Parameters, we therefore strip the numbers from the units.

This means that working with .value can silently alter the results by 40 or more orders of magnitude. This is not OK.

I have had to build in workarounds for this in EasyDynamics, storing an extra scipp variable to track the units properly,, and it'll get worse when I start working with normalized data, since the natural unit of the y axis is either 1/meV (for S(Q,omega)) or worse, mbarn/sr/f.u./meV (f.u. is formula unit) (for cross sections). mbarn gets stored as 1e-30m^2, and so absolute units will be wrong by a factor $10^{30}$. It's a constant nuisance in EasySpectroscopy, and it WILL trip us up in other libraries as well. It will also mess things up for users, because it happens silently and unpredictably. Making aliases of units only fixes units that we think of, and will not cover all use cases.

Idea

We should stop removing the number from the unit, so that we stop converting parameter values silently. The only acceptable time to do it is when the unit is only a number (like the 1cm/1m case). It should be easy enough to change the code to check if there's anything but units in the unit field.

[damskii9992]

Actually, I have a proposal to your proposal. The only issue I have with the conversion of units, is how they are displayed weirdly. So my counter-proposal is this:

We don't convert the unit to its base ever, leaving it as-is by scipp internally.
all methods which display or return the unit, converts only the returned unit to the correct format.

This way, it will look like our parameter is in units of "1/meV^2", while it is actually stored as that weird "1/J^2" internally.

[henrikjacobsenfys]

I'm not sure how to determine the right format? Should we store the unit the user assigned and print that? (re)use desired_unit?

[damskii9992]

Well, can just store the string that the user supplied and convert to that whenever we return the unit. This won't work for dependent parameters, obviously, as they have no user-supplied unit, but there we can simply use the desired unit, if it has been supplied.

[damskii9992]

Also, I had @jl-wynen report the non-working alias to the underlying units package:
llnl/units#466
Because it i really annoying that the unit aliases don't work properly . . .

[jl-wynen]

I'd like to clarify how units are encoded because that might help you understand why it behaves in surprising ways sometimes. Sorry for the spam if you already know all this.

Memory Representation

Units are always stored as powers of SI units plus a multiplier plus flags. Let's ignore flags for now. We can access the underlying representation using to_dict:

>>> import scipp as sc
>>> sc.Unit("m").to_dict()
{'__version__': 2, 'multiplier': 1.0, 'powers': {'m': 1}}
>>> sc.Unit("2m").to_dict()
{'__version__': 2, 'multiplier': 2.0, 'powers': {'m': 1}}
>>> sc.Unit("kg*m/s^").to_dict()
{'__version__': 2, 'multiplier': 1.0, 'powers': {'m': 1, 'kg': 1, 's': -1}}

(The version is Scipp-specific and we use it for serialisation to HDF5)

This also applies to non-SI units:

>>> sc.Unit("Å").to_dict()
{'__version__': 2, 'multiplier': 1e-10, 'powers': {'m': 1}}
>>> sc.Unit("meV").to_dict()
{'__version__': 2, 'multiplier': 1.6021766339999998e-22, 'powers': {'m': 2, 'kg': 1, 's': -2}}

The underlying LLNL/Units library is designed to be very efficient and only use a couple words to encode a unit. So it doesn't use a hash map like the output above but a compact C struct. The multiplier if a double.

Limitations

Each power only gets a few bits in the C struct, see https://units.readthedocs.io/en/latest/details/unit_base.html but note that we configure the library to use a larger 64bit base type which roughly doubles the bit size of each power.

This impacts which units can be encoded. In the configuration we use, the allowed ranges of exponents are quite large. E.g., we support meter to have powers [-256, 255]. But, e.g., mol is limited to [-8, 7].

I doubt you will ever run into this, but if you do, now you know why. :-)

Flags

Each unit carries a couple of flags. In Scipp, we mostly disallow them and raise an error if parsing a unit string would set flags. With some exceptions, e.g.,

>>> sc.Unit("degC").to_dict()
{'__version__': 2, 'multiplier': 1.0, 'e_flag': True, 'powers': {'K': 1}}
>>> repr(sc.Unit("degC"))
'Unit(K, e_flag=True)
>>> str(sc.Unit("degC"))
'°C'

This unit has multiplier=1 because it has the same scale as kelvin but it has an offset. The library will handle this automatically in operations. You don't super need to know this, but I mentioned it here for completeness.

String formatting

When you render a unit to a string, the underlying library matches the multiplier, powers, and flags to a long list of known units. It will try to find the best match in the list to determine how to render the unit. This works well for many units but it doesn't cover all possible cases. When there is no good match, the library falls back to an SI unit. That is how you get

>>> str(sc.Unit("1/meV^2"))
'3.89564355265760538e+43/J^2'

But keep in mind that this formatting works off of the stored multiplier, powers, and flags. The original string when you construct the unit is irrelevant:

>>> str(sc.Unit("Å"))
'Å'
>>> str(sc.Unit("10^-10 m"))
'Å'
>>> str(sc.Unit("10^10 Å"))
'm'

Unit aliases

Scipp's unit aliases wrap user defined units. When you define, e.g.,

sc.units.aliases[name] = name

you ultimately call

addUserDefinedUnit(name, unit_from_string(name));

So you create a new user defined unit. The Unit Aliases documentation shows the impact of this.

So aliases are inserted into the formatting table and can be matched when the library wants to render a unit. This can be the best match for some units but there is no guarantee.

[henrikjacobsenfys]

Thanks!