ENH: New DTypes for numpy - #6

@jorisvandenbossche from pandas might be interested

It should be noted that in my thoughts, I have ignored object dtypes a bit for now. Those may need a few additional hooks (similarly, I did not think about adding zeros or so logic). Overall, I think we may be fine with only allowing very specific object dtypes which cannot override some of the most basic things (such as how data is copied). I very much want object[specfic_python_type] to be a thing, and for some types that could have faster ufuncs (and possible some specifics defined on the dtype itself).

@seberg I just caught up on the numpy mailing list thread on this topic, am reading through the hackmd page you wrote. Which of these is the preferred place for discussion? Mainly I'm thinking of how pandas could use this.

Here is fine, there is also a PR at numpy, which is a bit more up to date than the HackMD in places, but maybe the HackMD is better and I should just keep them more or less in sync. You can certainly comment inside the hackmd , although that is a bit easy to miss, and thus slow maybe.

Email will work (me or mailing list), but maybe github PR/issue is a bit nicer. I suppose we could also just open a discussion issue on NumPy. I could also schedule another video conference meeting soon. (It will have to be within about the next week, after that I will be unavailable for a month).

Oh, we also have weekly community calls from BIDS which will happen in about 1.5 hours. We could extend/talk a bit more at the end (I expect it will not be super long today, but since these are so big discussions would prefer to do it at the end when everyone can excuse themselves).

EDIT: Meeting info https://hackmd.io/76o-IxCjQX2mOXO_wwkcpg?both

Sounds good. Will stick with here for now.

Some thoughts/questions on what would make this really useful on my end (likely overlap with what TomAugspurger has said):

1. Would this make the non-python parts of the numpy codebase easier/harder to navigate for those with limited C-fu?

2. Would/could cython recognize dtypes implemented with the C-API?
   a) See e.g. https://github.com/pandas-dev/pandas/blob/master/pandas/_libs/algos.pyx#L691 where we handle M8[ns] by passing it as i8 along with a timelike=True kwarg, so that it treats NaT appropriately. If cython understood M8 directly, we could simplify a bunch of paths like this. (Note: since this isn't actually about new dtypes, I'd also be ecstatic if there is already a way to do this I'm just not aware of)

1. Would this make the non-python parts of the numpy codebase easier/harder to navigate for those with limited C-fu?

Tricky question. I think it is a pretty even thing probably, however, the plan now is to use slots/__dunder_methods__ more. Plus the plan is to add a Python interface (although I am not focusing on that right now). I think once there is a python interface to add dtypes, that helps a lot as an entry point into understanding the C-API.

2. Would/could cython recognize dtypes implemented with the C-API?
   a) See e.g. https://github.com/pandas-dev/pandas/blob/master/pandas/_libs/algos.pyx#L691 where we handle M8[ns] by passing it as i8 along with a timelike=True kwarg, so that it treats NaT appropriately. If cython understood M8 directly, we could simplify a bunch of paths like this. (Note: since this isn't actually about new dtypes, I'd also be ecstatic if there is already a way to do this I'm just not aware of)

This is a pretty big issue, that I have not given much thought. The problem with these things is that the buffer interface (and C) do not define these types. We can tag on additional information, and I think we should probably make it easy for cython to ask for the "C loop function of the ufunc", but that still is going to be annoying.

It is an interesting issue, I am tempted to claim it is a bit orthogonal, but that does not help much. I do not mind (allowing) tagging on information to dtypes if that makes life easier. But I am not sure what that would need to be.

I am tempted to claim it is a bit orthogonal, but that does not help much

Totally fair, declaring things orthogonal is the best way to keep big things like this manageable.

A few more thoughts after reading through the hackmd document:

• Some of the casting stuff reminds me of pandas.core.dtypes.cast.maybe_promote (which is kind of a mess and subject of several active PRs).

• You mention CyberPandas a couple of times as an example of a pandas ExtensionArray (EA) subclass. The other examples I'm aware of are fletcher, geopandas, (WIP) pint-pandas, and several test cases in pandas.tests.extension (also internally Categorical, PeriodArray, DatetimeArray, TimedeltaArray, IntervalArray).

• I think we've been pretty happy with the dtype registry pattern.

• Not long ago I looked at the quarternion implementation as a reference to try to implement my my own np.dtype subclass and got absolutely nowhere. If making this easier is part of the goal, then please consider providing examples/recipes/etc.

• The use cases for downstream authors to use this have a lot of overlap with the EA use cases. What part of the ExtensionDtype API could be handled by the new numpy implementation you're proposing (backward compat and versioning notwithstanding...)

• Not sure how helpful I can be, but feel free to reach out. I'm going to need a stockpile of goodwill with numpy folks to get help with my numpy wishlist.

