Refine dynamic borrow checking to track data ranges.

[adamreichold]

Follow-up to #274 as discussed in #274 (comment). Should show that this can work, but still needs needs many unit tests to drive the range data structures into various stages...

[adamreichold]

Needs to be rebased after #302 lands...

[kngwyu]

To be honest, I don't think the current conservative strategy is not that bad. I don't imagine the situation where we want to handle slices this smartly.

[adamreichold]

To be honest, I don't think the current conservative strategy is not that bad. I don't imagine the situation where we want to handle slices this smartly.

I spent some time thinking about this: I do not personally have a use case where this part - supporting non-overlapping views, but not interleaved ones - would make a difference, but I would certainly want to support the use case of interleaving views. For example three-dimensional arrays which are made up of multiple two-dimensional layers that modified concurrently but as separate components is something that I handle regularly, for example as image data or multi-layered raster GIS data. So I see this primarily as a stepping stone towards handling this even more general case. (If we actually can handle that as it might be NP to prove absence of aliasing in the general case given only the data pointer, strides and dimensions.)

My main fear is that once this goes out and people start using PyReadwriteArray to write safe Rust code, is that we get somewhat angry bug reports because it sometimes works and sometimes not because some situations involving views are supported and some are not and it will not be immediately obvious at the call site what the problem is. So we should probably at least make the errors more informative, but just figuring out what the kind of aliasing situation is from data pointer, strides and dimensions is almost the same as solving the problem in a general manner.

At least there is no regression of functionality, i.e. PyRreadonlyArray and unsafe as_array_mut should work as they have until now...

[adamreichold]

So I see this primarily as a stepping stone towards handling this even more general case. (If we actually can handle that as it might be NP to prove absence of aliasing in the general case given only the data pointer, strides and dimensions.)

@kngwyu Thanks to the help of @Darksonn, I think we can significantly improve our approximation and are able to handle the "colors channels of an image array" case.

I am still pondering whether there actually exist situation which would run awful of this level of checking as I have not yet constructed an example that does not alias, but has neither disjoint data ranges nor fulfils the GCD criterion, i.e. that does alias only out-of-bounds.

Note that num-integer was already part of our dependency closure via ndarray.

The MSRV build failed because Iterator::reduce was stabilized in Rust 1.51. Will write that manually using Iterator::fold... |-\

[davidhewitt]

Using GCD is a neat trick; I had an inkling that there would be something like that but didn't have the time to think about it. TBH I'd hope that most cases will now be supported by this implementation 👍

[adamreichold]

TBH I'd hope that most cases will now be supported by this implementation +1

As identified in the unit test, cases where the step size does not divide the dimension along the axis which is sliced are not covered. Handling that would mean recursively solve the linear Diophantine equation implied by the strides of the two arrays.

Personally, I think for now I would like to get this into the wild as part of version 0.17 and then see what sticks as doing the above has significant engineering costs: BorrowKey would need to store all strides, so either Vec<isize> or [isize; 32] and the costly solver would need to run for every possible conflict in the same allocation for every borrow.

But I definitely need to write more unit tests for driving the BORROW_FLAGS through its state space, before marking this ready to be merged... 💦

[adamreichold]

@kngwyu David already had a look, but since we are not in a hurry to release 0.17 (I think we should stick to release in parallel to PyO3.), I would be glad if you could look into this as well now that it hopefully is in good enough shape to do that. Thank you!

I think the GCD-check-based solution is good compromise between effort and accuracy. I am not sure how good the tests are and would be glad if anyone had suggestions on what test cases might be missing.

[kngwyu]

Thank you for digging the problem this deep! I have two questions:

• I think I understand the 1-D case, but what does gcd(stride1, stride2, stride3, ...) mean?
• What do you think of adding a naive range checking using shape?

[Darksonn]

The gcd operator returns the greatest common divisor, that is, the largest integer that is a divisor of all of the numbers. It can be computed as gcd(stride1, stride2, stride3, ...) = gcd(gcd(gcd(stride1, stride2), stride3), ...)

[kngwyu]

Sorry my question was why we are computing gcd of all strides, but I think I understand now. GCD of all strides should be the minimum (pointer) distance to the next element in most cases where the GCD is equal to the smallest stride.

[adamreichold]

What do you think of adding a naive range checking using shape?

I am not sure what you mean exactly? We do have range checking based on the highest and lowest address that is part of the view, so computing shape[i] * strides[i] and adding to the lowest or highest offset depending on whether strides[i] is positive or negative, c.f. the data_range function. This is used as a first condition to reject conflicts without considering the case of interleaving views.

If you mean by checking the slice indices themselves which were used to create the views, e.g. : and 1:7:-3 in array[:,1:7:-3], we do not have access to those, only to the data pointer and the strides which resulted from them.

I think I understand the 1-D case, but what does gcd(stride1, stride2, stride3, ...) mean?

The idea here is that if there are two arrays A and B with data pointers P_A and P_B, dimensionalities D_A and D_B, dimensions N_{1,A}, ..., N_{D_A,A} and N_{1,B}, ..., N_{D_B,B} and strides S_{1,A}, ..., S_{D_A,A} and S_{1,B}, ..., S_{D_B, B}, they overlap if and only if the following linear Diophantine equation

P_A + S_{1,A} * i_{1,A} + ... + S_{D_A,A} * i_{D_A,A} = P_B + S_{1,B} * i_{1,B} + ... + S_{D_B,B} * i_{D_B,B}

has at least one solution that is within bounds, i.e. that fulfils 0 <= i_{k,x} < N_x for each 1 <= k <= D_a and x in {A,B}.

This equation has at least one solution if and only if

GCD(S_{1,A}, ..., S_{D_A,A}, S_{1,B}, ..., S_{D_B,B}) divides P_A - P_B

This is the condition that the code is checking and if it does not hold considers them non-conflicting. If it is fulfilled, then it considers them conflicting which is an approximation as all solutions could still be out of bounds (which is the case in the third example).

One possible extension to the above would be to actually compute these solutions (meaning to actually compute their bounds) by recursively solving the equation via intermediate equations with two-variables. But I think this would add considerable code complexity, runtime (on top of the around 500 ns this currently takes) and memory (we would need to store all strides and dimensions, not just the data range and the GCD of the strides) which is why I would like give this a try and see how many complaints it raises. (Especially since we can release a refinement that allows more cases to pass as a point release as it is not a breaking change.)

[adamreichold]

the minimum (pointer) distance to the next element in most cases where the GCD is equal to the smallest stride.

Indeed, and when it is not equal to the smallest stride, this fact can be used to handle interleaving views, e.g.

import numpy as np
a = np.zeros((10, 10))
a.strides # (80, 8)

v3 = a[:, ::2]
v3.strides # (80, 16)

v4 = a[:, 1::2]
v4.strides # (80, 16)

has a pointer difference of 8 but a GCD of all strides of 16 which implies that those views do not alias.

[kngwyu]

We do have range checking based on the highest and lowest address that is part of the view

Sorry, I just overlooked that part.

One possible extension to the above would be to actually compute these solutions (meaning to actually compute their bounds) by recursively solving the equation via intermediate equations with two-variables

Yeah, I think the GCD of all strides covers most cases

[kngwyu]

Thank you, pretty good job, but it seems like there is still a corner case

[kngwyu]

Thank you, I just added a few comments on the document.

[kngwyu]

Thank you for this difficult job, merge this at your own time.