What should be the calling convention for ufunc inner loop signatures?

Mentioned in a bunch of places, but would be valuable to have a single place to track this.

To be able to reuse an inner loop across multiple dtypes (ie across variable-length strings, structured arrays, etc), we need access to the dtype object within the inner loop.

Perhaps we can achieve this by adding a new npy_bool innerloop_takes_dtype flag to ufuncs, which when set indicates that the inner loop signature now also takes a PyArray_Descr **dtypes argument.

If it's possible it would be nice to avoid piling new flags and if statements into the ufunc code. I'm not sure if it's possible.

There was consensus a few years ago about hiding the ufunc struct internals from the C API. At least at that point, the only user was numba.

If that's still doable, then we can hopefully reimplement the current register-a-ufunc-loop function as one that takes an old-style function pointer, wraps it in a new-style function pointer, and then register that. That way the core ufunc code will be able to assume that all loops are new-style.

There was consensus a few years ago about hiding the ufunc struct internals from the C API

This gets a +1 from me. I hadn't considered the fact that we can wrap one into the other.

Another option... In C, is it safe to call a function pointer with more arguments that it expects? Won't it just put extra items on its stack that it can ignore?

Is there consensus that this is the right approach? Because my gut feeling is that it might not be.

Basically, what operation would you want the dtype for in a more generic loop where it is not better to do it either before or after the all of the inner operations? My gut feeling would rather be:

Run dtype specific Loop Selection/Preparation function before actual calculations:
- This can error out (say a unit error)
- It could in principle even return the inner loop
- It would set/create the new dtype object.
- It can return a void * object to be passed in as additional data to the ufunc, which could be, but does not have to be the dtype.
Possibly run another function after the inner loop is finshed (for example to allow error handling, since the inner loop should be able to unlock the GIL).

I find it a bit strange to specifically pass in the dtype object to the inner loop, because I think we need a "setup" step anyway, and that setup step can decide what it wants to pass in. Of course for some "generic" inner loops (wrapping python objects or so), they may end up passing in the dtype, but that would be implementation.

Not sure I make sense, but it seems to me that almost all logic derived from the dtype inside the inner loop would rather already be handled before calling the inner loop.

EDIT: Maybe to add to it. If dtypes are not created new (but reuse the old one) and we pass them in. Then it seems to me that tagging errors to them will be tricky in a thread safe manner in any case.

I'm also not sure we want to pass on dtype by default, since in general the inner loops should be simple and fast, and not have to deal with details like it. Even for structured dtypes, I'm not so sure I understand the general case; even for equality, if it weren't for being unable to override the ufunc, the current outside implementation would seem the proper solution (which suggests we just need to move some of the logic from ndarray.__eq__ to one of the setup functions in the np.equal ufunc).

For strings, the case seems clearer, but it does seem that the machinery is in principle present to deal with it: we should be able to use the void * for passing on, e.g., length information.

All this said, it does seem worth checking exactly where those helper functions get called, so that we can make sure, e.g., it is possible to return a custom ufunc loop or to create an output array with a custom dtype.

We probably should rework the ufunc loop APIs in general to support setup/teardown/dispatch in a more logical way. (For example, np.add on strings will need some way to compute the output dtype before the actual loop starts running.)

But we should still pass the dtype information to every loop. Forcing people to jump through hoops to access this basic information makes no sense to me. Writing np.equal on strings shouldn't require writing a special setup trampoline that just copies a pointer from one place to another.

Possible new ufunc loop signature:

