@eteq - thanks for working on this! Just out of curiosity, how long does the erfa.pyx source take to be converted to C and also compiled? Does it slow down the build process significantly?

Mapping from the arcane ERFA/SOFA function names to something more readable would be nice

This could be a lot of extra work - an alternative is to simply have verbose comments before an ERFA call to explain what we are calling and why. In any case, I agree this can be done later. +1 to having the vectorized ERFA interface not be part of the public API for now.

In liberfa/erfa-wrapping-experiments#6, @mdboom asked something that should probably be discussed here:

Do we need to use the all-in-one file? I've never personally understood the motivation for that -- it makes using an upstream ERFA have an extra step that maybe we don't need to do either.

Not necessarily, now that you mention it. I think the main reason originally was that it makes it much easier to tell cython what to link to (just erfa.c instead of a long list of .c files). But maybe that's less important with the wrapper that has to be auto-generated anyway?

@astrofrog - The process of converting erfa to .pyx is pretty fast - less than a second for me. The cythoing/compiling is a lot slower (~5-10 sec for me), but it wasn't too much longer than erfa_time.pyx, so it's mostly trading the one for the other.

an alternative is to simply have verbose comments before an ERFA call to explain what we are calling and why.

Yeah, that's a good point. Not sure how well we'd enforce it long-term, but it is a good short-term solution.

@eteq - right, I guess my thinking is, let's see first which functions we end up using, and once that's done then we can more easily figure out where to concentrate efforts for providing functions with better names (no point in doing it for all ERFA functions).

@astrofrog - ok, yeah, that sounds good.

Any opinion on whether we should try to automate the erfa.pyx as part of the build process vs. tracking it in source code?

@eteq - I think it'd be better to add the pyx source to the repository, it's less that can break, and it allows us (if needed) to treat special cases as needed. Regarding @mdboom's suggestion of not having a single C file (and I imagine by extension, not a single pyx file), I'm +0.5 on the suggestion since I'm not a fan of 17,000 line files.

