RFC: Adoption of ``std::mdspan`` like structures in C++ code about scipy HOT 40 OPEN

+1, but I'm afraid compiler/ecosystem support remains a blocker for the foreseeable future?

Very good point, though we are also in the tough spot for switching to anything newer than C++17 but do you have particular architectures in mind?

I regularly come back to wanting to have a NEP29-style document for our platform/arch/compiler support, perhaps like the tiered approach in Cpython or rust.

In any case, if Rust is good enough for the kernel, I have a hard time imagining that it couldn't be used in SciPy; C support is certainly broader, but are there actually any architectures these days with a C++ compiler that's unavailable for rust?

Aside from the unclear outcome with the various standardisation efforts, we're also signing up for a long wait before we can require C++20 (much less C++23/26), so while I'm sympathetic to the overall effort being proposed here, I do see the argument by @ilayn too, and tend to prefer we'd go for rust.

from scipy.

Also, @fancidev see #19601 and cupy/cupy#8140 for examples of this in action. I’d say we can’t use Boost directly because it doesn’t support CUDA. The hope is to have a single special function library that could be used by SciPy, CuPy, PyTorch, Jax, and other array libraries.

from scipy.

+1, but I'm afraid compiler/ecosystem support remains a blocker for the foreseeable future?

Very good point, though we are also in the tough spot for switching to anything newer than C++17 but do you have particular architectures in mind? I'm not as versed as you are regarding different architectures but I tried already with a few embedded targets with preliminary success.

I had a great experience improving perf with Rust at #14719.

Oh how did I miss that? That's a wonderful result and also gives me more confidence that I am not the only one getting affirmative results on the "Fearless Concurrency" claims

from scipy.

If what @steppi said about making special a separate codebase and hydrating all stakeholders from a common place, is going to happen soon, then noone probably will care what you did in the bowels of the library and just compile it. But switching to mdspan overall in SciPy seems to me a really bad investment regardless what the standard committee claims.

I think it depends to an extent on whether our EOSS 6 proposal is funded, otherwise it might be more difficult to secure the bandwidth, but I'd hope to put it in a separate codebase within one year, and and start doing outreach to get more stakeholders interested within 6 months.

At the moment I don't feel qualified to have an opinion on the questions of Rust vs future versions of C++ vs <???> and how exactly we should standardize on ways to work with arrays in low level code, but plan to take some time to educate myself on these matters.

from scipy.

As an outsider to some aspects of the discussion - in many of my own projects I let cython deal with all the hard work of multidimensional array wrangling (strides, etc). After all, it's a superset of Python so the syntax is pretty easy. Then I drop down into C/C++ for the actual scalar kernel computation once all the hard work has been done.

From past experience the more is in Python/cython (and less in C) the easier it is to maintain. e.g. With the TNC optimiser code the extension interface was originally written in C, calling into vendored C code. IMO writing Python extensions directly in C is painful, lots of weird function names (e.g. PyArg_ParseTuple) to understand, etc. This made it hard to maintain. Rewriting the interface in cython made it much more maintainable and easier to extend in the future.

from scipy.

First and foremost, IMO this RFC is not only for C++ experts! It may impact us mere mortals as well.

With the C++ momentum by @izaid and @steppi in scipy.special and other people, I think it'd be great to

• formulate a best-practices guide for C++ in SciPy, and
• put common building blocks to _lib or _build_utils for all submodules to use.

At least three things come to mind:

