shbhuk / barycorrpy

Python version of Barycorr
GNU General Public License v3.0
37 stars 6 forks source link

barycorrpy unit tests failed under new installation #31

Closed RussLaher closed 4 years ago

RussLaher commented 4 years ago

I installed python3 via Anaconda3-2019.10-Linux-x86_64.sh on TACC machine Lonestar5:

Operating System: SUSE Linux Enterprise Server 12 SP3 CPE OS Name: cpe:/o:suse:sles:12:sp3 Kernel: Linux 4.4.103-6.38-default Architecture: x86-64

which includes astropy-4.0, astroquery-0.3.10, barycorrpy-0.2.2.1, jplephem-2.12, and astroscrappy-1.0.8.

cd ~/barycorrpy/ export PATH=/home1/07040/rrlaher7/anaconda3/bin:$PATH

$ python -m unittest discover |& grep -i error ERROR: test_JDUTC_BJDTDB (barycorrpy.tests.test_barycorrpy.Barycorrpy_tests) AttributeError: 'numpy.ndarray' object has no attribute 'to_value' astropy.units.core.UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("d")' astropy.units.core.UnitConversionError: '' (dimensionless) and 'd' (time) are not convertible ERROR: test_flux_weighting (barycorrpy.tests.test_barycorrpy.Barycorrpy_tests) AttributeError: 'numpy.ndarray' object has no attribute 'to_value' astropy.units.core.UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("d")' astropy.units.core.UnitConversionError: '' (dimensionless) and 'd' (time) are not convertible ERROR: test_hip_id (barycorrpy.tests.test_barycorrpy.Barycorrpy_tests) AttributeError: 'numpy.ndarray' object has no attribute 'to_value' astropy.units.core.UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("d")' astropy.units.core.UnitConversionError: '' (dimensionless) and 'd' (time) are not convertible ERROR: test_hip_id_astropy_obs (barycorrpy.tests.test_barycorrpy.Barycorrpy_tests) AttributeError: 'numpy.ndarray' object has no attribute 'to_value' astropy.units.core.UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("d")' astropy.units.core.UnitConversionError: '' (dimensionless) and 'd' (time) are not convertible ERROR: test_inputs (barycorrpy.tests.test_barycorrpy.Barycorrpy_tests) AttributeError: 'numpy.ndarray' object has no attribute 'to_value' astropy.units.core.UnitConversionError: 'Unit(dimensionless)' is not a scaled version of 'Unit("d")' astropy.units.core.UnitConversionError: '' (dimensionless) and 'd' (time) are not convertible FAILED (errors=5)

Any suggestions on how to proceed?

RussLaher commented 4 years ago

Here is more version info:

Successfully installed astropy-4.0 Successfully installed barycorrpy-0.2.2.1

$ python --version Python 3.7.4

$ which python /home1/07040/rrlaher7/anaconda3/bin/python

RussLaher commented 4 years ago

I downgraded numpy on the TACC machine to 1.16.2 and reran the barycorrpy tests successfully.

shbhuk commented 4 years ago

Hi @RussLaher, Thanks for bringing this up. I will have to check what part of the Astropy Quantities routine is failing under the new Numpy version. I will attempt to fix this in the new version release - v0.3

shbhuk commented 4 years ago

Hi, I'm reopening this issue so that I can remember to look into this numpy update issue and figure out what is breaking..

njcuk9999 commented 4 years ago

Also getting this, maybe the traceback will help solve this faster (we don't really want to have to downgrade packages if possible)

command ran:

bkwargs = {'ra': 263.3741740568687, 'dec': -5.745214357217233, 'epoch': 2457205.9999942128,
 'px': 21.353919213862373, 'pmra': -41.93841615099209, 'pmdec': -95.48427622624564, 'lat': 19.825252, 'longi': -155.468876, 'alt': 4204.0, 'rv': 0.0, 'leap_update': False}

barycorrpy.get_BC_vel(2458849.04505247, zmeas=0.0, **bkwargs)

versions:

traceback:

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/barycorrpy/barycorrpy.py in get_BC_vel(JDUTC, starname, hip_id, ra, dec, epoch, pmra, pmdec, px, rv, obsname, lat, longi, alt, zmeas, ephemeris, leap_dir, leap_update)
    157                  zmeas=zm,
    158                  loc=loc,
--> 159                  ephemeris=ephemeris, leap_dir=leap_dir, leap_update=leap_update,**star_output)
    160         vel.append(a[0])
    161         warning.append(a[1])

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/barycorrpy/barycorrpy.py in BCPy(JDUTC, ra, dec, epoch, pmra, pmdec, px, rv, zmeas, loc, ephemeris, leap_dir, leap_update)
    189     ##### NUTATION, PRECESSION, ETC. #####
    190 
