Wrap erfa functions with ufuncs and gufuncs by mhvk · Pull Request #7502 · astropy/astropy

mhvk · 2018-05-25T03:37:18Z

This is an alternative to #7484, with the main difference being that for matrices and vectors now the more logical "generalized ufuncs" are used, which take care of the trailing dimensions of (3,) and (3,3), respectively. EDIT: improved performance is just like for #7484.

In this version, I decided on purpose not to similarly treat pv (trailing dimensions of (2, 3)) and integer decompositions of years, hours, and degrees (ints with trailing dimension of (4,)). I think those really are more logically treated as structured dtype, since their parts have different units - it is really bad to have to remember that pv[..., 1, :] is the velocity part, with units that can be au/day or m/s depending on context - this is what we have Quantity for! But it will be very difficult to have different units for different parts of an array; in contrast, it should be quite trivial to have different units for different parts of a structured dtype. And seeing pv['v'] it is much more obvious what one is looking at.

In the end, I found that surprisingly few changes were required to the code to make this work, and I feel all add to the clarity. Hence, I prefer this PR over #7484.

p.s. The only thing I like a little less is the work-around I need to do for those functions which produce vectors or matrices, but do not take them. This is because currently gufunc signatures cannot indicate that some dimension should be a fixed number (being addressed... EDIT: numpy/numpy#11175). This means I have to generate an empty output with the right dimensions. But this is unhandy for Quantity, so right now I work around that by making the output of a special class. Not too elegant, but works.

astropy-bot · 2018-05-25T03:37:29Z

Hi there @mhvk 👋 - thanks for the pull request! I'm just a friendly 🤖 that checks for issues related to the changelog and making sure that this pull request is milestoned and labeled correctly. This is mainly intended for the maintainers, so if you are not a maintainer you can ignore this, and a maintainer will let you know if any action is required on your part 😃.

Everything looks good from my point of view! 👍

If there are any issues with this message, please report them here.

mhvk · 2018-05-30T20:51:37Z

@astrofrog - I posted a gist with only a small number of erfa functions wrapped - but ones that show the things that are needed to deal with different cases: ~~https://gist.github.com/mhvk/964a26183f92ab3f928d39f2b3bc08ec~~ EDIT: https://gist.github.com/mhvk/a3c17f0db77b91151d50a0d278621d97

mhvk · 2018-05-30T21:00:34Z

Sorry, that should be https://gist.github.com/mhvk/a3c17f0db77b91151d50a0d278621d97, as I should at least include one function that needs the fix for when there are only output arguments with a core dimension of 3 (fk5hip).

mhvk · 2018-05-30T21:03:12Z

More general comment on gists: note how for most functions, the python wrapper does essentially nothing but combine status codes. The only ones where more is done is for signatures where the output has some core dimension (e.g., is a vector) and the input is not, since in that case the gufunc needs an output argument to know what the dimension should be.

It is actually very annoying that one cannot give a fixed size, but this is being fixed in numpy/numpy#11175

Using ufuncs instead of C routines and a custom iterator has the advantage that casting to arrays with the right dtype and iteration is taken care of automatically, in numpy C code. For small arrays, this means substantially faster execution. Another advantage is that classes can now override the ufuncs using __array_ufunc__ - this is important both for ensuring Quantity input is disallowed (or at least treated correctly rather than just cast to array), and that a new Distribution class can work with the erfa functions.

Before, contents of astrom and ldbody were kept as well.

Also simplify the viewing back.

Rest should be kept as structured arrays.

mhvk · 2018-06-14T16:44:41Z

I realize this is a difficult PR to review, but it would be nice to get this in - substantial performance increase, plus more to come with things that are enabled by this. @astrofrog - do you think you will have time to review? Or do you perhaps have a suggestion for someone else?

astrofrog

I've just reviewed this - I don't understand 100% of the C code, but overall this looks great to me. Just a few minor comments below.

I know that @lpsinger has experience with writing ufuncs - @lpsinger - would you mind taking a look at this? (note that you can view the generated code if you need, see previous comment by @mhvk). In particular, I would be curious to know if there is anything that is being done here that would preclude using OpenMP in future to speed things up?

astrofrog · 2018-06-17T20:31:32Z

-    except:
-        raise ValueError("{0} must be of trailing dimensions {1}".format(name, shape))
+class D3Fix(numpy.ndarray):
+    # We do this to help ensure Quantity happily lets this pass...


Could you provide more details in the comment here? I don't quite understand why this is needed.

Yes, how could one... I added a very detailed comment.

astrofrog · 2018-06-17T20:32:29Z

+
+
+def arrayify_inputs_and_create_d3_fix(inputs, core_dims, out_core_shape, out_dtype):
+    # Need to ensure we can get a shape


This comment isn't very clear - could you write a comment to describe more verbosely what this function does?

Now a proper docstring.

astrofrog · 2018-06-17T20:33:36Z

+# <--------------------------Actual ERFA-wrapping code------------------------>

 {% for constant in constants %}
 {{ constant.name }} = {{ constant.value }}


The generated code here has brackets that aren't needed:

DR2AS = (206264.8062470963551564734)

This appears to be the case only for floating point values (not integers)

The parentheses are in the erfam.h file. I don't think it hurts to just let them be.

astrofrog · 2018-06-17T20:38:48Z

+
+/* Extra/PVMergeExtract */
+void eraPav2pv(double p[3], double v[3], double pv[2][3]);
+void eraPv2pav(double pv[2][3], double p[3], double v[3]);


Could you add a comment to this file about why it is needed?

astrofrog · 2018-06-17T20:40:55Z

+                or (subsection == "AngleOps")
+                or (subsection == "SphericalCartesian")
+                or (subsection == "MatrixVectorProducts")
+                or (subsection == 'VectorOps')):


