---
title: New Datatype dependend UFunc Dispatching
---

[todo/meeting notes](https://hackmd.io/5aV5K49pT8GtxBcMgC7PYQ)

# NEP XX — New Datatype dependent UFunc Dispatching

:Author: Sebastian Berg <[e-mail]>
:Author: ...
:Status: Draft
:Type: Standards Track
:Created: 2019-07-17

Abstract
--------

We propose revising the universal function (ufunc) dispatching to better
support user defined data types (dtype) as well as simplify the implementation
of ufuncs for example for the numpy internal string dtypes.
In the new model all calculation logic is moved into `UFuncImpl` objects
which are specific to the input and output dtypes.
An additional dispatching step will be used to decide what the correct
`UFuncImpl` is for the given call.

Motivation and Scope
--------------------

The current dispatching system has a setup step which can for example
decides which output dtype to use, and which is used to find the correct
dtype signature for the calculation.
However, this system is not extensible by extension dtypes [1] and due
to the necessary complexity has deterred the implementation of ufuncs for
strings (flexible dtypes) for which the length of the output string has to be
found.

Furthermore, the new system should allow for optimization.
The current inner loops (the functions handling the low level calculations)
often include multiple branches to optimize execution and checking whether
or not an error occured is not as specific as it could be.

This NEP whishes to address all of this by making the calculation and ufunc
execution dtype specific.
It does not suggest to implement any behaviour change with respect to the
current data types.

Detailed description
--------------------

The situation will be described in terms of the simple ufunc call;
similar things apply for all other ufunc methods (e.g. ``reduce``, ``outer``).

The current implementation can be described by the following simplified
pseudo code:

class UFunc:
    def __call__(input_arrays, **kwargs):
        self.handle_array_ufunc_dispatching(input_arrays, **kwargs)

        dtypes = tuple(arr.dtype for arr in input_arrays)
        dtypes = dtypes + (None,)  # can include output dtypes

        # Find the exact dtypes used:
        loop_dtypes = self._resolve_types(dtypes)

        # Find the matching loop which handles the calculation:
        inner_loop = self._find_inner_loop(loop_dtypes)

        self.__execute_loop(input_arrays, loop_dtypes, inner_loop)

    def __execute_loop(input_arrs, loop_dtypes, inner_loop):
        """
        Handle the outer iteration and call the inner loop, uses
        `np.nditer` to allocate and cast input arrays.
        """
        iter = np.nditer(input_arrays + (None,),
                         op_dtypes=loop_dtypes, casting="same_kind",
                         flags=["external_loop", ...])
        for arrs in :
            inputs = arrs[:-1]
            output = arrs[-1]

    def _resolve_types(dtypes):
        """
        Find the correct dtypes for all operands (including the output)
        and will often raise an error if there is no possible inner loop.

        This slot is specialized for many ufuncs, for example datetimes require
        that this step checks and fixes the datetime unit (ns, ms, ...).
        For simple math functions, this will usually return
        ``np.result_type(*(dtypes + out_dtypes))``.
        """
        return loop_dtypes

    def _find_inner_loop(self, loop_dtypes):
        """
        Find the correct inner loop. This is often an exact match to all
        dtypes, but can also be any loop for which all casts can be done
        safely.
        """
        return correct_inner_loop

    def register_inner_loop(self, dtypes, loop):
        """
        It is possible to register new inner loops for specific dtypes.
        """
        pass

which allows ufunc specific behaviour only by modifying the ``_resolve_types``
and ``_find_inner_loop`` slots (which is used internal).
New inner loops can be registered, and inspected using for example
``np.add.types``.
Although it is slower and limited, the above does accept user defined dtypes.

This NEP suggests to change this into the following:

class UFunc:
    def __call__(input_arrays, **kwargs):
        self.handle_array_ufunc_dispatching(input_arrays, **kwargs)

        dtypes = tuple(arr.dtype for arr in input_arrays)
        dtypes = dtypes + (None,)  # can include output dtypes

        ufunc_impl = self.__resolve_ufunc_impl(dtypes)
        
        return ufunc_impl(input_arrs)

    @memoize  # This function should be memoizable (but may not need to be)
    def resolve_ufunc_impl(dtypes):
        resolver = find_best_registered_resolver(dtypes)
        return resolver(dtypes)

    def register_resolver(dtypes, resolver):
        """
        Register a new resolver function for the specified dtypes.
        
        Parameters
        ----------
        dtypes : tuple of dtypes
        resolver : callable
            A function which returns a new `UfuncImpl` or sets an error.
        """
        pass


class UFuncImpl:
    dtypes = (Float64, Float64, Float64)  # specific dtype classes

    def __call__(input_arrays):
        pass


class LegacyUfuncImpl(UFuncImpl):
    def __init__(self, dtypes, inner_loop, identity=None):
        self.dtypes = dtypes
        self.inner_loop = inner_loop
        self.identity = identity

    def __call__(self, input_arrays):
        return self.__execute_loop(input_arrs, loop_dtypes, inner_loop):

    __execute_loop = Ufunc.__execute_loop  # Use largely same implementation.

Here ``LegacyUfuncImpl`` would be as close as possible to the current ufunc
with the exception of only implementing a single type signature.
New subclasses of ``UFuncImpl`` could change the way ``__execute_loop`` works.
This also means that the C defined inner loops can work differently.
Specifically, we suggest to implement a new ``InnerLoopUfuncImpl`` which
implements works the same as the legacy one but expands functinality with
respect to:
* including an error indicator return value in the inner loop allowing for
early returns.
* allowing for specialized setup/teardown functions to allow to:
- optimize error checking (integers do not need to check floating point
errors)
- easier pass data into the inner loop (e.g. string dtypes need
to pass string length information)
- allocate and clean up working memory if the ufunc requires it
- possibly specialize the inner loop function.

Since in principle a ``UFuncImpl`` can replace the complete logic of the
ufuncs, it is not necessary to cover all use cases in the first implementation.
New ``UFuncImpl`` subclasses can easily be added later on.

#### Discussion:

To begin with, we will limit the accepted subclasses of ``UfuncImpl`` to
allow us to gain experience before settling for a fully public API which would allow users to do arbitrarily complex things.

### ``UfuncImpl`` Resolution

The above text is intentionally unclear about how the exact resolution and
finding of the best resolver function will be implemented.
This mainly depends on the decision for the dtype hierarchy and can easily
be extended in followup steps.

There are a few possible solutions with regards to the dtype resolution.
The biggest difficulty is that we need to handle value based promotion,
this means that for some dtype inputs a simple mapping from dtype classes
of the input (and provided dtypes) to the output dtypes is not generally
possible. This will require dynamically created dtypes, which
will be hidden from the user (**TODO:** It may be that these need to be
dtype classes, but maybe dtype instances are good enough, which would be nicer).

Dispatching has two modes, which may overlap:

* An exact matching type signature, when all types are given.
* A type resolution function, which matches also when not all types are given;
if more then a no single resolution function clearly matches best.
This will support an abstract type hierarchy, where dtypes can fall into
at least one category, such as `NumericalDType`.
(*TODO:* Depending on the type hierarchy and complexity we allow, this
might been a single matching resolution function). One thing to keep
in mind is that we have ufuncs that have both `OO->O` and `OO->?` defined.
- Type resolution may fall back (as it does now) to a common type operation
and repeat the resolution step with the fully specified signature.
- Type resolution may also return a new ``UfuncImpl`` object; it must return
an identical one for the same input.

### Usage

**TODO:** Add an example of a ``UFuncImplWrapper`` implementation which
can be used for example for Units, or other dtypes which can reuse existing
implementations.

Related Work
------------

**TODO:** This section should list relevant and/or similar technologies, possibly in other
libraries. It does not need to be comprehensive, just list the major examples of
prior and relevant art.

Implementation
--------------

**TODO:** This section lists the major steps required to implement the NEP. Where
possible, it should be noted where one step is dependent on another, and which
steps may be optionally omitted. Where it makes sense, each step should
include a link to related pull requests as the implementation progresses.

Any pull requests or development branches containing work on this NEP should
be linked to from here. (A NEP does not need to be implemented in a single
pull request if it makes sense to implement it in discrete phases).

Backward compatibility
----------------------

If necessary breaking downstream libraries which directly modify
the current UFunc struct may be in scope. To the best of our knowledge
this should only affect the astropy and numba projects.
Generally, the proposed solution implements a superset of features based
on the old proposal.
All API functions are planned to be supported, although they may be
deprecated.

It may turn out to be easier to fix all (public) downstream uses rather
than ensuring backward compatibility.
If this is the case, and the usage does not affect a huge audience,
priority will be given to moving the project forward rather than ensuring
compatibility at all cost.

Alternatives
------------

**TODO: should discuss alternatives.**

Discussion
----------

**TODO:** This section may just be a bullet list including links to any discussions
regarding the NEP:

- This includes links to mailing list threads or relevant GitHub issues.

References and Footnotes
------------------------

.. [1] Downstream users could modify the ``TypeResolution`` slot directly in
principle, but this is neither used nor a reasonable API. Such modification
is used by Astropy, however, only for its own ufuncs.

.. _Open Publication License: https://www.opencontent.org/openpub/

Copyright
---------

This document is publish under the Open Publication License [_Open].
(**TODO:** maybe public domain is OK/works?)