vallis/libstempo

Some units support bugs

Opened this issue · 1 comments

When loading a tempopulsar with units=True, some functionality does not work well or is wrong.

  1. psr.designmatrix() fails when a pulsar has units=True, with the error message pasted below. This seems to stem from the fixunits=True bit of code. With fixunits=False the designmatrix is correctly returned.
  2. The units are not correct for all tempo2 parameters. For instance, ELONG and ELAT are actually in radians in Tempo2, yet libstempo assumes they are in degrees.
---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:982](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:982), in Quantity.to_value(self, unit, equivalencies)
    981 try:
--> 982     scale = self.unit._to(unit)
    983 except Exception:
    984     # Short-cut failed; try default (maybe equivalencies help).

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1159](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1159), in UnitBase._to(self, other)
   1157         return self_decomposed.scale [/](https://file+.vscode-resource.vscode-cdn.net/) other_decomposed.scale
-> 1159 raise UnitConversionError(f"'{self!r}' is not a scaled version of '{other!r}'")

UnitConversionError: 'Unit("deg")' is not a scaled version of 'Unit(dimensionless)'

During handling of the above exception, another exception occurred:

UnitConversionError                       Traceback (most recent call last)
File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:1350](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:1350), in Quantity.__float__(self)
   1349 try:
-> 1350     return float(self.to_value(dimensionless_unscaled))
   1351 except (UnitsError, TypeError):

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:985](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:985), in Quantity.to_value(self, unit, equivalencies)
    983 except Exception:
    984     # Short-cut failed; try default (maybe equivalencies help).
--> 985     value = self._to_value(unit, equivalencies)
    986 else:

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:891](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:891), in Quantity._to_value(self, unit, equivalencies)
    889 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    890     # Standard path, let unit to do work.
--> 891     return self.unit.to(
    892         unit, self.view(np.ndarray), equivalencies=equivalencies
    893     )
    895 else:
    896     # The .to() method of a simple unit cannot convert a structured
    897     # dtype, so we work around it, by recursing.
    898     # TODO: deprecate this?
    899     # Convert simple to Structured on initialization?

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1195](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1195), in UnitBase.to(self, other, value, equivalencies)
   1194 else:
-> 1195     return self._get_converter(Unit(other), equivalencies)(value)

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1124](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1124), in UnitBase._get_converter(self, other, equivalencies)
   1122                 return lambda v: b(converter(v))
-> 1124 raise exc

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1107](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1107), in UnitBase._get_converter(self, other, equivalencies)
   1106 try:
-> 1107     return self._apply_equivalencies(
   1108         self, other, self._normalize_equivalencies(equivalencies)
   1109     )
   1110 except UnitsError as exc:
   1111     # Last hope: maybe other knows how to do it?
   1112     # We assume the equivalencies have the unit itself as first item.
   1113     # TODO: maybe better for other to have a `_back_converter` method?

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1085](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/core.py:1085), in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
   1083 other_str = get_err_str(other)
-> 1085 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")

UnitConversionError: 'deg' (angle) and '' (dimensionless) are not convertible

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
[/Users/vhaasteren/research/projects/combine-ptas/combine-pulsars.ipynb](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/combine-pulsars.ipynb) Cell 20 line 1
----> [1](vscode-notebook-cell:/Users/vhaasteren/research/projects/combine-ptas/combine-pulsars.ipynb#Y133sZmlsZQ%3D%3D?line=0) psr_epta.designmatrix()

File libstempo/libstempo.pyx:1754, in libstempo.libstempo.tempopulsar.designmatrix()

File [~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:1352](https://file+.vscode-resource.vscode-cdn.net/Users/vhaasteren/research/projects/combine-ptas/~/miniconda3/envs/pint/lib/python3.10/site-packages/astropy/units/quantity.py:1352), in Quantity.__float__(self)
   1350     return float(self.to_value(dimensionless_unscaled))
   1351 except (UnitsError, TypeError):
-> 1352     raise TypeError(
   1353         "only dimensionless scalar quantities can be "
   1354         "converted to Python scalars"
   1355     )

TypeError: only dimensionless scalar quantities can be converted to Python scalars

Note that units support is becoming more important for IPTA projects also, where for some projects the design matrices of PINT and Tempo2 need to match. Currently, the default PINT units are not the same as the Tempo2 units for all parameters. Enterprise cannot compare units without libstempo support