Thanks for taking the time and reading it!

• Some of the casting stuff reminds me of pandas.core.dtypes.cast.maybe_promote (which is kind of a mess and subject of several active PRs).

Can you point me to this? Casting is tricky with flexible dtypes, but even if it is not perfect, I am currently fairly happy about the state. Although, I am considering adding fast-path methods for the can_cast check (i.e. users do not have to implement them, but can).

To be honest, I do not remember what Jeff Reback told me about it, but I cannot actually find how casting is currently defined in pandas?

• You mention CyberPandas a couple of times as an example of a pandas ExtensionArray (EA) subclass. The other examples I'm aware of are fletcher, geopandas, (WIP) pint-pandas, and several test cases in pandas.tests.extension (also internally Categorical, PeriodArray, DatetimeArray, TimedeltaArray, IntervalArray).

Thanks, will be useful to expand on these things a bit, before the end...

• I think we've been pretty happy with the dtype registry pattern.

This is a global type -> dtype mapping for coercion, or the python code to create a new extensiondtype/array? Just to be sure I got a clear idea, could you point me to the documentation for that?

• Not long ago I looked at the quarternion implementation as a reference to try to implement my my own np.dtype subclass and got absolutely nowhere. If making this easier is part of the goal, then please consider providing examples/recipes/etc.

Yeah, we need that. There is a wide range of use cases, from just creating something like quaternions (which is practically a specialized structured array) to completely new dtypes in C that also need scalars. It would be nice to make (common) things like quaternions super simple.

I think even currently, for someone who knows the python C-API well, the current dtype API is a bit limited. Mind you, I actually do not think it is super difficult, but there is not much docs to get in easy and that is assuming you are one of the few people who feels at easy in writing python types in C from scratch...

• The use cases for downstream authors to use this have a lot of overlap with the EA use cases. What part of the ExtensionDtype API could be handled by the new numpy implementation you're proposing (backward compat and versioning notwithstanding...)

While I would somewhat prefer to limit at least modification of the metaclass for starters. The plan is to give users full flexibility to define their DType (classes) and instances as they like.

Which would mean that even if NumPy dtypes are missing some fields of your ExtensionDTypes, you could still make them full fledged NumPy dtypes. I.e. in this case, probably mainly tag on construct_array_type and na_value, etc.

The extension dtype class will be a DTypeMeta instance:

class ExtensionDtype(..., metaclass=np.DTypeMeta):
     pass

arr = np.array([1, 2, 3, 4], dtype=ExtensionDType())

And the instance of that will be a dtype instance that can be attached to a numpy array. I do not really think that either of that interfere with pandas at all. There may be a bit of a funny thing if you map to self as the numpy array dtype.

One main difference may be that I currently did not have any aspiration for a construct_from_string method. I do not quite love it personally. You could still keep that around in Pandas though, no matter what we decide. (One main issue I still have is how to handle the .npy format saving, that would work straight forward with this construct_from_string hook).

Can you point me to this? Casting is tricky with flexible dtypes, but even if it is not perfect, I am currently fairly happy about the state. Although, I am considering adding fast-path methods for the can_cast check (i.e. users do not have to implement them, but can).

The function I was thinking of is maybe_promote which does upcasting slightly differently from numpy (Note: the function is a mess and subject of several current PRS).

can_cast checks are probably most analogous to Block._can_hold_element. ATM ExtensionArray/ExtensionDtype don't handle this directly, but probably should.

The other thing that comes to mind is lib.infer_dtype, though that is less about casting than inference.

I think we've been pretty happy with the dtype registry pattern.

This is a global type -> dtype mapping for coercion, or the python code to create a new extensiondtype/array? Just to be sure I got a clear idea, could you point me to the documentation for that?

This is a global registry so that pandas is aware of any user-defined ExtensionDTypes, e.g. registering IntervalDtype. (This may not actually be necessary now that I think of it. Going to open an issue about that...)

ExtensionArray (EA) does not have a well-defined notion of a scalar type, so there is no type->dtype mapping. ExtensionDtype does have a _construct_array_type method to get the corresponding ExtensionArray subclass.

The docs for EA authors are here. Definitely let me know if things are unclear.

I currently did not have any aspiration for a construct_from_string method. I do not quite love it personally

Yah, I'm not a big fan of that one either. So it goes.

I have to spend a bit more time with the pandas extension types maybe. Ijust tried to try out pint-pandas, but it seems not to work using pandas 0.25.1 and pint 0.9?

The thing I would like to see is simply how pint currently handles that Series * Series gets the correct new unit, but I have difficulty finding the source in pint (or the obvious hook in pandas)?

