Add fallback when resize! does not work

[cako]

I have been thinking about the issue #24 (using resize! within the composition operator) and have come to the conclusion that users should not be expected to always provide "clean" arrays which can be resized. Take my case for example, I had to essentially change my function by adding a couple of copy commands to ensure resize! was working. I was only able to do this after guidance. In fact, I have a more complex example (which I won't share because it is pretty long and specific) that I simply cannot get to work with the current code. Copying data both at the input and at the output do not work. I have another (somewhat contrived) example demonstrates how this can happen. It can also be used for code coverage if the commit is accepted.

So either we can accept that some people's code will mess with the array in such a way that it cannot be resized, or we make them accept the fact that the library is not for them. I don't think the latter is beneficial for anyone!

On the other hand, I fully appreciate that resize! is much faster than array creation. I benchmarked resizing an array of size 1000 to several different sizes, ranging from 100 to 100_000. The cost of only the resizing is in orange, while the cost of allocating from scratch is shown in green.

Clearly allocating from scratch grows linearly, while resizing may eventually be linear as well, but it is offset by how much has already been allocated.

So I propose a very simple solution that I think will make everyone happy: we resize when we can, and when we can't, we allocate. I have added a simple try/catch block to do this. As shown in the figure above, this "fix" (blue line) retains decent performance for larger arrays, and is still much better than allocating for smaller arrays. It also has the obvious benefit of not breaking anyone's code!

The commits are very simple, simply wrap all resize! calls with a try/catch block which in the catch statement allocates an array.

[codecov]

Codecov Report

Merging #31 into master will decrease coverage by 0.96%.
The diff coverage is 80%.

@@            Coverage Diff             @@
##           master      #31      +/-   ##
==========================================
- Coverage   98.59%   97.63%   -0.97%     
==========================================
  Files           7        7              
  Lines         142      211      +69     
==========================================
+ Hits          140      206      +66     
- Misses          2        5       +3