The inner parentheses are not needed here

I left them as I found them.

But talking about those: shall just take out the whole clause? We're compiling nearly everything anyway, so there is little point in terms of trying to keep it small or so.

I'll do that as a second commit, so it would be easy to revert.

astrofrog · 2018-06-17T20:41:33Z

+    # pvmatarr = np.array([pvmat]*3, order='F')
+    # astrom4 = erfa.apcs13(jd1, jd2, pvmatarr)
+    # assert astrom4.shape == (3,)
+    # np.testing.assert_allclose(astrom4['em'], 1.010428384373318379)


Was this temporarily commented out?

Some of those tests no longer make sense, since the pv vectors are not input as (..., 2, 3) arrays. I'll keep what does work still.

Also uncomment tests that were not passing in an earlier stage.

astrofrog

Looks good to me know. I'll leave it up to you as to whether to wait for another review.

mhvk · 2018-06-19T00:40:24Z

OK, thanks. @jwoillez, @lpsinger - if you plan to review, could you let me know? I think it is actually all OK, so if I don't hear in the next few days, I'll merge, so I can go on and use the erfa routines in coordinate transformations.

Cadair · 2018-06-22T15:28:31Z

Merging this into master appears to have broken astropy.coordinates.get_body_barycentric:

>>> get_body_barycentric("earth", Time("2018-01-01T00:00:00"))

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-2-d4041245237a> in <module>()
----> 1 get_body_barycentric("earth", Time("2018-01-01T00:00:00"))