• how to link LAPACK from C++ (this is gh-20002)
• how to properly pass complex-valued arrays from python to C++ and back (scipy.special developed some means? but it's awkward to use in other subpackages)
• how to grab a view onto a numpy array from C++.

For the latter, std::mdspan seems to be the recommended way. So my suggestion would be to put the backcompat implementation from gh-20320 to _build_utils or some such. One other thing from past experience: we need a way to toggle bounds checking at compile time, this saves long hours of debugging.
So an access_policy with bounds checking would be a great addition to what is being introduced in gh-20230.

One other relatively minor sticky point in C++ vs python glue is signed / unsigned indices: C++ containers use size_t while python indexing is ssize_t (I believe?). A best-practices recommendation of how to convert from one to the other would be great.

from scipy.

Hi all! As has been mentioned, @steppi and I are on a long-standing quest to greatly improve SciPy's special. Our goal is to get all of the numerical functions in special currently mixed between C, C++, Cython, Python, and F̶o̶r̶t̶r̶a̶n̶ into a single header-only C++ library. That library will then be able to provide special function implementations to all the other array libraries, and indeed it already is being used in CuPy. We are doing this by converting the special function implementations into NumPy ufuncs with C++ kernels. There are some tools we need to do this (and those will exit SciPy when this header-only library becomes its own submodule), so another question is what SciPy wants in general.

As part of our efforts, we've come across all the things discussed above. I'm happy to help if useful. Some thoughts are below.

For complex numbers, we use std::complex<float> and std::complex<double>. You can directly cast a pointer to a npy_cfloat and a npy_cdouble to them, then operate on them in a convenient fashion.

The other thing that came up is how to handle functions that operate on subarrays. In our case, this arose from functions like the associated Legendre polynomials that return a (m + 1, n + 1) array for each value z in the input. Really all that is needed is a simple class that wraps a pointer, a shape, and strides, no memory ownership and no reference counting. One can either write their own, or use the upcoming standard std::mdspan from C++23 that is exactly intended for this use case. It was inspired from the idea of NumPy views, but without the reference counts. CuPy is an example of a library that wrote their own, the CArray above, but they did that before std::mdspan was accepted into the standard.

Note that to do views without memory allocation, the number of dimensions must be known at compile-time. I don't think this is a great restriction, as most functions in SciPy that need to drop to C++ somehow know if they are operating on 1D, 2D, 3D, or whatever arrays.

Most of the special work has been tools to help with ufuncs, I think what's needed here is actually way simpler. Again, happy to help if that is useful.

from scipy.

For complex numbers, we use std::complex and std::complex. You can directly cast a pointer to a npy_cfloat and a npy_cdouble to them, then operate on them in a convenient fashion.

I would be really happy if that works out of the box but I am highly suspicious of it to be honest. I just need to see to believe it.

The other thing that came up is how to handle functions that operate on subarrays. In our case, this arose from functions like the associated Legendre polynomials that return a (m + 1, n + 1) array for each value z in the input. Really all that is needed is a simple class that wraps a pointer, a shape, and strides, no memory ownership and no reference counting.

That's the sales pitch but not that the arrays must be C-order contiguous or strided etc hence we need to babysit everything just as we do for regular C flatbuffers. Otherwise you have to copy things to fresh arrays. So this problem is not solved by C++, this is already the default behavior. All it brings is the md part over the regular span. So it will be as annoying as Cython memoryviews or malloc'ed arrays to convert back and forth between NumPy arrays.

Hence my not so enthusiastic comments about this new class. It is accepted to standard but still miles away from a NumPy or Fortran like convenience and indexing/slicing support and safely manipulating them. To be honest, I don't even know what the benefit is for us compared to regular std::span.

Also who writes these pages I will never understand; https://en.cppreference.com/w/cpp/container/mdspan

from scipy.

That's the sales pitch but not that the arrays must be C-order contiguous or strided etc hence we need to babysit everything just as we do for regular C flatbuffers. Otherwise you have to copy things to fresh arrays. So this problem is not solved by C++, this is already the default behavior. All it brings is the md part over the regular span. So it will be as annoying as Cython memoryviews or malloc'ed arrays to convert back and forth between NumPy arrays.

Hence my not so enthusiastic comments about this new class. It is accepted to standard but still miles away from a NumPy or Fortran like convenience and indexing/slicing support and safely manipulating them. To be honest, I don't even know what the benefit is for us compared to regular std::span.

Sure, that's definitely why it's good to have a discussion! I'm here mostly to provide information about what we've done in special.

I will note that std:mdspan is meant to support different layouts, and it already has compile-time customization for row-major (rightmost contiguous), column-major (leftmost contiguous), and arbitrary strides. If you want something more exotic, it can also be customized. It's also been backported to CUDA by NVIDIA.

It's primary purpose is to establish a standard pattern for views on multidimensional arrays in C++. One can always create their own thing instead.

from scipy.

Yes sure, I'm not holding you accountable or ambassador of C++ 😃 It is just what it is.

from scipy.

Yes sure, I'm not holding you accountable or ambassador of C++ 😃 It is just what it is.

Personally, I do think it's the right way to go. These things are helping a lot in special, so perhaps they can also help elsewhere.

from scipy.

In case it is relevant, I know the National Labs here implemented https://github.com/kokkos/mdspan, and the PyKokkos Python bindings to Kokkos are still being developed at UT Austin (I was helping them implement ufuncs before I ran out of bandwith).

from scipy.

In case it is relevant, I know the National Labs here implemented https://github.com/kokkos/mdspan, and the PyKokkos Python bindings to Kokkos are still being developed at UT Austin (I was helping them implement ufuncs before I ran out of bandwith).

Not only is that relevant, the Kokkos repository is actually what we're using in special!

from scipy.

Kokko's/CUDA's implementations seem worth looking at, does anyone know whether their licenses are compatible with SciPy's?

from scipy.

Our goal is to get all of the numerical functions in special currently mixed between C, C++, Cython, Python, and F̶o̶r̶t̶r̶a̶n̶ into a single header-only C++ library.

May I ask how is this going to be related to those special math functions already available in boost?

from scipy.

Our goal is to get all of the numerical functions in special currently mixed between C, C++, Cython, Python, and F̶o̶r̶t̶r̶a̶n̶ into a single header-only C++ library.

May I ask how is this going to be related to those special math functions already available in boost?

We'll use Boost where we can, but Boost doesn't support CUDA among other things. SciPy has an enormous collection of special functions, they are valuable in and of themselves. We discussed this several months ago, see #19404

from scipy.

Thanks both for the reply and references! I find this comment quite informative and answers all my questions. Looking forward to have those functions reusable in C++ (even if by copy pasting)!

from scipy.

Personally, I do think it's the right way to go. These things are helping a lot in special, so perhaps they can also help elsewhere.

I have nothing against C++ as it is being used in many parts, I would not use for anything but that's anectodal.

But if I consider numerical work, I am not really convinced about C++/F90 in terms of portability. So far there are two, that stood the test of time; one is C and the other is F77. Both C++ and F90 are standing on the laurels of their "less well-mannered cousins". The rest is just lots of some successful and some broken codebases with mixed feelings. I know Sandia labs went all in with C++ which is perfectly fine for them with their team fully behind that decision. A lot of others went with Modern Fortran, whatever that means today. But we are not a C++ codebase. We just need to get by with maintainable code + performance if possible. Much to my astonishment after so many decades, Reference LAPACK is trying to get some C++ code into the majestic Fortran codebase that we all love with their home-made matrix classes Reference-LAPACK/lapack#991 So nobody is agreeing on a proper array structure and we are in 2024. So this makes no sense to me and has all the signs of sunken cost for the future.

For linear algebra work, currently mdspan does not really offer anything yet if compared to, say even in C++ ecosystem, Eigen, and that's not smooth sailing at all this page being my favorite.

For special it might be OK since there is not that many array operation requirements but if you want to do a bit more tricky ops then mdspan is just an unnecessary abstraction that will continuously get in the way.

However there is std::linalg is in the works (C++26 and onwards) and who knows when we will be able to switch to it, if we ever stay with C++. Then things can get really nice if it works.

But if I accept that level of complexity, I would definitely switch to Rust instead which is already much better than this. So not sure where the benefit of mdspan lies. If you don't want the whole bulky Eigen dependency and still use mdspan locally inside special that's obviously a different requirement for keeping a lean header only module. I totally get that.

Like we discussed in the very early times on another RfC, C++ stuff is only OK in my opinion if we stay away from the bad parts of C++ and these data structures look like the bad parts of C++ to me with unnecessary abstractions without the actual convenience we need that Fortran or NumPy (or even Eigen) offers.

If what @steppi said about making special a separate codebase and hydrating all stakeholders from a common place, is going to happen soon, then noone probably will care what you did in the bowels of the library and just compile it. But switching to mdspan overall in SciPy seems to me a really bad investment regardless what the standard committee claims.

from scipy.

I would definitely switch to Rust instead

+1, but I'm afraid compiler/ecosystem support remains a blocker for the foreseeable future? I had a great experience improving perf with Rust at #14719. I'm with you on the pain of C++, to me it is like the saying about regular expressions where you then have two problems, the actual problem and dealing with C++.

That said, if folks are willing to maintain it I guess I've mostly stayed out of the way. The C++ adoption seems to be happening in NumPy too. We used to have a rule about order of preference being this for maintenance reasons: Python, Cython, C, C++, Fortran. I believe that was in the SciPy 1.0 paper, but then there were various discussions around allowing C++ more frequently.

from scipy.

I think the important question now is then: would converting from C++ (using similar features to the current work) to Rust be significantly more work than converting directly from C / Cython?

from scipy.

Every rewrite takes a substantial effort, so converting first to C++ and later to Rust certainly needs more effort than one rewrite.

That said, the one argument I see for the C++ side is GPU support, not least since Nvidia is strongly betting on C++. I don't know what Rust's story for heterogeneous computing is, but given that it's all LLVM under the hood, I kinda expect it to be... not that bad (certainly given the time we'd wait for universal support of newer C++ standards)?

