reshape should be consistent re copying or sharing data

[dmbates]

Currently the documentation for reshape states:

julia> help(reshape)
Loading help data...
Base.reshape(A, dims)

   Create an array with the same data as the given array, but with
   different dimensions. An implementation for a particular type of
   array may choose whether the data is copied or shared.

Copying or sharing should be consistent across array types. Perhaps use reshape! for a version that shares storage and save reshape for a version that copies the array's contents.

[JeffBezanson]

See #2279. I'm against calling it reshape!, since a sharing reshape introduces no mutation, impurity, or side-effect of any kind. I believe I would be ok with making reshape always shared though, and have it simply unavailable for types that can't implement it efficiently. However, the defense of the current behavior is that it is both efficient and generic for purely-functional code. A possible compromise is to keep reshape as-is, but add reshape_shared for code that uses mutation and absolutely depends on data being shared (theory being that mutating code needs to bear all responsibility for and ugliness resulting from mutation).

[denizyuret]

A possible use case for reshape is neural network optimization code. I am trying to reimplement Ng's deep learning tutorial code in Julia (http://ufldl.stanford.edu/wiki/index.php/UFLDL_Tutorial). The optimizer (e.g. NLopt) expects a flat parameter vector. The neural net code is more readable if parameters are interpreted as several matrices. One would like to view the same data in both ways without duplicating, large models are memory bound. When I try to interpret a vector as several matrices I seem to get a "copy":

c = rand(100)
d = reshape(c[1:10], 2, 5)
d[1,1] = 0 # does not effect c[1]

If I try to interpret the whole vector as a matrix I seem to get a "share":

d = reshape(c, 10, 10)
d[1,1] = 0 # does effect c[1]

It seems the culprit is c[1:10]:

d = c[1:10]
d[1] = 1 # does not effect c[1]

Which can be solved using subarray:

d = sub(c, 1:10)
d[1] = 1 # does effect c[1]

However not with reshape:

d = reshape(sub(c, 1:10), 2, 5)
d[1,1] = 0 # does not effect c[1]

Of course I do not know if data duplication happens when I define d or when I try to write on d. If it is the latter, than memory cost would not be an issue with read-only code.

[JeffBezanson]

I don't believe subarrays can be reshaped in-place in general, but it is possible in 1d. That raises the question of whether it's ok for reshape of a subarray to copy the data "sometimes".

[tknopp]

Is there a technical reason that a function reshape! that mutates the array header of its argument cannot be implemented?

[ivarne]

When the shape is part of the type, you can't change the shape without changing the type. If the type of arrays are prone to change, all code that dispatches on array shape will have to do runtime checks to see if someone changed the shape/type of the array. That is exactly the problem we have that makes global variables so slow.

You could of course make reshape! able to change an Array{T,2} from a 1x10 array to a 10x1 array, but that would be rather limiting. Changing size of arrays in place, might also make it harder to have automatic bounds check removal in the future.

[tknopp]

Ah yes, thanks for the explanation! reshape! does not makes sense then.

[JeffBezanson]

I do think we should make reshape always share data. The first change for that is easy: remove the copying fallback. The harder case is, you guessed it, sparse matrices. The sparse method is the only type-specific one in Base that copies data. Should we deprecate that to a combination of copy! and similar? cc @ViralBShah @StefanKarpinski

[tkelman]

An extension of the generic lazy view could maybe be made to work with a different output size? The copying reshape on sparse matrices could be useful in a few places but for consistency it makes sense to give that a different name.

[timholy]

The sparse method is the only type-specific one in Base that copies data.

Yes, but we have

julia> A = rand(10,10);

julia> B = sub(A, 1:8, 2:7);

julia> @which reshape(B, (48,))
reshape(a::AbstractArray{T<:Any,N<:Any}, dims::Tuple{Vararg{Int64}}) at abstractarray.jl:202

which points to the copy version.

Ref #10507.

[JeffBezanson]

Yes, we would need to remove that copy implementation, so reshaping SubArrays becomes an error because it doesn't obey the contract. Then we'd need #10507 to get the function back.

[StefanKarpinski]

I wonder if we shouldn't rethink the whole reshape API more radically. After all, this is really the problem child here. Without it, everything else composes fairly nicely.

[timholy]

@JeffBezanson, I see, I hadn't thought about forcing the user do an explicit copy during an intermediate phase.

For writing generic code, it seems possible people will define a myreshape which doesn't copy for Arrays and does for everything else, but perhaps that's still progress.

[JeffBezanson]

One use case I see a lot is reshape(1:n, i, j). Clearly that doesn't need to be a reshape, but we should really have an array constructor that directly accepts an iterator of values and a shape. Hard to believe we don't have that already --- do we?

[timholy]

Not that I'm aware of.

Over at #15449, reshape(1:n, i, j) is working as a 2D view of the range. If the user never tries writing to it, it's actually nicer than using an Array; the performance is very good (1:n is LinearFast, so you're OK), and you never get cache misses during usage.

However, as soon as you try writing, you're hosed. Unless you first say copy(reshape(1:n, i, j)).

[timholy]

Worth pointing out that this construct is more common in test code than anywhere else.

[JeffBezanson]

That sounds great, and a good example of why this issue is needed. We can fully separate a no-copy reshape of a range from making a dense mutable range.

[JeffBezanson]

Here's an idea: we could give collect the same arguments as similar, and allow it to accept an element type and shape. We have collect(T, itr), but that would be deprecated to collect(itr, T).

[timholy]

👍 Do you throw an error if prod(sz) != length(iter)? Or does it truncate if shorter, and error only if longer?

Since you're thinking about this, another thought: should we consider operators and symmetry? In operations like A+B, allocating the output with similar(A, ...) is obviously asymmetric. We try to handle that now by letting the element type take account of features of both A and B, often making use of promote_op.

However, the "right" type to allocate may depend on both A and B in ways that the element type doesn't capture. For example, for ::Tridiagaonal + ::Tridiagonal, the right type to allocate is Trididagonal, but for ::Tridiagonal * ::Tridiagonal the right type is BandedDiagonal. This is not a distinction that can be conveyed by an operation on the element type.

A potential API (which could be parallel to the existing API) is similar(op, (A, B), sz). The advantage of this API is that one can provide specialized methods. The harder part is knowing how to come up with reasonable fallbacks. Presumably, falling back to Array would be the best choice, esp. if any of the inputs are Arrays.

[andreasnoack]

Way too late here but after having had some fights with ReshapedArrays wrapping SparseMatrixCSC and DArray lately, I really don't think it it is reasonable to have memory sharing to be part of the AbstractArray contract. It's basically to expect that all array types are like Array. ReshapedArray masks the layout of the underlying array and makes it hard to benefit from specialized methods so in generic code you'd basically would have to copy whenever you have done a reshape, dropdims, vec etc.