~/Git/astropy/astropy/coordinates/solar_system.py in get_body_barycentric(body, time, ephemeris)
    362     """
    363     return _get_body_barycentric_posvel(body, time, ephemeris,
--> 364                                         get_velocity=False)
    365 
    366 

~/Git/astropy/astropy/coordinates/solar_system.py in _get_body_barycentric_posvel(body, time, ephemeris, get_velocity)
    241 
    242         body_pos_bary = CartesianRepresentation(
--> 243             body_pv_bary['p'], unit=u.au, xyz_axis=-1, copy=False)
    244         if get_velocity:
    245             body_vel_bary = CartesianRepresentation(

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

pllim · 2018-06-22T15:31:21Z

@Cadair , best to open a new issue then!

lpsinger · 2018-06-26T20:37:37Z

Whoa. My mind is blown by using jinja2 to generate ufunc wrappers. Your Numpy-fu exceeds my own.

mhvk · 2018-06-27T03:03:01Z

@lpsinger - the jinja2 stuff is really just an adaptation of what @jwoillez came up with... But it is pretty amazing!

mhvk added time coordinates Performance erfa Refactoring labels May 25, 2018

mhvk added this to the v3.1 milestone May 30, 2018

mhvk changed the title ~~Wrap erfa functions with ufuncs, version 2~~ Wrap erfa functions with ufuncs and gufuncs May 30, 2018

mhvk mentioned this pull request May 30, 2018

Wrap erfa functions with ufuncs #7484

Closed

mhvk mentioned this pull request May 30, 2018

ENH: Generalized ufunc signature expansion for frozen and flexible dimensions numpy/numpy#11175

Merged

mhvk requested review from astrofrog and jwoillez May 30, 2018 21:04

mhvk force-pushed the erfa-wrapped-by-gufuncs branch from c9a55c0 to df37662 Compare May 30, 2018 21:08

mhvk added 14 commits May 31, 2018 11:34

Do not pass quantity to erfa in time.

4aa2101

Adjust coordinate docs for the change in error message.

cb475af

Add changelog entry.

e8d2de3

Simplify erfa_generator, ensuring only real arguments are kept.

63f2ce7

Before, contents of astrom and ldbody were kept as well.

Give the special dtypes more descriptive names.

0a55c5c

Also simplify the viewing back.

Add also the dtype for pvdpv, in case we ever compile all of erfa.

877f431

Write any multi-d argument with a gufunc.

f027737

Temporarily create all functions

202f1e2

Ensure gufuncs are used for vectors and matrices.

49a3eea

Rest should be kept as structured arrays.

Ensure core dimensions are 3 even if it requires blatant hacks.

4ebeb09

Less hacky way to deal with outputs with core dimensions.

414cbab

Fixes to Time and Coordinates because of structured erfa output.

af8106f

Changelog entry fix for this PR number

43ffab7

Only compile erfa routines actually needed.

370a036

mhvk force-pushed the erfa-wrapped-by-gufuncs branch from df37662 to 370a036 Compare May 31, 2018 18:15

astrofrog requested changes Jun 17, 2018

View reviewed changes

mhvk added 2 commits June 18, 2018 11:54

Document some of the more obscure parts better.

f0f6f4a

Also uncomment tests that were not passing in an earlier stage.

Compile all of erfa.

912da17

astrofrog approved these changes Jun 18, 2018

View reviewed changes

mhvk mentioned this pull request Jun 22, 2018

BUG: better detect memory overlap when calling gufuncs numpy/numpy#11381

Closed

mhvk merged commit 91fc183 into astropy:master Jun 22, 2018

mhvk deleted the erfa-wrapped-by-gufuncs branch June 22, 2018 11:17

Cadair mentioned this pull request Jun 22, 2018

get_body_barycentric is broken #7591

Closed

mhvk mentioned this pull request Jul 2, 2018

Make erfa functions ufuncs instead of using custom wrappers? #6865

Closed

mhvk mentioned this pull request Jul 19, 2018

Unknown ufunc apio13 #7665

Closed



		def arrayify_inputs_and_create_d3_fix(inputs, core_dims, out_core_shape, out_dtype):
		# Need to ensure we can get a shape

Conversation

mhvk commented May 25, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

astropy-bot Bot commented May 25, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented May 30, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk commented May 30, 2018

Uh oh!

mhvk commented May 30, 2018

Uh oh!

mhvk commented Jun 14, 2018

Uh oh!

astrofrog left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

astrofrog left a comment

Choose a reason for hiding this comment

Uh oh!

mhvk commented Jun 19, 2018 • edited by astrofrog Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Cadair commented Jun 22, 2018

Uh oh!

pllim commented Jun 22, 2018

Uh oh!

lpsinger commented Jun 26, 2018

Uh oh!

mhvk commented Jun 27, 2018

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

mhvk commented May 25, 2018 •

edited

Loading

astropy-bot Bot commented May 25, 2018 •

edited

Loading

mhvk commented May 30, 2018 •

edited

Loading

mhvk commented Jun 19, 2018 •

edited by astrofrog

Loading