Repository navigation
Transforming a CartesianRepresentation gives the wrong answer for a non-contiguous float matrix #15503
Description
Description
Passing a non-contiguous float matrix to CartesianRepresentation.transform() gives the wrong answer because the relevant lower-level function (erfa.ufunc.rxp()) accesses the matrix as if it were contiguous. (Curiously, this problem does not occur for a non-contiguous integer matrix, presumably because rxp() is given a float copy, and hence contiguous version, of the original matrix.)
Expected behavior
The correct answer
How to Reproduce
Create a matrix with a spurious extra column, and then slice the column away to get an identity matrix that is non-contiguous.
>>> matrix = np.array([[1, 0, 0, 5], [0, 1, 0, 6], [0, 0, 1, 7]], dtype='float')[:, :3]
>>> matrix
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
>>> matrix.flags
C_CONTIGUOUS : False
F_CONTIGUOUS : False
OWNDATA : False
WRITEABLE : True
ALIGNED : True
WRITEBACKIFCOPY : FalseTransforming any CartesianRepresentation by this identity matrix should return the same representation, but it does not, because rxp() uses the values in the hidden column:
>>> from astropy.coordinates import CartesianRepresentation
>>> CartesianRepresentation([1, 2, 3]).transform(matrix)
<CartesianRepresentation (x, y, z) [dimensionless]
(1., 8., 12.)>Versions
Windows-10-10.0.19045-SP0
Python 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:34:57) [MSC v.1936 64 bit (AMD64)]
astropy 5.3.4
Numpy 1.26.0
pyerfa 2.0.0.3
Scipy 1.11.3
Matplotlib 3.8.0
Activity
ayshih commentedon Oct 18, 2023
If I copy the matrix to make a contiguous version, I get the correct answer:
ayshih commentedon Oct 18, 2023
Incidentally, using
@instead ofrxp()works fine:mhvk commentedon Oct 19, 2023
Ouch. That's definitely a bug in pyerfa, so I copied your example to liberfa/pyerfa#123
mhvk commentedon Oct 19, 2023
I'm releasing a new pyerfa version for this - it affects anything that takes in matrices - and will make a PR with your example as a test case so that also on the astropy side we guard against this.
Thanks so much for reporting this, @ayshih, it must have been a very tricky one to uncover...
mhvk commentedon Oct 19, 2023
See #15511 for the regression test.
ayshih commentedon Oct 19, 2023
Heh, yeah, it was a wild trip narrowing down why refactored code was giving completely different answers, especially because at one point I mistakenly ruled out non-contiguity as a factor because I sanity checked with sliced integer matrices, where the bug doesn't actually happen.
Thanks for fixing this promptly!