from scipy.

After reading the great discussions, it seems the main debate is whether the complexity of C++ justifies its adoption in SciPy.

It appears to me the key question is whether the C++ code to be introduced will be used exclusively by SciPy, or do you foresee it to be adopted by other projects? The latter goal is certainly much more ambitious, more complex, and more demanding, and for that goal C++ may be a suitable choice because of its wide adoption in industries.

from scipy.

or do you foresee it to be adopted by other projects

I'll let @steppi answer himself, though widespread reusability is certainly a goal AFAICT (both in terms of libraries that can leverage this, as well as architectures that this can run on, in particular GPUs).

I think no matter the implementation choice, we realistically need to provide the API in multiple flavours (either through a FFI or wrappers), e.g. C-style, C++-style, and potentially Rust-style (depending on what the ecosystem ends up doing, but certainly if the implementation itself were to be in Rust).

from scipy.

compiler/ecosystem support remains a blocker for the foreseeable future?

So if we drop the kokkos implementation somewhere into scipy, just as gh-20230 does, do we expect issues from the compile/build side?

Realistically, I don't think we could/should/would maintain it in any meaningful sense, it'd be just there as a black box. Would that cause problem?

from scipy.

It's been really interesting to watch this conversation. And I think the point about the Linux kernel adopting Rust is a good one.

