array_t: overload operator() for multi-dim indexing

[fnrizzi]

Description

Currently, array_t supports direct subscripting via:

template<typename... Ix> const T& at(Ix... index) const;
template<typename... Ix>       T& mutable_at(Ix... index); 

both of these always check preconditions on the rank (ndim) and bounds, and are methods named with a specific name.

This PR proposes to overload operator():

template<typename... Ix>  const T& operator()(Ix... index) const;
template<typename... Ix>        T& operator()(Ix... index);

and checking preconditions only in debug mode.

The main advantage of this PR is that it would make it easier to write generic functions like this:

template<class T>
void foo(T & matrix){
  matrix(0,1) = 454.;
}

void caller_cpp(){
  Eigen::MatrixXd A(10,11);
  foo(A);
}

using eigen_mat_layleft_t = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
using eigen_ref_mat_layleft_t = Eigen::Ref<eigen_mat_layleft_t>;

void caller_py(eigen_ref_mat_layleft_t & mat){
  foo(mat);
}

This is especially useful considering the interoperability of pybind11 and Eigen.

Suggested changelog entry:

TBD

[rwgk]

(FYI: I just clicked the Approve & Run button for GitHub Actions without actually looking at the code. Will look asap.)

[fnrizzi]

@rwgk i added a test too, which I am building locally now - will push once i have it

[rwgk]

@rwgk i added a test too, which I am building locally now - will push once i have it

Sounds good, definitely add tests first.

After the new tests are place:

The code duplication is now 4x, which is not good. It's essentially only one line, but it's non-trivial. I'd try hard to have only 1x, even if the total number of lines of code stays the same (or even if it increases); the quality of implementation is more important.

*(static_cast<const T *>(array::data()) + byte_offset(ssize_t(index)...) / itemsize());

Question: Is it guaranteed from the context that byte_offset(ssize_t(index)...) is an integer multiple of itemsize()?

Why I'm asking: If that's not guaranteed, I'd look into adding something like

assert(byte_offset(ssize_t(index)...) % itemsize() == 0);

Note for completeness: if NDBUG is defined, that's a no-op.

Not sure: maybe just that calculation in a new helper function will do. Duplicating array::data() + elem_offset(index...) would seem fine to me.

