smooth!() causes unexpected array size change

[kura-okubo]

Hi Tim,

smooth!() in src/tools.jl might cause unexpected array size change for 1d Array as follows:

using SeisNoise
t = collect(0.0:0.05:60.0)
A = sin.(0.2*pi.*t)

# first time
smooth!(A);
len1 = length(A)
# second time
smooth!(A);
len2 = length(A)
# third time
smooth!(A);
len3 = length(A)

OUTPUT:

julia> println((len1,len2,len3))
(1208, 1215, 1222)

I found that the movingaverage() could cause this issue; the modification of A::AbstractArray seems to affect the returned value (although it shouldn't be without !). So this might be Julia-related issue.

For the moment, I modified them as follows:

function smooth_debug!(A::AbstractArray, half_win::Int=3)
    if ndims(A) == 1
        return movingaverage_debug(A,half_win)
    end

    Nrows, Ncols = size(A)

    for ii = 1:Ncols
        A[:,ii] .= movingaverage_debug(A[:,ii],half_win)
    end
    return nothing
end

function movingaverage_debug(U::AbstractArray, half_win::Int=3)
    #NOTE: avoid to rewrite A
    A = deepcopy(U)
    #----------#
    prepend!(A,A[1:half_win])
    append!(A,A[end-half_win:end])
    N = length(A)
    B = zeros(eltype(A),N)
    window_len = 2 * half_win + 1
    s = sum(A[1:window_len])
    B[half_win+1] = s
    for ii = half_win+2:N-half_win
      s = s - A[ii-half_win] + A[ii+half_win]
      B[ii] = s
    end
    B ./= window_len
    return B[half_win+1:end-half_win-1]
end

And it's solved:

B = sin.(0.2*pi.*t)
# first time
smooth_debug!(B);
len1 = length(B)
# second time
smooth_debug!(B);
len2 = length(B)
# third time
smooth_debug!(B);
len3 = length(B)

OUTPUT:

julia> println((len1,len2,len3))
(1201, 1201, 1201)

Can you manage that until we figure out the call by reference in Julia?

Kurama

[tclements]

Thanks for pointing this out. I think the smooth function as currently implemented assumes a 2D array. I added a new version of the smooth function to the GPU branch a54a642. This removes the dependence on movingaverage and should work for any dimension. It's also faster (for 2D arrays) and works on the GPU! Since this is a bug, I'll push this commit to master.

Here is the new version of smooth:

function smooth!(A::AbstractArray, half_win::Int=3, dims::Int=1)
           T = eltype(A)
           window_len = 2 * half_win + 1
           csumsize = tuple(collect(size(A)) .+ [i==1 for i in 1:ndims(A)]...)
           csum = similar(A,T,csumsize)
           csum[1,:] .= zero(T)
           csum[2:end,:] = cumsum(A,dims=dims)
           A[half_win+1:end-half_win,:] .= (csum[window_len+1:end,:] .- csum[1:end-window_len,:]) ./ window_len
           return nothing
end

smooth(A::AbstractArray,half_win::Int=3, dims::Int=1) =
             (U = deepcopy(A);smooth!(U,half_win,dims);return U)

Now it should keep the same size

smooth!(A);
len1 = length(A)
# second time
smooth!(A);
len2 = length(A)
# third time
smooth!(A);
len3 = length(A)

println((len1,len2,len3))
(1201, 1201, 1201)

[tclements]

I added the new version of smooth to the master branch. Let me know if you have any problems.

[tclements]

This is fixed in the newest branch. Closing.