I suppose there should also be a method to achieve Series(..., dtype="pint[ft]").astype("pint[m]"), which I feel I should see (and probably is part of maybe_promote?

One main difference, as far as I understand it, is that I want to make sure everything is defined as slots on the dtype and completely separate from the container, which is something pandas does not currently do/achieve? However, that does not mean that pandas cannot keep supporting it, I think.

There is another thing in that it would be awesome to tag on dtype specific (u)functions to arrays. It seems pandas has (weird?!) hooks to do so. But I think that is again a whole new story/NEP.

Ijust tried to try out pint-pandas, but it seems not to work using pandas 0.25.1 and pint 0.9?

Sorry if I was unclear, AFAICT pint-pandas is still under construction. I've been following it intermittently in the hopes that it would come together. I may end up implementing a units EA for pandas myself, but would rather wait and avoid re-implementing lots of pint-based stuff that they have already thought through.

The thing I would like to see is simply how pint currently handles that Series * Series gets the correct new unit, but I have difficulty finding the source in pint (or the obvious hook in pandas)?

If I were to write this, I would do something like:

class UnitsDtype(ExtensionDtype):
    def op(self, other_dtype, binary_op):
        # binary_op here is e.g operator.add, operator.mul, ...
        # this is going to validate that this operation is allowed (e.g. not adding kilos+meters)
        # and return the appropriate result dtype, e.g. if dividing meters by seconds, we need a dtype for m/s
        ...

from pandas.core.arrays import ExtensionArray, ExtensionOpsMixin

class UnitArray(ExtensionArray, ExtensionOpsMixin):
    def __init__(self, data, dtype):
         self._ndarray = data   # <-- should be numeric dtype
         self.dtype = dtype  # <-- assume already validated etc to UnitDtype

    @classmethod
    def _create_arithmetic_method(cls, op):
        def method(self, other):
             other_dtype = getattr(other, "dtype", None)  # <-- probably needs to be smarter for e.g. datetime objects
             new_dtype = self.dtype.op(other_dtype, op)
             new_data = op(self._ndarray, other)  # <-- probably need to unwrap other if it is also a UnitArray
             return type(self)(new_data, new_dtype)

        method.__name__ = op.__name__
        return method

    _create_comparison_method = _create_arithmetic_method

UnitArray._add_arithmetic_ops()
UnitArray._add_comparison_ops()

(I'm not a huge fan of the _add_arithmetic_ops/_add_comparison_ops pattern, but so it goes). You can see something like this happening in pandas.core.arrays.integer.IntegerArray.

I haven't looked at all, but assume that pint implements something akin to UnitsDtype.op here.

I actually think a lot of EA subclasses (not just units) can be defined with a similar pattern. A large majority of extant EA subclasses are thin wrappers around np.ndarrays.

I want to make sure everything is defined as slots on the dtype and completely separate from the container

slots as in __slots__? If so, i'm missing something and don't see how that factors in. Or do you just mean you want the action to take place on the ExtensionDtype and not on the EA? I take it the example above dispatching to self._ndarray wouldn't work for what you have in mind?

There is another thing in that it would be awesome to tag on dtype specific (u)functions to arrays. It seems pandas has (weird?!) hooks to do so

Can you expand on this? Some of our EA subclasses implement __array_something__, but I forget which versions implement what.

Thanks, that is very helpful. I see the existence of op is interesting. So the way it works is a lot like __array_ufunc__ so to say, I think. In that you give a hook to handle all of the operand (in numpy world the ufunc). And you have to register your new functions with the ufunc rather than the object. Basically, I think what the extension array does is already somewhat possible nowadays in NumPy, and is what astropy uses for its Quantities (i.e. using __array_ufunc__), except that you cannot call it a "data type", and thus cannot use astype or our coercion functions.

Yes, everything should be on defined either:

• On the (Extension)DType
• On the universal function

Dispatching to the "_ndarray" method will work, because you will be able to say:

float64_loop = np.multiply.resolve((np.Float64, np.Float64, np.Float64))

def adapt_unit_func((in_unit1, in_unit2, out_unit), casting="safe"):
    # Do magic to find the correct out_unit, or check if its OK.

unit_loop = WrappedUFuncLoop(float64_loop,
    signature=(Unit64, Unit64, Unit64),
    adapt_dtypes=adapt_unit_dtypes))

np.multiply.register(unit_loop)

(I did not think too much about the syntax, this is not how it is supposed to look like in the end).

You would also be able to register resolver functions, since your Unit DType will want to see what numpy makes of float32 + float64 to decide what it should make of unit32 + unit64.

Can you expand on this? Some of our EA subclasses implement array_something, but I forget which versions implement what.

Oh, I just meant arr = np.arra([1, 2, 3, 4], dtype=Unit[m]) and then be able to do arr.to_si()...

And you have to register your new functions with the ufunc rather than the object

Just to clarify on the "have to" part of this: the example I gave used ExtensionOpsMixin, which is what has the _create_arithmetic_method classmethod. That is useful primarily in the (very common) case where the arithmetic method in question is a template-able wrapper around an np.ndarray method. For DatetimeArray/TimedeltaArray/PeriodArray, we define the arithmetic methods directly in pandas.core.arrays.datetimelike.DatetimelikeArrayMixin.

So to be clear: if you wanted to write an EA and use some completely different pattern to define these methods, none of that would violate the EA interface.

unit_loop = WrappedUFuncLoop(float64_loop

This looks vaguely familiar since I'm looked through the numpy code a few times for datetime-parsing-related stuff, but I've never really gotten a grasp on it. Is there some bite-sized piece of the numpy code that I could/should read to be able to discuss this more competently?

Basically, I think what the extension array does is already somewhat possible nowadays in NumPy

Now I'm curious about how numpy defines arithmetic/comparison operations. Consider DatetimeArray (DTA). In principal, op(DTA, scalar) should behave just like a vectorized version of op(Timestamp, scalar). But actually dispatching to the Timestamp implementation would be incredibly non-performant. Am I understanding correctly that op(ndarray[datetime64], scalar) is doing something akin to dispatching to the np.datetime64 implementation?

and thus cannot use astype or our coercion functions

For calling ndarray[custom_dtype].astype(np.int64) it makes sense that the custom_dtype author would define that. Are you saying that they would also define (how? where?) how to do ndarray[np.int64].astype(custom_dtype)?

And you have to register your new functions with the ufunc rather than the object

Just to clarify on the "have to" part of this: the example I gave used ExtensionOpsMixin, which is what has the _create_arithmetic_method classmethod. That is useful primarily in the (very common) case where the arithmetic method in question is a template-able wrapper around an np.ndarray method. For DatetimeArray/TimedeltaArray/PeriodArray, we define the arithmetic methods directly in pandas.core.arrays.datetimelike.DatetimelikeArrayMixin.

So to be clear: if you wanted to write an EA and use some completely different pattern to define these methods, none of that would violate the EA interface.

Yes, but I think it is all pretty similar in our world. Your interface is very similar to __array_ufunc__, I think. __array_ufunc__ can also override anything. the difference with new dtypes in NumPy will be that the container will be a true ndarray. Unlike what astropy.units or pandas ExtensionArrays does.

unit_loop = WrappedUFuncLoop(float64_loop

This looks vaguely familiar since I'm looked through the numpy code a few times for datetime-parsing-related stuff, but I've never really gotten a grasp on it. Is there some bite-sized piece of the numpy code that I could/should read to be able to discuss this more competently?

This pattern (or anything like this) does not really exist right now. What we have is __array_ufunc__ which is a bit like your extension array op dispatching (although I do not know how close it is, with the dispatching rules and all).
Inside __array_function__ you could do input.view(...) and then call the original ufunc again.

Basically, I think what the extension array does is already somewhat possible nowadays in NumPy

Now I'm curious about how numpy defines arithmetic/comparison operations. Consider DatetimeArray (DTA). In principal, op(DTA, scalar) should behave just like a vectorized version of op(Timestamp, scalar). But actually dispatching to the Timestamp implementation would be incredibly non-performant. Am I understanding correctly that op(ndarray[datetime64], scalar) is doing something akin to dispatching to the np.datetime64 implementation?

Those are all ufuncs, and containers have a the __array_ufunc__ hook, which allows them to intercept the call to np.equal(obj1, obj2) if they so wish. They can do whatever they want inside that call. I suppose you may have to define the normal __eq__ operator as well.
In this sense, DTA has to define the vectorized version itself of course (although it can fall back to using other numpy functions where possible).
That is a bit different for the new dtypes, which would allow you to define a fast inner-loop (possibly in C), and does all the broadcasting etc. for you.

and thus cannot use astype or our coercion functions

For calling ndarray[custom_dtype].astype(np.int64) it makes sense that the custom_dtype author would define that. Are you saying that they would also define (how? where?) how to do ndarray[np.int64].astype(custom_dtype)?

Of course the custom_dtype author has to define both directions. I am currently opting for an operator approach __can_cast_from_other__ and __can_cast_to_other__ (names and some details up for change of course). The other possibility would be to do np.cast.register(...), but I think I prefer the slot/dunder method approach, and it actually is subtlety different from other ufuncs which would have that registration approach.