Impacted Files	Coverage Δ
src/identitymap.jl	100% <ø> (ø)	⬆️
src/composition.jl	95.83% <80%> (-4.17%)	⬇️
src/wrappedmap.jl	87.5% <0%> (-12.5%)	⬇️
src/linearcombination.jl	100% <0%> (ø)	⬆️
src/LinearMaps.jl	97.91% <0%> (+0.09%)	⬆️
src/functionmap.jl	100% <0%> (+5.55%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 164f67e...09ccffe. Read the comment docs.

[dkarrasch]

@cako , thanks for working on this. This is not a final opinion, but one thing I don't like about the catch-try block is that it may silently switch to the catch-case, because there is something wrong with the code that the user would want to (or should) know. Hardcoding error catching like this has a somewhat unpleasant flavor.

I acknowledge, though, the following (mis-)behavior:

x = rand(10)
y = reshape(x, 2, 5)
resize!(x, 8)

fails for the same reason. So x is aware of the fact that its data is shared with y. In our context, this is bad since, in a way, a linear map that calls reshape of its argument mutates its argument, which the user (and LinearMaps.jl) don't expect. Let's hear what @Jutho has to say.

[Jutho]

reshape is indeed a very common use case to trigger this, and I agree this is a bug.

I think something like this can detect whether an Array is shared, and thus, if resize! will fail

_flags(A::Array) = unsafe_load(convert(Ptr{UInt16}, pointer_from_objref(A)), 9)
_isshared(A::Array) = !(_flags(A) & 0x4000 == 0x0000)

I am however a bit worried if this depends on the architecture (i.e. 32 vs 64 bits), in particular the value 9 in the first function will probably have to be different on 32 bit. So probably it needs to be

if Sys.WORD_SIZE == 64
 _flags(A::Array) = unsafe_load(convert(Ptr{UInt16}, pointer_from_objref(A)), 9)
elseif Sys.WORD_SIZE == 32
 _flags(A::Array) = unsafe_load(convert(Ptr{UInt16}, pointer_from_objref(A)), 5)
else
 error("Unknown word size")
end
_isshared(A::Array) = !(_flags(A) & 0x4000 == 0x0000)

But I am not sure and don't have a 32 bit machine to test this.

[dkarrasch]

That sounds like a plan: detecting the specific reshape problem and adjusting the code for this case, instead of a generic try-catch-block. Let's see what the experts say in the push!/pop! issue...

[Jutho]

Another way out is that we pass view(dest, :) to the user's linear map. This can never lead to shared data, but upon reshaping it creates a ReshapedArray. Depending on what the linear map does, this might also have a small performance cost. Furthermore, this assumes that the user did write a function that accepts AbstractVector arguments instead of just Vector.

[cako]

Thanks for the comments, I think we are making progress!

@dkarrasch I think in general it is definitely better if we can use something other than try/catch, but in this case it might be warranted. See the end of this comment for my reasoning!

@Jutho let me see if I got this right: every Julia Array has a flag which says if it will be resizable (shared) or not. Since Julia Arrays are implemented using C structs, we have to somehow get the flag from it. So, first you get the Ptr which points to that array in memory, and then convert it to Ptr{UInt16} with

convert(Ptr{UInt16}, pointer_from_objref(A))

Now, this pointer is the memory location of jl_array_t as defined in julia.h.

This is where I start getting confused, in 32 bits, void* and size_t should have size 4 bytes each. In 64 bits, they should have 8 bytes each, placing the memory location of jl_array_flags_t for the respective architectures at 9 and 17. However, you have indices 5 and 9 for that memory location. Why is this? After that, you load jl_array_flags_t using unsafe_load and check the if shared bit indicator flag at 0x4000 is on.

With regard to speed, I ran your solution against the previous benchmarks, and this is what I got:

The conclusion is that the both try/catch and flag-checking have virtually the same performance of a bare resize!. The try/catch block, however, slightly penalizes performance on very small arrays (<100). So, in terms of performance, I think either option is fine.

In terms of "cleanness" although it is a matter of opinion, I quite like the try/catch block. Although in general all-inclusive try/catch blocks are not great as @dkarrasch has pointed out, in this case I think it make sense to use it. The alternative solution provided by @Jutho hit a snag as pointed out by @chethega in issue 24909 (namely, we don't know if it is possible to estimate STORE_ARRAY_LEN at runtime). In addition, the try block has only one function (resize!) operating on an internal array to LinearMaps.jl. The user does not know it is there, and does not care. The only thing the user cares about is that the resulting array is returned. If that array cannot be resized, the user will immediately realize it is because of the action of one of their functions. If they care about it, they can try and fix it!

Let me know what you think. On a separate note, if you'd like, I can start writing my example as a test and add the commit to this pull request.

[dkarrasch]

Regardless of how we're going to fix this issue, I think adding the test would be good. This is a generic problem that we would want to be able to solve.

Unfortunately, there is not much to be learned from the implementation of resize!, the error is thrown from the inside of some ccall. It is not thrown when you have arrays with shared data that happen to keep the size. For these cases we need to make sure not to add significant overload. Is that what your simulations are showing, @cako? Are you taking some good old "standard" case and measure runtime with and without try-catch etc.? Could you perhaps share your benchmarking code in the gist?

[Jutho]

This is where I start getting confused, in 32 bits, void* and size_t should have size 4 bytes each. In 64 bits, they should have 8 bytes each, placing the memory location of jl_array_flags_t for the respective architectures at 9 and 17. However, you have indices 5 and 9 for that memory location.

@cako, the numbers 5 and 9 are measured in units of sizeof(UInt16), because the pointer is a Ptr{UInt16}. I convert the pointer to a Ptr{UInt16} because I know from the source code of array.c that the flags field in the array struct is 16 bits wide, with the isshared flag being on bit 15 out of 16, see here: https://github.com/JuliaLang/julia/blob/master/src/julia.h#L150

But I agree that this approach is too fragile and clearly advised against.

The try .... catch solution is certainly fine with me, I think the only other solution is giving the user views into the array. Or a unsafe_wrap(Array, pointer(dest), ...), but that's again becoming "unsafe"!

[cako]

@dkarrasch cool, I will see if I can clean up/simplify the test and commit it. As for the benchmarks, you can find them here. Basically I write functions which allocate a size 1_000 array, then resize, with optional try/catch or flag testing. I benchmark that function, then benchmark only allocating the size 1_000 array, and subtract the latter from the former. That gives me the benchmarks I used for the graphs above. I've run it a few times and I reliably get slighly worse performance for small arrays (<1000), and essentially the same performance for larger arrays.

@Jutho and @dkarrasch Perhaps a compromise is catch only the shared data error:

try
    resize!(dest, size(A.maps[n], 1))
catch err
    if err == ErrorException("cannot resize array with shared data") # use `==` instead of `isa` because `resize!` throws a general exception
        dest = Array{T}(undef, size(A.maps[n], 1))
    else
        rethrow(err)
    end
end

[dkarrasch]

Yes, I think that's a good compromise. In my own (unrelated to this package) code, I used try-catch-blocks (perhaps in a bad way), and had difficulties of finding mistakes after changes in the code, and just got Nothing. So that's my own personal bad experience with it. That's why I think it is good to be as specific as possible about the error that one wants to catch (as you suggested), in order to not cover up other issues which the user should consider and solve.

[dkarrasch]

Thanks @cako for the nice tests. I made some minor changes. The transpose/adjoint of a CompositeMap is a CompositeMap again, and therefore gets handled by the A_mul_B!, and not the A*_mul_B!s. I hope this makes coverage look better.

[dkarrasch]

@Jutho: While @cako was working on his branch, our master branch diverged slightly (in runtests only). How do I make sure that our latest tests and @cako's new tests both survive the merge? Is that what rebasing is about? With the union of the tests, we should be able to keep excellent coverage.

[Jutho]

Since Github still states: "This branch has no conflicts with the base branch" I think this can still safely be merged in master and all edits (commits) will survive. But we can also try to rebase, this shows up even on Github, if you fold out the "Squash and merge" drop down menu.

[dkarrasch]

I wasn't sure what is meant by "base branch", if it's our current master or the one where @cako branched off. I'll simply squash, merge, and tag a new fix release 2.2.2 then.

[Jutho]

Alternatively, @cako could rebase locally on the latest master and then push (maybe push -f for forced) here

[cako]

@dkarrasch Thanks for fixing the transpose/adjoints, I couldn't figure out how to do it!

I can do the rebase locally, but only later today!

[dkarrasch]

No worries. It's a bit tricky (and in fact, I was believing up until recently that those methods are not needed at all ;-) ).

The rebase is mainly to have all tests together, let code coverage run again and get good grades.