nin
nout
dtypes (array of dtype object pointers with nin+nout entries)
data (array of char*'s with nin+nout entries)
outer_count: how many times the core operation should be performed (single integer)
outer_strides: how far to step between each core operation (array of integers with nin+nout entries)
core_ndims: the dimensionality of each core argument (array of integers with nin+nout entries, all zero for scalar ufuncs)
core_shapes: the shapes of each core argument (use ndims to decode; can be NULL if core_ndims is all-zero)
core_strides: the strides within each core operation (same layout as core_shapes)
void *per_call_user_data: arbitrary data passed from the loop setup function
void *per_loop_user_data: arbitrary data passed when registering the loop function.

Here's a use case that might clarify things. At some point we'll want to expose the loop machinery to python, so people can write loops implemented in python, for prototyping and for all the stuff people use vectorize for.

So: given a proposal for the C loop signature, how hard would it be to write a single C loop function that can be registered with any ufunc, and turns around and invokes a python function instead? Obviously this is much easier if the loop has information on dtypes, so it can convert the raw data pointers into ndarrays to pass to python.

Good point, all flexible dtypes will need it in some form or another, that is definitely true. It is just, for downstream it seems to me that improving setup/teardown is more important. For the support of flexible dtypes, passing the dtypes would be the most convenient thing for sure.

So I am not against it, it will hardly hurt and is indeed more convenient in many cases. Too many flags (or typically unused args) may also not be convenient though, so before that happens, I am not sure if a minimal solution may not be just as well. But, we will need an iteration on the inner loop API anyway probably, so there is probably no point to trying to make due with the current signature.

I am sure it is already on everyone's list, but one other API change is that a stop iteration/error integer return value may be a good idea.

@njsmith

@njsmith well, but you need to write a single C-loop+setup func, I think we could live with that single loop /setup func being inconvenient to write :).

Anyway, I am not against it. I just thought it might be focusing on the wrong thing. But that is probably just me not being up to the full picture.

I am sure it is already on everyone's list, but one other API change is that a stop iteration/error integer return value may be a good idea.

Good point! Loops can raise exceptions already, so we might as well make things less weird and have them return a -1 to signal that, like every other python C API. And it would let us stop manually checking for exceptions after every call.

Re: error handing, the other thing to think about is the floating point error flags. Currently we clear the flags before every call, and check them afterwards. This is pretty expensive on some CPUs, so it would be nice to be able to skip it for things like integer or string operations where it's irrelevant.

I guess the simplest way to do this would be to declare that new-style loops are responsible for doing their own checking for floating point errors, and providing a simple API to help them do this (e.g. call PyArray_FloatErr_Start at the beginning of your loop, and PyArray_FloatErr_Check at the end).

Oh man, there's also the GIL handling to think about. What even is the current behavior? Something like, loops are called with the GIL dropped, unless there's a special flag set on one of the dtypes?

But the basic argument for needing a dispatch table in the ufunc object is very simple. Suppose you're writing a third-party library and want to define a custom ufunc that works on floats, like scipy.special.erf. Presumably you aren't going to go mutate numpy's float64 dtype to make its __dtype_ufunc__ method know about scipy.special.erf. So the ufunc has to have a dispatch table in it, at least for built-in types. And if you want to avoid having "two different notions of dtype", then that should be how all dispatch is done, rather than using one dispatch mechanism for built-in dtypes and a different one for user-defined dtypes.

This may be the fundamental design tradeoff. Either:

We have two different types of ufunc dispatch: "for dtypetype" and "for dtype given a dtypetype"
We have two different types of ufunc dispatch: defined on ufuncs and defined on dtypes.

I don't really like either of this options but it seems like we'll have to pick. Either way, in practice we will be able to reserve the second type of dispatch for unusual cases.

For what it's worth, I agree that all things being equal putting dispatch in ufuncs makes more sense. I would try to save dispatch on dispatch on dtypes only for new dtypes / old ufuncs.

The biggest virtue of supporting defining dispatch on dtypes might be that you could potentially define generic implementations of every ufunc, e.g., by casting to another dtype and then casting back. For example, this would probably be the easiest way to implement a variable length string dtype in NumPy, by doing all computation on object dtype arrays. (That said, it would probably work nearly as well to expose a generic override for implementing the inner loop of "every built-in arithmetic operation".)

Also, what if I want to define my own ufunc that works on datetime64?

You would have to define casting rules for handling inputs with different dtypes. The difference is that you would do this by providing a "default dtype resolver" function, rather than by specifying a dtypetype implementation + "dtypetype specific dtype resolver".

datetime64 is actually an easy case: you can rely on NumPy's builtin casting rules for datetime64/timedelta64 to give the same units for everything, e.g., by making use of a protocol version of np.result_type. Or if you're pandas, you can only bother to implement your functions for datetime64[ns].

I guess a virtue of "dispatching on ufuncs" is that you can get this sort of casting for free, so much of the time you don't need to write the "dtypetype specific dtype resolver." But it's probably not that much more difficult to make a generic "dtype resolver that uses casting."