--> 191     r_pint, v_pint = PINT.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
    192 
    193     r_eci = r_pint[0]  # [m]

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/barycorrpy/PINT_erfautils.py in gcrs_posvel_from_itrf(loc, toas, tts)
     86 
     87     # Get dX and dY from IERS B in arcsec and convert to radians
---> 88     dX = np.interp(mjds, iers_tab['MJD'], iers_tab['dX_2000A']) * asec2rad
     89     dY = np.interp(mjds, iers_tab['MJD'], iers_tab['dY_2000A']) * asec2rad
     90 

<__array_function__ internals> in interp(*args, **kwargs)

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity.py in __array_function__(self, function, types, args, kwargs)
   1496             function_helper = FUNCTION_HELPERS[function]
   1497             try:
-> 1498                 args, kwargs, unit, out = function_helper(*args, **kwargs)
   1499             except NotImplementedError:
   1500                 return self._not_implemented_or_raise(function, types)

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity_helper/function_helpers.py in interp(x, xp, fp, *args, **kwargs)
    812     from astropy.units import Quantity
    813 
--> 814     (x, xp), _ = _quantities2arrays(x, xp)
    815     if isinstance(fp, Quantity):
    816         unit = fp.unit

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity_helper/function_helpers.py in _quantities2arrays(unit_from_first, *args)
    348     # as we want to allow arbitrary unit for 0, inf, and nan.
    349     try:
--> 350         arrays = tuple((q._to_own_unit(arg)) for arg in args)
    351     except TypeError:
    352         raise NotImplementedError

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity_helper/function_helpers.py in <genexpr>(.0)
    348     # as we want to allow arbitrary unit for 0, inf, and nan.
    349     try:
--> 350         arrays = tuple((q._to_own_unit(arg)) for arg in args)
    351     except TypeError:
    352         raise NotImplementedError

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity.py in _to_own_unit(self, value, check_precision)
   1337             try:
   1338                 as_quantity = Quantity(value)
-> 1339                 _value = as_quantity.to_value(self.unit)
   1340             except TypeError:
   1341                 # Could not make a Quantity.  Maybe masked printing?

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity.py in to_value(self, unit, equivalencies)
    732             except Exception:
    733                 # Short-cut failed; try default (maybe equivalencies help).
--> 734                 value = self._to_value(unit, equivalencies)
    735             else:
    736                 value = self.view(np.ndarray)

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/quantity.py in _to_value(self, unit, equivalencies)
    663             equivalencies = self._equivalencies
    664         return self.unit.to(unit, self.view(np.ndarray),
--> 665                             equivalencies=equivalencies)
    666 
    667     def to(self, unit, equivalencies=[]):

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/core.py in to(self, other, value, equivalencies)
    985             return UNITY
    986         else:
--> 987             return self._get_converter(other, equivalencies=equivalencies)(value)
    988 
    989     def in_units(self, other, value=1.0, equivalencies=[]):

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/core.py in _get_converter(self, other, equivalencies)
    916                             pass
    917 
--> 918             raise exc
    919 
    920     def _to(self, other):

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/core.py in _get_converter(self, other, equivalencies)
    902         try:
    903             return self._apply_equivalencies(
--> 904                 self, other, self._normalize_equivalencies(equivalencies))
    905         except UnitsError as exc:
    906             # Last hope: maybe other knows how to do it?

/scratch/bin/anaconda3/envs/aperoenv/lib/python3.7/site-packages/astropy/units/core.py in _apply_equivalencies(self, unit, other, equivalencies)
    886         raise UnitConversionError(
    887             "{} and {} are not convertible".format(
--> 888                 unit_str, other_str))
    889 
    890     def _get_converter(self, other, equivalencies=[]):

UnitConversionError: '' (dimensionless) and 'd' (time) are not convertible

Going into the offending function it seems that it is the np.interp function:

np.interp(mjds, iers_tab['MJD'], iers_tab['dX_2000A']) * asec2rad

where for me:

In [59]: mjds                                                                                                                                                
Out[59]: array([58848.54505247])

In [60]: iers_tab['MJD']                                                                                                                                     
Out[60]: <Quantity [37665., 37666., 37667., ..., 58748., 58749., 58750.] d>

In [61]: iers_tab['dX_2000A']                                                                                                                                
Out[61]: <Quantity [0.0e+00, 0.0e+00, 0.0e+00, ..., 2.0e-06, 1.1e-05, 1.9e-05] arcsec>

So your iers_tab['MJD'] has units "d" and mjds has no units...hence the UnitConversionError I'm not sure how this worked in the past... unless now astropy is pushing units to iers_tab['MJD'] when it wasn't before or np.interp was stripping units in the past... (I would be more likely to believe astropy changed... but it downgrading numpy worked it would suggest the latter).

you guess you have two solutions (1) add units to mjds or (2) remove units from iers_tab['MJD']

njcuk9999 commented 4 years ago

i.e. for me this worked in my local barycorrpy files: In PINT_erfautils.py:

    dX = np.interp(mjds, iers_tab['MJD'].value, iers_tab['dX_2000A'].value) * asec2rad
    dY = np.interp(mjds, iers_tab['MJD'].value, iers_tab['dY_2000A'].value) * asec2rad

and

    xp = np.interp(mjds, iers_tab['MJD'].value, iers_tab['PM_x'].value) * asec2rad
    yp = np.interp(mjds, iers_tab['MJD'].value, iers_tab['PM_y'].value) * asec2rad

but again I'm not sure how you make older versions compatible -- that depends why it worked previously

shbhuk commented 4 years ago

Hi, I am not sure what is going on with the numpy updated versions, it is breaking some other packages of mine too.. I'll need to check

shbhuk commented 4 years ago

Okay, so I ran two different virtual envs, where I experimented with the astropy and numpy versions. Case I: Astropy 4.0, Numpy 1.16.4 = Worked Case II: Astropy 4.0, Numpy 1.18.1 = Fails Case III: Astropy 3.2.3, Numpy 1.16.4 = Worked Case IV: Astropy 3.2.3, Numpy 1.18.1 = Worked

So its combination of the latest Numpy and Astropy updates that are causing this. Very interesting..

astrowright commented 4 years ago

I might post a bug report to the astropy folks—I suspect they're aware of the problem and can recommend a solution.

On Sun, Jan 19, 2020 at 4:31 PM Shubham notifications@github.com wrote:

Okay, so I ran two different virtual envs, where I experimented with the astropy and numpy versions. Case I: Astropy 4.0, Numpy 1.16.4 = Worked Case II: Astropy 4.0, Numpy 1.18.1 = Fails Case III: Astropy 3.2.3, Numpy 1.16.4 = Worked Case IV: Astropy 3.2.3, Numpy 1.18.1 = Worked

So its combination of the latest Numpy and Astropy updates that are causing this. Very interesting..

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/shbhuk/barycorrpy/issues/31?email_source=notifications&email_token=AASHFHJDDFSXAXDVW47YGELQ6TBEFA5CNFSM4KEIRGVKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJK5EAA#issuecomment-576049664, or unsubscribe https://github.com/notifications/unsubscribe-auth/AASHFHJPVYCBZGS2OMPJMCTQ6TBEFANCNFSM4KEIRGVA .

--


Jason T Wright Professor of Astronomy & Astrophysics Deputy Director, Center for Exoplanets and Habitable Worlds NEID Project Scientist http://neid.psu.edu/ he/him/his https://sites.psu.edu/astrowright/ https://sites.psu.edu/astrowright/ @Astro_Wright

http://twitter.com/Astro_Wright

shbhuk commented 4 years ago

So, @njcuk9999 iers_tab['mjd'] has had units of days in both Astropy 3.2.3 and 4.0, so that's not the issue. The numpy documentation for np.interp does not seem to have changed.

What baffles me is that the combination of Astropy 3.2.3, Numpy 1.18.1, seems to work. I have checked the inputs with this sample script

from astropy.utils.iers import conf,IERS_A, IERS_A_URL, IERS_B, IERS_B_URL, IERS
from astropy.utils.data import download_file
from astropy.time import Time
import astropy.table as table
import numpy as np
import astropy

# Load IERS table
iers_b_file = download_file(IERS_B_URL, cache=True)
iers_b = IERS_B.open(iers_b_file)
IERS.iers_table = iers_b
iers_tab = IERS.iers_table

mjd_astropy = iers_tab['MJD']

JDUTC = Time(2458000, format='jd', scale='utc')
JDTT = Time(2458000, format='jd', scale='tt')

toas = JDUTC
if toas.isscalar:
    ttoas = Time([toas])
N = len(ttoas)
mjds = np.asarray(ttoas.mjd)
PM = iers_tab['PM_x'].value

# Check version and interpolation
print(astropy.__version__,':', np.__version__)
print(mjds, mjd_astropy, PM)
print(np.interp(mjds, mjd_astropy.value, PM))
print(np.interp(mjds, mjd_astropy, PM))

I'm not sure what exactly is causing this emergent error, however for now I will include an update with v0.3 where I shall modify the code to use just the values from the Astropy quantities for the interpolation in PINT_erfautils

shbhuk commented 4 years ago

As Jason @astrowright suggested, I think the issue might be this update under Astropy -

Astropy 4.0, Units: : Point 3 image

https://github.com/astropy/astropy/pull/8808