search

SNEWS2 / snewpy

A Python package for working with supernova neutrinos
https://snewpy.readthedocs.io
BSD 3-Clause "New" or "Revised" License
24 stars 17 forks source link

Units simplification bug in the RateCalculator #267

Closed Sheshuk closed 9 months ago

Sheshuk commented 9 months ago

Bug in astropy.units https://github.com/astropy/astropy/issues/15200 sometimes leads to units not simplified after the rate calculation, resulting in units like cm^2/cm^2 or h/s.

Sheshuk commented 9 months ago

Looks like the problem comes from the presence of various systems in astropy.units. Centimeters are both defined in SI and CGS, and they are imported in astropy.units.__init__.py: https://github.com/astropy/astropy/blob/b74412ce6359650a6ba3a8800ab865b3f7e31a4a/astropy/units/__init__.py#L28-L34 Since cgs is imported last, astropy.units.cm is always astropy.units.cgs.cm.

We have four variants to define a centimeter:

  1. u_cgs = u.cgs.cm
  2. u_si = u.si.cm
  3. u_def = u.cm is always u.cgs.cm
  4. u_str = u.Unit("cm") - this is a result of parsing the string and a lookup in the units registry. And this is where result is unstable from run to run.

The registry is a dictionary and can be accessed by

registry = u.get_current_unit_registry().registry

This registry is defined by the function u.set_enabled_units - by default it loads many unit systems. We can limit this lookup registry, for example, setting u.set_enabled_units(u.si), but this doesn't help the fact that units provided by user will come from the main astropy.units namespace and can be from different unit systems

Sheshuk commented 9 months ago

Probably the solution for our case is just to call Quantity.decompose function after we do the calculation, so the result is in SI units