I also think it's worthwhile to consider the different topics that have brought up in this thread. The original post was "can we have a standard way of handling arrays in C++ code". Then there was some discussion about "Rust is pretty good, can we use it instead". Rather than get into the points about C++ versus Rust, perhaps there should be some serious thought about how to get the best practices of each language for the whole project. It's also worth noting that Python extension modules are written in C, hence there is a natural connection to the C / C++ ecosystem.

In a very real way, it seems SciPy has always been a multilanguage project: C, C++, Fortran, Python, maybe Rust too at some point. And that also evolves with time. For various good reasons, Fortran is now being taken out of the codebase. Maybe at some point Rust will become a large part of it. Maybe SciPy will even branch off from just Python, like how IPython became Jupyter. In the meantime, I would think the point isn't choosing C++ versus Rust, but writing the best C++ or the best Rust or the best whatever.

So, back to the original question. I would think the need for a data structure that reinterprets memory as a multidimensional array is clear. That's exactly what NumPy is for Python. This data structure really only needs to do the memory reinterpretation, which basically means providing indexing for different memory layouts or making subviews. If you are working in C++, you will need something in C++. If you are working in Rust, you will need something in Rust. If we don't provide this, people are forced to mix up things like pointer indexing with the actual numerical algorithms. That has been one of the problems in special.

