Composition bug when chaining 4 LinearMaps

[cako]

I ran into a bug when composing 4 maps using LinearMaps 1.0.4. When running the following script in Julia 0.6.3, I get ERROR: cannot resize array with shared data.

import LinearMaps
dz, dx = 15, 15
X = 0:dx:1000; nx = length(X)
Z = 0:dz:1000; nz = length(Z)

ricker(to, f) = (1 - 2pi^2 * f^2 * to.^2) .* exp.(-pi^2 * f^2 * to.^2)
rick = ricker(Z-Z[div(nz-1,2)+1], 0.01); 

function conv_fwd(mod, wavelet, nz, nx)
    return hcat( (conv(reshape(mod, nz, nx)[:,ix], wavelet)[div(nz-1,2)+1:div(3nz-1,2)] for ix in 1:nx)... );
end
function conv_adj(mod, wavelet, nz, nx)
    return hcat( (conv(reshape(mod, nz, nx)[:,ix], wavelet[end:-1:1])[div(nz-1,2)+1:div(3nz-1,2)] for ix in 1:nx)... );
end

C = LinearMaps.LinearMap(x -> conv_fwd(x, rick, nz, nx)[:],
                         x -> conv_adj(x, rick, nz, nx)[:],
                         nz*nx, nz*nx)
x = ones(nz*nx)
(C*C*C*C*x)[10] == 0.888757975075548

The culprit is a resize! in composition.jl (line 93). I propose to simply create a new array at each composition multiplication (in a commit soon to follow). It is true this is not the most efficient solution, however, people using LinearMaps will tend to have very expensive operators such that the cost of allocation is negligible compared to the cost of multiplication.

[dkarrasch]

Hi @cako , that was a tough one. Generally, this functionality is tested, so it came as a big surprise to me. The issue seems to be here with your function definition. As a MWE, take

src = ones(nz*nx)
dest = conv_fwd(src, rick, nz, nx)[:]
resize!(src, 20)

which yields

ERROR: cannot resize array with shared data

on my machine. So, applying your function to a src variable does something to it such that it cannot be resized afterwards. Regardless of whether that is meaningful, it is strange; after all, you would like to be able to do whatever you want to src. There are three options here:

1. This is a Julia 0.6 bug, then you may wish to find an even simpler MWE and report it there. My gut feeling is, however, that this was intentionally.
2. You can reshape(copy(mod), nz, nx). That solves the problem on my machine and Julia 0.6.3.
3. You work with Julia 0.7 and our brandnew LinearMaps update, where it works as you wrote (after using DSP for conv).

For that reason, I would not like to give up the semi-inplace resize!, which effectively facilitates a as-low-as-possible-memory matvec multiplication. Hope this helps.

[cako]

Hi @dkarrasch, awesome, thanks for looking into this.

Using your solution 2 for me is absolutely fine, so I'll stick to that!

[Jutho]

Having a collaborator is already paying off on day one ! 🥇

[dkarrasch]

Thanks. On my end, I wasn't sure if that was a wise decision, on day one. :-))))))))) But my Ph.D. student helped me, so just in case, I'll pass the issues over to him and keep all merits and medals. ;-)