@astrofrog - it would be just one .pyx file the way it works right now even if we use the one-function-per-.c file method. It would be no additional work to split out the c files (that's what the original code assumed anyway), but it would be more to make different .pyx files.

I think actually that just one erfa.pyx is better because there'll be less Cythoning overhead (less boilerplate "create a new module" stuff going on), and it's all going to end up in one place in python anyway.

For the .c files it doesn't really matter in the end except that it's less work to maintain to just copy erfa over instead of dealing with the script that smooshes them together. (And a little bit more overhead in setup_package).

@eteq - ok, that makes sense. In that case, I agree with the separate C file/single pyx file approach. Fewer steps is better.

Yes -- we definitely want only a single .pyx. Each .pyx has all the overhead of initializing a Python extension module (both in file size and import time). Whereas the separate .c files has no real overhead in the final product.

Alright, I just updated this to use the multi-file version of ERFA instead of the single-file, as well as a fix that at least allows everything to import fine with erfa_time removed. I pushed this up such that we should get two travis builds. The first one is with erfa_time.pyx still in place, and the second one includes the final commit, which remove erfa_time.pyx and replaces it with a python module that just has some aliases.

From local testing, I think most of the problems in erfa_time are due to the change that array outputs now use an array dtype rather than having multi-dimensional arrays. If that's really all, the easiest option might be to just put in a few converter functions to work around that for now.

I'll try to get the check_errors functionality in ASAP. Once that's in, in principal this is good to go as we can always move forward temporarily with the old erfa_time.pyx.

Hmm, Travis ran both tests on the version that doesn't have erfa_time.pyx. I think maybe I really confused it somehow by putting two up at the same time? Anyway, the tests for the second-to-last commit pass for me locally. Will try to fix this up more in a bit.

@eteq, all - Sorry for the silence on this.

From local testing, I think most of the problems in erfa_time are due to the change that array outputs now use an array dtype rather than having multi-dimensional arrays. If that's really all, the easiest option might be to just put in a few converter functions to work around that for now.

This is the part that I wasn't too happy about. But as you suggest, I wouldn't wait for this to be fixed (at least by me) to move forward. This is even more critical for the matrix operations available in ERFA/SOFA.

The documentation at http://docs.scipy.org/doc/numpy/reference/c-api.iterator.html#NpyIter_AdvancedNew seems to indicate that the broadcasting rules can be tweaked. Ideally, we would want for a C function with signature (float[3], float[2], float) -> float to broadcast all dimensions except the rightmost one.

The record aspect seems not that bad, as it separates the dimenions of the number of items from what is contained within them. And with this, broadcasting already works as I would expect it!

In [33]: erfa.gst00b((24512345.+np.arange(2.))[:, np.newaxis], 0.123456789+np.arange(10.)*0.01)
Out[33]: 
array([[ 4.29258181,  4.35558571,  4.41858961,  4.48159351,  4.54459742,
         4.60760132,  4.67060522,  4.73360912,  4.79661302,  4.85961692],
       [ 4.30978663,  4.37279053,  4.43579443,  4.49879833,  4.56180224,
         4.62480614,  4.68781004,  4.75081394,  4.81381784,  4.87682174]])

Of course, a change is only a .view away! If we are happy with the record dtypes for now, I can make the required changes for Time. This will almost automatically solve #3039 as well...

But first some debugging... I found at least one error due to wrong casting. As this would be for liberfa generally, I'll raise this at liberfa/erfa-wrapping-experiments#6

On what to include in tar balls: I'm somewhat in favour of just including the actual source, not an intermediate product, i.e., go all the way from the .c files.

@mhvk - so just to be clear, you are saying we keep only the initial code in the repository, and only the final wrappers in the release tar balls? Nothing should get generated for users who download stable releases (like for Cython, where we include only the final .c files)

yes, sorry for being confusing, I would put only the initial code in the repository, not intermediate products such as the .pyx file. But I should add that I'm not speaking from a lot of experience, so don't take my opinion too seriously!

@mhvk - In theory I like your suggestion of not having the .pyx files. but, the problem is that it's not possible to do that without mucking around with the setup scripts. So my thinking is to store the .pyx file in source for now until someone gets around to automating the process (in a later PR).

If we don't do that, we have to tell anyone building from source to cd into the relevant directory, run the generator script, and then go back, which will likely be a big pain.

Now, for the changes I just pushed up: there's quite a bit there. What I've added now is:

1. Various cleanup, including a changelog entry
2. The status codes are now checked for non-0 values. If they are negative, that's an error (raises an exception), while positive issues a warning. That seems to be the convention in the ERFA code, although I can't guarantee that for sure.
3. I removed the record/structure stuff by simply grabbing the relevant field before returning the value (for outputs), and some recasting of arrays for the input.

With these changes, I think this is ready to go (modulo the next paragraph).

@mhvk - At first I thought the same, that maybe leaving these as structured arrays is a good thing. But after playing with it a bit, I started running into oddities, like it not being possible to do arr.view(dtype(['fi', 'd', (3,3))) - for some reason that doesn't seem to work if the field has more than one dimension (numpy just gives an opaque exception). Also, it fails in various non-intuitive ways if you slice/stride your array, or have something that's fortran-ordered (not impossible given the presence of fits files floating around...) So I concluded it's better to just use the final dimensions for vectors and the like, which is my 3 above. But that could be tweaked further or backed out if you don't like it...

Also, this PR currently retains erfa_time.pyx. But https://travis-ci.org/eteq/astropy/builds/38687174 has a build with it removed (from my erfa-cython-auto-no-erfatime branch)

@eteq - just to make sure I understand, does the fact the build without erfa_time.pyx fails mean that not all of astropy is using the new wrappers? Is the plan to do this in a separate PR?

@eteq - OK on leaving pyx in for now; not a big deal.

On the record arrays: recasting them, isn't Fortran order even more problematic if you decide to use the final dimenion? In a record array, the values in the record are stored right next to each other, so I would have thought that for Fortran order you'd need to use the first dimenion(s) to ensure you have the same view of the memory. For matrices, it may be even less obvious what to do, since the matrix itself is presumably always stored in C order...

With record arrays recast, Time should work better, and indeed I see the number of errors has been reduced substantially! Some of the failures that do occur, however, seem to indicate real problems. E.g., the doc test calculating UT1 gives one entry that is off by half a second.

013     >>> from astropy.time import Time
014     >>> t = Time(['2012-06-30 12:00:00', '2012-06-30 23:59:59',
015     ...           '2012-06-30 23:59:60', '2012-07-01 00:00:00',
016     ...           '2012-07-01 12:00:00'], scale='utc')
017     >>> t.ut1
Differences (unified diff with -expected +actual):
    @@ -1,3 +1,3 @@
    -<Time object: scale='ut1' format='iso' value=['2012-06-30 11:59:59.413'
    - '2012-06-30 23:59:58.413' '2012-06-30 23:59:59.413'
    - '2012-07-01 00:00:00.413' '2012-07-01 12:00:00.413']>
    +<Time object: scale='ut1' format='iso' value=['2012-06-30 11:59:59.913' '2012-06-30 23:59:59.413'
    + '2012-07-01 00:00:00.413' '2012-07-01 00:00:00.413'
    + '2012-07-01 12:00:00.413']>

Currently, is a validation done of the python wrappers? (e.g., run a version of the erfa test file)

But perhaps this suggests starting to use the wrapped erfa routines should indeed be separated from getting the package to compile properly.

@eteq - sorry, with the many files it seemed hard, but just looking at the commit doing record->regular array (eteq@25acb27), I see you do some tests on fortran order. Just to be sure, maybe add a test with matrix input/output?

Though I'm still not completely sure the record array approach isn't better. Do you have an explicit example of what goes wrong with it?

@mhvk - good idea on the matrix input/output - will add something for that and push it up soon (along with some fixes that will hopefully get the tests passing).

Currently, is a validation done of the python wrappers? (e.g., run a version of the erfa test file)

I am not running the full erfa test suite mostly because translating it to python would be a painful act: They aren't as consistent as the docs are, so it's probably not easy to automate the translation, so it would need to be done by hand for every function. There are some "spot checks" where I compare a few functions here and there to what the erfa tests say the answer should be (and they are passing). Will add some for ut1 and other places where time is being suspicious and see what happens.

Though I'm still not completely sure the record array approach isn't better. Do you have an explicit example of what goes wrong with it?

I agree with you in theory, but in practice a bunch of awkward things jumped out once I tried to actually use it this way:

• Converting to/from vectors isn't too bad, but for matricies it is difficult:

>>> dtv=np.dtype([('fi0','d',(3,))])
>>> dtm=np.dtype([('fi0','d',(3,3))])
>>> arrvec = np.array([[1,2,3], [4,5,6]], dtype=float)
>>> arrmat = np.array([np.eye(3)]*2)
>>> arrvec.shape
(2, 3)
>>> arrvec.view(dtv)
array([[([1.0, 2.0, 3.0],)],
       [([4.0, 5.0, 6.0],)]], 
      dtype=[('fi0', '<f8', (3,))])
>>> arrmat.view(dtm)
ValueError: new type not compatible with array.

• Even with vectors, view can do unexpected/surprising things silently:

>>> dtv=np.dtype([('fi0','d',(3,))])
>>> arrvec = np.array([[1,2,3], [4,5,6]])
>>> arrvec.view(dtv)
array([[([5e-324, 1e-323, 1.5e-323],)],
       [([2e-323, 2.5e-323, 3e-323],)]], 
      dtype=[('fi0', '<f8', (3,))])

It's not obvious, but the problem is that arrvec's dtype here is int, but we cast it to a double. The current scheme automatically takes care of that by using numpy's casting rules, but view won't do that. Similar problems happen for fortran ordering or striding, along the lines about what you were worried about in the first comment, @mhvk .

• Inputs to some of the ERFA functions would have to be structured array dtypes. But fields have to have a name, and you have to know the name to pull out the "real" array. (here, internally, it's always "fi0", but that's not guaranteed if users can provide their own dtypes with any name)
• Generically, I find them more difficult to work with because this sort of thing isn't really what structured arrays were meant for. To create a new one without the view problems above, you need to do a three step process of creating the dtype, creating an empty array, and then populating it from another array. It's much easier to just do np.eye(3) or whatever and manipulate that with all the tools numpy provides.

It may be there's some clever trick to get around all of that, but this just seemed easier.

@eteq - yes, view is not always obvious, though there are of course solutions:

In [10]: arrmat.reshape(-1, 9).view(dtm)
Out[10]: 
array([[([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],)],
       [([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],)]], 
      dtype=[('fi0', '<f8', (3, 3))])

(fastest dimension has to match the total length of the field)

And one can get the name via dtype.names:

In [46]: arrmat.reshape(-1, 9).view(dtm)[dtm.names[0]]
Out[46]: 
array([[[[ 1.,  0.,  0.],
         [ 0.,  1.,  0.],
         [ 0.,  0.,  1.]]],


       [[[ 2.,  0.,  0.],
         [ 0.,  2.,  0.],
         [ 0.,  0.,  2.]]]])

Some more tricks at http://wiki.scipy.org/Cookbook/Recarray

But as you say, it may not be worth it to stick to the record arrays; as long as it is always the same set of indices that iterate over input and output, it should be fine (and I think I convinced myself there is no fortran issue here...).

Rebased on master and it like the tests have passed, so this should be set to merge!

Once this is in I can add an issue about cleaning up time's aliases and making it multidimensional. @taldcroft - before hitting the merge button I just want to make sure you're ok with this as far as time is concerned. I think all others have signed off, so feel free to hit the merge button if you're happy with this as-is.

@eteq - sorry for the delay, I've been at a meeting. I just merged this. 👏

Fantastic, thanks @taldcroft !

🎉

This is very cool, guys. Good work. I'll be using this stuff immediately.

Now I can just:

from astropy import erfa

So... 😄