Roadmap: Next milestone for numpy-quaddtype (post v1.0.0)

[SwayamInSync]

For the next iteration we can focus on upgrading the dtype API at NumPy and other performance related issues. Listing down some known ones to track progress.

Dtype API upgrades (requires changes at NumPy)

• string parsing of user-defined dtypes
  • Tracking Issues: ENH: add a hook to the DType API for parsing dtype strings numpy#30699
  • PR:
• Exposing hook for np.dot & fixing np.einsum support for user-defined dtypes
  • Tracking Issue: ENH: Make functions in dot family and related to support user-defined data types numpy#30793
  • PR:
• Paths depending on copyswap or copyswapn leading to segfault
  • Tracking Issue: segfault: QuadPrecision(1).byteswap() #13 , np.place segfaults for quaddtype #7 , ENH: dtype API needs an API slot for a single-element copy numpy#24859
  • PR:
• Using string form in array printing for user-defined dtypes
  • Tracking Issue: QuadPrec scalar repr isn't ideal for display in arrays #29 , BUG: Array printing bypasses user dtype's formatting silently corrupts precision. numpy#30796
  • PR:
  • Extra context: Subclassing dtype to the native dtype (like quaddtype < np.floating) seems to give good formatting but not upto the entire precision, that has to be explicitly set via np.printoptions(...) which might not be ideal as we want the printed output to be in native precision (or any other properties) of the dtype by default
• Allowing np.linalg submodule to allow user-defined dtypes routine registration
  • Tracking Issue:
  • PR:

Updates/Fixes in numpy-quaddtype

• Supporting WASM build
  • PR: BUILD: Add Pyodide CI and build recipes #66
• casting to quaddtype emits uneccesary warnings
  • Tracking Issues: casting to quaddtype emits uneccesary warnings #25
  • PR:

[SwayamInSync]

I think we can create separate tracking issues for the numpy API fixes in numpy repo and initiate some forward discussion on each and the potential fixes, so that contributors from both numpy and user-dtypes side can come forward as per the interests.

cc: @ngoldbaum @seberg

[ngoldbaum]

Please do open numpy issues for improvements for the user dtype system.

[mdhaber]

Thanks for working on this! I'm excited about the potential for use in SciPy, especially the elementwise iterative methods like scipy.differentiate.derivative and scipy.integrate.tanhsinh (quadrature), which are written in terms of the array API. How feasible is support throughout NumPy's array API functions? I thought differentiate.derivative might "just work", but I ran into trouble with input validation - np.isdtype(..., "real floating") didn't recognize a QuadPrecDType yet. If it would be helpful, I can try to find pain points in some representative SciPy features.

[ngoldbaum]

If it would be helpful, I can try to find pain points in some representative SciPy features.

That would be great, if only to figure out where to focus effort in NumPy.

I think there are open issues over there about better fleshing out DType hierarchies and adding ways to register dtypes in a hierarchy, but right now there isn't anything like that.

[SwayamInSync]

Thanks @mdhaber for looking into it and definitely it'll be great to get know of more limitations.

but I ran into trouble with input validation - np.isdtype(..., "real floating") didn't recognize a QuadPrecDType yet.

I checked isdtype and the docstring stays "This function only supports built-in NumPy's data types". One recommendation might be here to use issubdtype

In [5]: np.issubdtype(QuadPrecDType(), np.floating)
Out[5]: True

but I feel we should also update isdtype function as well inside NumPy

[mdhaber]

One recommendation might be here to use Issubdtype

Ah, but that is not part of the array API. I could add a special case temporarily as a workaround, but it's best to minimize that.

[mdhaber]

It was not too hard to get scipy.differentiate.derivative and scipy.integrate.tanhsinh working with a QuadPrecDType; the only real incompatibility was xp.isdtype in a few places (as discussed above). A few other things I noticed:

• A quad precision linear solve would be very helpful
• It looks like str of a quad dtype does not (at least by default) use scientific notation? (e.g. print(np.finfo(QuadPrecDType()).smallest_normal) prints all ~5000 digits, mostly leading zeros.)
• I'm seeing some warnings from NumPy for valid arithmetic (e.g. np.log(np.asarray(1e-300, dtype=QuadPrecDType())**2))

Otherwise, I think it will mostly "just work" with these and other features that could really benefit from very high working precision. Thanks!

[SwayamInSync]

Thanks @mdhaber . This might be an important issue, I just checked and np.linalg module only supports hardcoded _real_types and _complex_types and directly calls the LAPACK routines via _umath_linalg gufuncs.

This is inconsistent with how ufunc works, we need to make _umath_linalg gufuncs support the same loop registration for dtypes that regular ufuncs do. This might require big refactoring of gufunc infrastructure (but should be good plan for long term)

[mdhaber]

Yeah, I imagine linalg support is a tall order, but thought solve was worth mentioning. If we ever get a special function extension into the array API standard (data-apis/array-api#725), that would be the next thing I'd ask for.
@steppi any chance that xsf could have quad precision implementations of some functions?

[SwayamInSync]

I think linalg support is worth solving in NumPy as well, I will by this weekend raise all the found issues over NumPy and we can start working towards them.

@steppi any chance that xsf could have quad precision implementations of some functions?

The external implementations of routines might not work with our dtype as to ensure cross-compatibility the quaddtype is defined differently over different platforms but this won't be a pain point as its straightforward to implement such operations within the package and then register them in NumPy to use.

[steppi]

If we ever get a special function extension into the array API standard (data-apis/array-api#725), that would be the next thing I'd ask for.
@steppi any chance that xsf could have quad precision implementations of some functions?

From a numerical perspective, I know that the main thing we'd need to do that would take more than minimal work is to fit some quad specific mini-max polynomial approximations. We could use Boost's or Cephes' implementation of Remez-exchange for this. I don't know enough about numpy-quaddtype yet to know how defining ufuncs would work for the quad dtype, but expanding dtype support in xsf by fitting more minimax approximations is planned.