I'm also wondering about naming. It's unfortunate that the safe alternative (at()) is relatively easy to pin-point, but the potentially dangerous one is very difficult to pin-point without sophisticated tooling (I'm used to looking at millions of dependencies). Therefore I'd lean toward avoiding operator() for this purpose. I'd look around a little for a name already established elsewhere. If there is nothing obvious, I'd use something like elem_reference() (long == good for something potentially dangerous).

Finally, the bounds-checking functions could/should call the non-bounds-checking counterparts. Similar idea as before: clean abstractions & the quality of implementation matter most (in the long run).

[fnrizzi]

Agree on fixing things like having mutable_at call this etc. But wanted to first discuss the choice of overloading operator() below

Therefore I'd lean toward avoiding operator() for this purpose. I'd look around a little for a name already established elsewhere. If there is nothing obvious, I'd use something like elem_reference() (long == good for something potentially dangerous).

Changing the name defies the purpose of this, because operator() is something Eigen already has (and most multi-dim containers do too for indexing like Blaze https://bitbucket.org/blaze-lib/blaze/wiki/Getting_Started or Armadillo https://arma.sourceforge.net/docs.html#element_access) so that is why i was proposing this.
If a specific name like you suggest is chosen, the generality of my foo function goes away :)

[rwgk]

Changing the name defies the purpose of this, because operator() is something Eigen already has (and most multi-dim containers do too for indexing) so that is why i was proposing this.

Oh, got it.

In that case, what do you think about

#ifndef NDEBUG
bounds check here
#endif
compute the reference here

I think that could stave off a lot of debugging misery, considering that the out-of-bounds access may be triggered by complicated generic code.

[fnrizzi]

so you are suggesting doing this?

template <typename... Ix>
T &operator()(Ix... index) 
{
#ifndef NDEBUG
bounds check here
#endif
   return *(static_cast<T *>(array::mutable_data()) 
                 + byte_offset(ssize_t(index)...) / itemsize());
}

if so, I agree on this

[rwgk]

so you are suggesting doing this?

Yes.

[fnrizzi]

agreed, let me modify the PR

[rwgk]

I think I misunderstood the existing code (silly in hindsight).

This is just checking the number of dimensions:

• pybind11/include/pybind11/numpy.h

  Lines 1096 to 1098 in 2fb3d7c

  	if ((ssize_t) sizeof...(index) != ndim()) {
  	fail_dim_check(sizeof...(index), "index dimension mismatch");
  	}

But that's not what I would think of as a bounds check. That one happens here:

• pybind11/include/pybind11/numpy.h

  Line 981 in 2fb3d7c

  check_dimensions(index...);

I.e. if we also want to eliminate the bounds check, we cannot use bytes_offset(), but we have to use detail::byte_offset_unsafe().

This makes it a slightly more tricky refactoring, but the work you did already is going a good way towards that.

I'd use the NDEBUG idea for the bounds check as well.

[fnrizzi]

@rwgk my bad for writing a description that was imprecise. My origianl main intent was to add operator() and when i said bounds was imprecise. I fixed the description now.

I added checks for the rank (ndim) and bounds and I believe i also added proper tests.

[rwgk]

I could only glance through quickly right now. What I saw looks great. I'll look carefully in the morning.

[fnrizzi]

CI failure is unrelated (failed to fetch boost)

[EthanSteinberg]

@fnrizzi This looks great, thanks for the addition!

One important note is that we technically already had this feature, but it was a poorly (IMHO) designed thing involving proxies. See https://github.com/pybind/pybind11/blob/master/docs/advanced/pycpp/numpy.rst#direct-access

What I think we should do is "deprecate" those proxies, remove them from the documentation and explicitly mark them as deprecated and encourage people to use this newer, much simpler and easier to understand API. What do you think @rwgk?

At the very least, let's remember to add this new feature to the documentation https://github.com/pybind/pybind11/blob/master/docs/advanced/pycpp/numpy.rst once you and @rwgk think it's ready.

[rwgk]

One important note is that we technically already had this feature, but it was a poorly (IMHO) designed thing involving proxies.

Oh ... I wasn't aware of that TBH. (I was just looking at this PR from a general best-practice point of view. I never actually used numpy.h myself in a serious way.)

What I like about the proxy approach is that it makes the potentially dangerous access explicit: I can go and look for pybind11 numpy unchecked in a large codebase and probably pin-point most use cases quickly. — I just tried that search and in our gigantic codebase I only found 58 hits, i.e. it's definitely working as a way to quickly narrow down where to look. BTW, this kind of thing is what I spend a lot of my time on.

Now, it's evidently (this PR) not obvious that the feature exists for someone coming in new. Could we help with that, but keep the established proxy? For example, could this work?

    template <typename... Ix>
    const T &operator()(Ix... index) const {
        static_assert(!std::is_same<T, T>::value, "Please use pybind11::array_t::unchecked<>");
    }

Similarly for mutable_unchecked<>.

I think that would be far better than introducing a second way to get direct access. Because, if there are two ways, both will get used, the "old" way will take a decade+ (no kidding) to get flushed out 90%, and a virtually infinite amount of time to disappear completely. All the while we have to maintain two solutions instead of one, people coming in new will go huh, then find themselves ~forced to use sometimes one sometimes the other depending on the environment, etc.

[fnrizzi]

I actually used those before but I remember having some issues. Obviously now I don't remember off the top of my head the issues...
but that was one of the reasons I resorted to adding operator().

Btw, the proxies IMMO adds an additional layer that one has to consider.
If the main concern is about safety, then we still check in debug mode.

[fnrizzi]

Btw I actually would have proposed a multi index subscript operator [] that would be intuitive and also equivalent to python syntax.

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p2128r6.pdf

I think that would probably be my preferred solution but that is also not super common.

[rwgk]

To explain where my thinking is coming from:

• I find that my life is consumed by mopping up little messes created by divergence of some sort of another. This is because I'm looking at a gigantic code base (Google's mono repo), and I'm responsible for keeping everything working as I evolve core features. (Note for completeness: I can make changes anywhere, although they all need owner approval from other engineers.)

• Looking at the world as a whole, my situation is probably very unusual. For me as one person. But I believe not actually for the world in aggregate, I just experience a common problem more concentrated. The churn caused by random little divergences left and right is very real also for the world at large, only usually the pain is more spread out, less personal.

With that background, coming back to this PR, I think we need to balance:

• What do we gain by introducing a 2nd way (divergence) to get direct access?

vs

• The churn we create as a side effect?

If we had a green field, I'd decide against the proxies, for sure. But we do not have a green field. The proxies are established, they have something akin to the first mover advantage. IMO opinion we need a very strong, very convincing argument to deprecate the proxies and accept the churn.

I actually used those before but I remember having some issues. Obviously now I don't remember off the top of my head the issues...

Could you please try to reproduce and explain here?

[rwgk]

Re https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2021/p2128r6.pdf

Thanks for pointing out, I wasn't aware of this before.

I want to share a long-term experience, FWIW.

A very long time ago I implemented an "array family" library, with all bells and whistles, what was possible with C++98. What's called "mdspan" in the referenced paper I covered under array "accessors":

https://github.com/cctbx/cctbx_project/tree/master/scitbx/array_family/accessors

Later it turned out, the sophisticated accessors (flex_grid.h is the most abstract, c_grid.h more special/common) where hardly ever used in C++ (c_grid.h more than others). I think that's because ultimately the only good reason to put up with C++ is to maximize performance. Consequently, almost every single time, I'd find myself unrolling the loops to cut out the abstraction penalties.

[EthanSteinberg]

@rwgk I think you are right here. Seems best to avoid churn by keeping the proxies, but adding those helper asserts to guide people towards the proxies if they need that feature.

[fnrizzi]

I understand all these viewpoints , but I like to point out that the proxies are (1) not intuitive (2) a deviation from any "linear algebra/container like" thing, and (3) add an additional layer making generic code (see my example at top) even more verbose. So it is not just a matter of divergence, but usability too.

Speaking as a user of pybind11, I personally prefer having my fork with operator () than using the proxies.

[rwgk]

I'm a believer in numbers and reason. Personal preferences are are a rare luxury for greenfield situations.

(see my example at top)

Sorry I'm failing to make the connection, maybe because I never wrote original Eigen- or pybind11/numpy.h-based code. The code currently in the PR description doesn't seem to involve pybind11/numpy.h at all. Could you please provide an example (A) with pybind11 as is, and (B) the equivalent with operator() for direct access?

This is to estimate the percentage of code size reduction enabled by operator() for direct access, in a typical situation (we have to assume there is usually some amount of code around the lines involving operator()).

[henryiii]

What about operator[], only if C++23+ is used? In C++23, this supports multiple args, and seems like the ideal operator in the future for this sort of thing. It also is clearly different than calling an object and triggering __call__. Adding () now then [] later seems like more churn. If we could integrate with the new mdspan nicely, that would be a benefit, and inline with pybind11's attempts to integrate well with the C++ standard library.

The proxies allow you to get checked or unchecked access, and allow you to select mutable or immutable access, which is nice. Given the official way to interact with ND arrays is now std::mdspan in C++23, which itself is a proxy, I don't think this is unusual at all in C++.

(PS: personally, I'd love some way to iterate over every element of arbitrary ND arrays - basically to allow making custom ufuncts in C++ easy. I recently found myself needing that to implement float -> int conversion that didn't round -.x to 0.)