Cases like new ufuncs for a custom units dtype will be tricky either way.

We have two different types of ufunc dispatch: "for dtypetype" and "for dtype given a dtypetype"

Sorry, maybe I'm being tired, but I don't know what these two types are :-). Can you explain?

We have two different types of ufunc dispatch: "for dtypetype" and "for dtype given a dtypetype"

Sorry, maybe I'm being tired, but I don't know what these two types are :-). Can you explain?

"for dtypetype": This is what currently have in NumPy, i.e., the loops associated with each ufunc. These can't handle metadata.

"for dtype given a dtypetype": This is the new thing we would need, to pick the specific dtypes/loop (e.g., (string[2], string[3]) -> string[5]) based on the dtypetypes ((string, string) -> string)) and the ufunc (np.add)

The "dispatch" part (i.e., picking which loop to use) is identical between those: in both cases you look at the input [1] dtypetypes and find the loop that matches (possibly after casting).

The difference comes after you pick which loop you want to use, and only depends on the output dtypetypes. If one or more of those dtypetypes requires some parameters to instantiate, then the loop needs a mechanism to specify those parameters. OTOH if all the output dtypetypes don't require any parameters (they're effectively singletons), then this mechanism collapses into triviality.

Calling this "two different types of ufunc dispatch" feels weird to me, because one is a special case of the other, and anyway, the difference happens after the dispatch phase :-). And anyway... the reason we're forced to add the concept of a "loop setup function" is to handle output dtypedtypes that require parameters to instantiate, but once we have the concept, it probably has other uses, even for non-parametric dtypetypes (e.g. precomputing things).

[1] I think in some cases currently, if an out= argument is given, we might also take its dtype into account when picking which loop to use, which means that passing out= arguments can change the output. If we can get rid of this I think that would be good.

The "dispatch" part (i.e., picking which loop to use) is identical between those: in both cases you look at the input [1] dtypetypes and find the loop that matches (possibly after casting).

I'm not sure we can pick loops merely based on dtypetype. Some parametric dtypes might correspond to different inner loops even with the same dtypetype, e.g., consider dtypes like Masked or Quantity that could wrap any number of low-level dtypes.

I guess this might not be a meaningful distinction if the loop has access to dtype metadata, e.g., so it can call the ufunc again on NumPy arrays if necessary.

I was going to reply to your first paragraph, but your second paragraph already said what I was going to say :-)

Though actually, this does suggest: maybe our ufunc inner loop signature should also receive the ufunc object? Similar to how methods receive self?

Maybe we need a pure python implementation we can play with to see how this works. It should be simple without too many decision points.

One of the detriments of the current mechanism(s) is that the decisions are made far away from where they are used, making debugging difficult. Being able to stop in a debugger and walk back up a linear stack to see how choices were made is a definite plus.

Though actually, this does suggest: maybe our ufunc inner loop signature should also receive the ufunc object? Similar to how methods receive self?

Sure, though I suppose you could also easily put this into void *per_loop_user_data if desired.

In practice I think most new dtypes will be implemented via high level Python interfaces (if possible), so it doesn't matter as much what we do for the lower level C interface.

Maybe we need a pure python implementation we can play with to see how this works.

It seems to me like we're trying to design a fully-fledged function-overloading system here, and trying to cover many of the parametric cases that C++ templates would cover. This seems like it might be in peril of falling outside the scope of numpy, although I'll admit that custom dtypes lose significant value without it.

Yes, we are designing a fully-fledged function-overloading system in a sense. I doubt we can avoid that, I am not even quite sure how much we already have? It would be nice to keep it very simple. I think it might be ok to remove implicit upcasting (do not cast "ee->e" to "ff->f", since there is no overlap at all; we can add/wrap the loops manually for backcompat).

Considering the use case of making progress with flexible types, maybe we can side step most of this discussion if we can settle the signature change.

The idea of a self could be neat. For the python side it seems convenient, but that does not mean it has to exist on the C-side. I am not sure there is much of a use case on the C-side, since the inner loop should be simple and usually unique to a single loop.

Though... arrays with missing values might be a problem here, now that I think about it :-/