So what is the solution for C++? Either write your own or use what has now been standardised. The benefits of using something standardised now is it'll be right there in the upcoming standard and it can interact with other code in the ecosystem. A potentially good way to go would be to work out best practices for C++, but that wouldn't mean C++ necessarily just appears all over the ecosystem. I guess it would still have to be discussed case by case. A separate issue could discuss working out best practices for Rust and potential adoption.

from scipy.

compiler/ecosystem support remains a blocker for the foreseeable future?

So if we drop the kokkos implementation somewhere into scipy, just as gh-20230 does, do we expect issues from the compile/build side?

Realistically, I don't think we could/should/would maintain it in any meaningful sense, it'd be just there as a black box. Would that cause problem?

It is backported to C++14. It has the same issues as other third-party code, you are reliant on someone else having done a good job! To me, it seems quite good. Also to NVIDIA, who have adopted it in their CUDA standard library.

The other option is to just write your own class like the CuPy CArray mentioned earlier, in which case you are just creating the same thing (which is not necessarily bad, but that's what it is).

from scipy.

I would think the point isn't choosing C++ versus Rust, but writing the best C++ or the best Rust or the best whatever.

This is my main argument indeed that C++ after 40 years has no best practices. Just lots of people saying the others are doing it wrong and defining newer and newer things. Note that Eigen is also a standard tool de facto accepted by everyone. So is all in between implementations for span and vector and so on. The language is void of a proper array structure that's why I think we should not accept experimental features of C++ even though they are backported, compilers are the important detail not the standard.

from scipy.

Out of curiosity, is there a procedure in SciPy developers to resolve design choices if they fail to arrive at a consensus? Like a vote or something?

from scipy.

or do you foresee it to be adopted by other projects

I'll let @steppi answer himself, though widespread reusability is certainly a goal AFAICT (both in terms of libraries that can leverage this, as well as architectures that this can run on, in particular GPUs).

@fancidev . I think for special, where we just need to have scalar kernels to plug into ufuncs and gufuncs, it makes a lot of sense for the scalar kernels to be a separate library that could be shared by different projects, and it doesn't really introduce any added difficulty.

But in a lot of cases, I think it makes more sense to try to write as much as possible using Python + Array API standard, offloading performance critical parts to compiled code. Perhaps, with careful forethought, we could factor things in such a way that these performance critical parts could be shared between array libraries. That sounds very difficult though, much more difficult than in special where what should go into a compiled kernel, and how these compiled kernels are used in Python, are very clear and well defined.

from scipy.

Out of curiosity, is there a procedure in SciPy developers to resolve design choices if they fail to arrive at a consensus? Like a vote or something?

Technically we have BDFL, though he's inactive at the moment. Seems to be voting-based mostly, though usually seems to reach an informal agreement between devs before it gets as far as actually needing to count votes. This may help a bit: https://scipy.github.io/devdocs/dev/governance.html

from scipy.

A few thoughts on the actual proposal here:

• Vendoring a well-tested thing like std::mdspan (leaving aside the licensing issue I commented on in gh-20320) seems clearly preferable to me to reinventing the wheel here,
• It's a lot of code, and we don't have a great history with directly vendoring code and then keeping it up to date; we may want to stick this in a separate git submodule for traceability. If the current mdspan.h is the result of running make_single_header.py from https://github.com/kokkos/mdspan (is it?), that could even be done dynamically to avoid committing anything that's not directly from upstream.

from scipy.

Re the preferred programming languages discussion: in most cases the answer is "it depends". There are pros and cons to everything. The most clear change we've had over the last 1-2 years: new Fortran 77 code is no longer acceptable. Beyond that, there are trade-offs. C++ seems like the best of a not-so-great lot at least for special.

Cython (or Pythran) are the right tool in certain circumstances where pure Python doesn't cut it, but certainly not always. They're easier to maintain for not-too-complex code for the average SciPy maintainer. And Cython itself is fairly well-maintained. However, we also know there are significant downsides to Cython: too large binaries, long compile times, a lot of extra complexity in the build system (and poor build system support in general), fused types and const support are not very mature, niche language so little general-purpose tooling support, and the code generation for our Cython bindings and exported cython_special/blas/lapack APIs are some of the most complex and difficult to understand/modify parts of the code base we have. I'm pretty sure that the notes in http://scipy.github.io/devdocs/dev/roadmap-detailed.html#use-of-cython are still correct; Cython isn't really scalable to the size of the whole SciPy code base.

tl;dr it depends

from scipy.

Re: cython vs C++ : I think we could do with a worked example. scipy.special gives one (and an excellent one at that!). However, it's somewhat special : it seems one can use numpy (g)ufuncs + scalar kernels. What would be helpful I think is to take a cython extension, port it to C or C++, including python bindings. It would be interesting to see what the binary size would be if we go from programming in Cython to programming in C/C++ and only wrapping in Cython. Does this solve the binary size issue? I don't know and I doubt it can be just guessed.

Likewise for Rust: if we are to seriously consider it, let's first look at a concrete example.

from scipy.

Re: cython vs C++ : I think we could do with a worked example. scipy.special gives one (and an excellent one at that!). However, it's somewhat special : it seems one can use numpy (g)ufuncs + scalar kernels. What would be helpful I think is to take a cython extension, port it to C or C++, including python bindings. It would be interesting to see what the binary size would be if we go from programming in Cython to programming in C/C++ and only wrapping in Cython. Does this solve the binary size issue? I don't know and I doubt it can be just guessed.

Likewise for Rust: if we are to seriously consider it, let's first look at a concrete example.

I think this is a very good point! In the case of special, things are very well-defined. All the C++ functions are (g)ufunc kernels, and literally all the language bridging is done by (g)ufunc calls. That is exactly the whole point of the special work, to get things into exactly that model.

For other situations where things don't neatly fit into the (g)ufunc model, it's a very good question. I personally try to stay in pure Python as much as I can, but acknowledge there are cases where that doesn't work.

from scipy.

Back to mdspan. I tried a small example in https://github.com/ev-br/scipy/commits/splrep_repro_with_mdspan/.
Here f0585b9 adds a small amount of C++ with homemade array wrappers, 874533d drops in the mdspan implementation from gh-20320 and 766cdaa replaces homemade wrappers with mdspan.

Overall the result looks nice, and using mdspan seems easy --- until it isn't.

Minor hiccups:

• indexing of 2D arrays is with operator() not operator[]. For 1D arrays these seem interchangeable, for 2D arrays the compiler complains with something about C++20. My take is round brackets are OK for the foreseeable future.
• Enumerating dimensions is not rank-checked. For a 2D array .extend(8) compiles and gives the same result as .extent(1). So it seems that for N-dim array .extent(M >= N) gives the size of the last dimension?

Larger things are layouts and accessors. In this example I control the arrays, so I know they are always C-contiguous.

• declaring a strided layout with known strides seems possible but I don't (yet) undestand how.
• No idea how to switch on compile-time boundschecking. I guess one needs to supply an AccessorPolicy, but it's very much not clear how.

from scipy.