Make Postpic aware of Units in a proper way

[Ablinne]

From my research there are a lot of packages available that allow unit-aware calculations in python. Most introduce a new class Quantity or similar, which represents a quantity with a value and a unit.
They are also able to just handle bare units without values.
Most are able to also wrap a numpy ndarray inside a Quantity.
The most stable and mature packages seem to be

• pint
• astropy.units
• sympy.physics.units

current versions of pint also support wrapping numpy addays.

Pint
Pint is lightweight and an arch linux package is available in the community repo. It wraps numpy arrays by subclassing them. I see three options of how we could use pint:

• Just use the units by themselves in the .unit member of Field and Axis
• Use Quantity wrapped arrays in .grid_node and .matrix and get rid of the .unit member
• combine this with Field should return copies instead of in-place editing. #75 and instead of subclassing np.ndarray, we subclass pint's Quantity

astropy.units
I didn't test ist, is not packaged for arch.

sympy.physics.units
Large package, takes 400ms to import

[Ablinne]

Turns out pint does not subclass arrays. Quantitys contain arrays, just like the Field contains an array. Is it meaningful to put a container in a container? or should these things be somehow next to each other on the same level? Would subclassing be easier? How should we do this according to the numpy API? Numpys manual does talk about subclassing...

[Ablinne]

• Before numpy 1.11 it was impossible to modify the input to ufuncs before they are called. This means that, if unit conversion needs to be done beforehand, the ufunc will actually be called twice. No proper way around it.
• Starting with numpy 1.11 there is a new interface, __numpy_ufunc__/__array_ufunc__ (renamed in 1.13) which makes this easy and possible. This is meant for both, subclasses of ndarray and containers that contain an ndarray. The old interface was meant for subclasses, but accidentally works also for containers.
• Pint does not yet use this new interface yet and does calculations twice.
• In a container class it is necessary to override __add__ and so on, in a subclass not. IMHO this is a big plus for a subclass. Using nested containers means duplicated effort.
• We could have an ndarray subclass with units called e.g. ndquantity (possibly using the units part of pint) and another subclass Field of this that adds the axes and stuff. But this would still have the problem of doing computations twice
• We could add this functionality only for numpy 1.11 and newer, while leaving the old behaviour intact for older versions.

[Ablinne]

• So far ufuncs are not supported at all. Units are represended as strings and not really handled at all
• One would need to make sure that in numpy 1.11 and newer the .unit attribute is never really used, while it must still be around for older numpys. .unit should be a property that either returns the value of ._unit_compat, which is what we currently use, or yields the proper unit of the underlying object. A setter could be used to trigger in-place unit conversion.
• two possibilities to implement the subset of the ufuncs which are also operators: implement __add__ and use that implementation from __array_ufunc__ or the other way round: implement the logic within __array_ufunc__ and implement __add__ in terms of np.add().
  • second option comes for free with a generic implementation of __array_ufunc__, but then we will also have to implement ufuncs for older numpy which will have the problem of calculating stuff twice for older numpy and at that point we must just live with that.
  • first option splits responsibilities: some ufuncs implemented in one places, other ufuncs implemented in other places but will mean least implementation effort, because we only implement stuff for the new numpy

[Ablinne]

• a mix of both options might look like this: implement the logic in a _operate method that is used to apply operators like +, -, *, /, >, <, ==, ... and contains the logic. The __add__ and so on just use these (old and new numpy), as does __array_ufunc__. Then the implementation of everything is just in 2 places, namely in _operate for operators and in __array_ufunc__ for the remaining ufuncs.

[Ablinne]

• instead of an attribute _unit_compat we could just have a dummy class ndquantity_compat which is just an array with an additional .unit-string member, similar to the example InfoArray from the numpy docs. This will be used in older numpy, but the full-fledged implementation is used in newer numpy. Then Field.unit would be inherited from ndquantity_compat.unit or ndquantity.unit or, if Field is still a container, we just need to return Field.matrix.unit.

[skuschel]

One more option to assess: There is a package called python-quantities. It is not packed in arch, but its pure python, so pip can easily install it. The Quantity is a np.ndarray subclass. They get around the double computation by overriding pretty much all operator overloading functions -- but what is the point of subclassing then? In addition if you do

import quantities as pq
a = 2*pq.meter
b = 2*pq.ft

then a+b works, but np.add(a,b) breaks because its not the same unit. However, that seems to be a problem we can live with. Have a look, the quantities code seems to be clearly written and well organized:
https://github.com/python-quantities/python-quantities/blob/master/quantities/quantity.py

[Ablinne]

That seems nice! And using __array_ufunc__ it would be easy to just repair the ufuncs for compatible units for python > 1.11 in python-quantities! They work around the double computation by rejecting not matching units in array_prepare, __array_ufunc__ on the other hand could just perform the necessary conversions (by just deferring the work to __add__ and sorts where appropriate).

[skuschel]

https://nrc-cnrc.github.io/MetroloPy/ is linked on the Scipy page https://www.scipy.org/topical-software.html . Maybe this is another option to assess.