For units, it's not a big deal, because units behave differently in different operations (e.g. add has to check unit compatibility, multiply creates new units, etc.), and part of defining a units package is making a list of which operations you want to support and writing code to support each of them. So that fits well with the model of attaching loops to each of the ufuncs you want to support.

But for missing values, you want to essentially support all ufuncs that have support for the wrapped dtype. And the details are complicated – I'm not sure a simple __dtype_ufunc__ even solves it, once you consider cases like a ufunc that takes two arguments and both of them are passed as different missing-value dtypes, and that if the inputs are missing-value dtypes then the output has to be too.

I think the ufunc machinery will have to have some hard-coded knowledge of missing values in any case (e.g., we'll want ufuncs to support a generic skipNA= kwarg). Maybe the thing to do is worry about these details later, when figuring out missing value extensions specifically, rather than trying to make the basic dtype machinery generic enough to handle missing values as a special case.

I have started to work on a draft here #16723

The current plan is much like some of the posts above, although details around how to "reuse" loops best seem a bit complicated to me.

Promotion: I am tending towards multiple-dispatching like promotion, we could also have a list of promoters, etc. but you need some way to "order" them probably. The biggest advantage of proper/multiple promoters, is that a Promoter could dynamically create new loops.
For the general loop signature, I think we can limit ourselves to "strided" loops for now, possibly with allowing loop specialization (e.g. if one of the strides is 0 or contiguous; which are things that casting currently does). I think the normal (and likely only) inital way should be to just give a single strided loop though.
On the actual signature change:
- I think it is clear we want a return value. I am considering using a nonzero return to break the loop. A negative return value would indicate that a Python error is already set.
- I find it hard to give a reasonable set of things that should be passed in, so I am considering to rather pass in a structure with the following fields:
  1. ufunc
  2. descriptors *
  3. threadstate (so loops can safely set an error without dancing around it in setup/teardown)
  4. static loop data (current void *, or maybe some other Python object)
  5. void *per-run loop data
  6. "self"? I would not call it a true self, this is not really an instance or class. But there is some information, such as the identity which lives on the UFuncImpl/Loop. Some notion of "self" (or at least additionally passed in object) probably needs to exist for Python, and it being special to some degree may be necessary for easier refcounting.
As you can see that list is long and I am not even sure that there may not be something I missed. A struct at least could grow in theory. Plus, many numeric loops actually do not need any of the information (because we can handle the floating point error management), so this means they do not have to list 5-6 unused parameters.
- I think we need setup/teardown as I mentioned before. But, I am leaning towards a simple independent flag for floating point error handling (it is pretty special after all). In fact, while I think we should have setup/teardown, and should expose it eventually. I am happy not to rush that.

The signature decisions have overlap with casting (although casting does not have a clear "ufunc" that can be passed in). Casting is almost a ufunc, but not quite (possibly they could use the same but just pass NULL for ufunc, and reconsider later).

I think most of this has been incorporated into NEP 43 which was merged (in draft form).

This is interesting from a historical perspective, but for now there is no action item here with the current state.

Navigation Menu

Repository navigation

What should be the calling convention for ufunc inner loop signatures? #12518

11 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

What should be the calling convention for ufunc inner loop signatures? #12518

Description

Activity

njsmith commented on Dec 9, 2018

eric-wieser commented on Dec 9, 2018

seberg commented on Dec 10, 2018

mhvk commented on Dec 10, 2018

njsmith commented on Dec 10, 2018

njsmith commented on Dec 10, 2018

seberg commented on Dec 10, 2018

seberg commented on Dec 10, 2018

seberg commented on Dec 10, 2018

njsmith commented on Dec 10, 2018

njsmith commented on Dec 10, 2018

11 remaining items

shoyer commented on Dec 11, 2018

njsmith commented on Dec 11, 2018

shoyer commented on Dec 11, 2018

njsmith commented on Dec 11, 2018

shoyer commented on Dec 11, 2018

njsmith commented on Dec 11, 2018

mattip commented on Dec 11, 2018

shoyer commented on Dec 11, 2018

eric-wieser commented on Dec 11, 2018

seberg commented on Dec 11, 2018

njsmith commented on Dec 12, 2018

seberg commented on Jul 1, 2020

mattip commented on Oct 15, 2020

seberg commented on Sep 30, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions