Assigning an array to a view

General Usage
1

I am trying to use views to write generic code that can work either with a single vector as an input, operating over numbers, or taking a matrix and operating over vectors (columns of the matrix). Since this is a pretty common need, I am also interested in hearing about how to write code for this kind of thing. But the question is really about a problem I got trying to use views.

The idea is to write a function taking the input and output that I treat as one-dimensional collections of elements. Then for vectors, it all just works. For matrices, I first create the output matrix and then vectors with views for each column of the matrices.

The problem I have is that I can’t assign from a vector to the view, even though I can assign from a view. What am I missing here?

a = rand(2, 5)
b = Array{eltype(a)}(undef, size(a)...)
av = [view(a, :, k) for k in 1:size(a, 2)]
bv = [view(b, :, k) for k in 1:size(b, 2)]

av[2] - av[1]
bv[1] = av[2]
bv[1] = av[2] - av[1]    # doesn't work
bv[1] = rand(2)    # doesn't work

Here’s a shorter example on the terminal

julia> b = zeros(3, 5)
3×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> bv = [view(b, :, k) for k in 1:size(b, 2)]
5-element Array{SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}:
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.0]
 [0.0, 0.0, 0.0]

julia> bv[1] = rand(3)
ERROR: MethodError: Cannot `convert` an object of type Array{Float64,1} to an object of type SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}
Closest candidates are:
  convert(::Type{T<:AbstractArray}, ::T<:AbstractArray) where T<:AbstractArray at abstractarray.jl:14
  convert(::Type{T<:AbstractArray}, ::LinearAlgebra.Factorization) where T<:AbstractArray at /home/user/src/julia/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/factorization.jl:46
  convert(::Type{T}, ::T) where T at essentials.jl:154
  ...
Stacktrace:
 [1] setindex!(::Array{SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}, ::Array{Float64,1}, ::Int64) at ./array.jl:767
 [2] top-level scope at REPL[78]:1

julia>
2

The problem is that eltype(bv) is SubArray{… i.e. a view, while rand(5) is an array. If you had bv = Vector[view(b, :, k) for k in 1:size(b, 2)] then you could do this, as eltype(bv) == Vector{<:Any}.

You can also copy into the views, bv[1] .= rand(3) or copyto!(bv[1], rand(3)). This changes b.

2 Likes
3

OK that fixes it, but the broadcast means the code will not work anymore with a vector of base types like Float64. Is there a way to make an assignment for both cases?

4

I guess you can either have an array of Vectors (if this is what you need to insert), or a array whose eltype is suffiently broad to take both views and arrays:

bc = [b[:, k] for k in 1:size(b, 2)]

bv = Any[view(b, :, k) for k in 1:size(b, 2)]

(What I said above with Vector[view(... actually converts to vectors, I believe.)

In both cases replacing these, like bv[2] = rand(3), won’t alter b.

If the vectors are small, then your problem may also be a candidate for StaticArrays.

5

The idea is that I want to read from and write to the original matrix columns. I can’t copy the columns into a new array of arrays. Assigning is supposed to change the original matrix. By making it an Array of Any, when we assign it replaces the view in bv with the new array instead of updating the original matrix b.

6

First of all, I should warn you that this code:

bv[1] = av[2]
bv[1] = av[2] - av[1]    # doesn't work
bv[1] = rand(2)    # doesn't work

wouldn’t do what you wanted, even if it did work. Assigning to bv[1] won’t actually modify your output matrix b at all. It just replaces one element of bv with some other value.

Can you perhaps explain a bit more about why your code needs to operate on both vectors and matrices? Using a matrix as a way to store a collection of vectors is a common pattern in languages like Matlab and python/numpy, but it’s much less common in Julia because it’s often just as efficient to use a vector-of-vectors in Julia.

In particular, if your data consists of a large number of relatively small columns (where each column has less than, say, 100 elements), then a Vector of SVectors from StaticArrays will generally be extremely efficient and easy to work with.

For example:

julia> using StaticArrays

julia> function f!(b, a)
         b[1] = 2 * a[1]
       end
f! (generic function with 1 method)

julia> a = [SVector(1, 2), SVector(3, 4)]
2-element Array{SArray{Tuple{2},Int64,1,2},1}:
 [1, 2]
 [3, 4]

julia> b = similar(a);

julia> f!(b, a)
2-element SArray{Tuple{2},Int64,1,2}:
 2
 4

julia> b
2-element Array{SArray{Tuple{2},Int64,1,2},1}:
 [2, 4]
 [0, 0]
7

On the other hand, if you absolutely must write code that is generic over vectors and matrices, treating a matrix as a collection of columns, then I would suggest taking advantage of multiple dispatch to separate out the behavior that must be different between the two implementations. For example, you might write your top-level function something like this:

function foo!(b, a)
  for i in my_indices(a)
    a1 = my_element(a, i)
    # compute b1 however you see fit
    my_store!(b, a1, i)
  end
end

Then you can define my_indices, my_element, and my_store!to operate differently on vectors and matrices:

my_indices(x::AbstractVector) = eachindex(x)
my_indices(x::AbstractMatrix) = 1:size(x, 2)

my_element(x::AbstractVector, i) = x[i]
my_element(x::AbstractMatrix, i) = @view(x[:, i])

my_store!(y::AbstractVector, x, i) = y[i] = x
my_store!(y::AbstractMatrix, x, i) = y[:, i] .= x

This is a really powerful pattern in Julia: write the high-level logic with as much common code as you can, and then use helper functions to implement the details that vary for different types of inputs.

3 Likes
8

Indeed! It was already doing the same thing as I said in the previous message… But then with .= it does update the original matrix.

9

Thanks, that looks more like it!

The idea is that the data has to sit in these matrices at some point, and if I put it in intermediate arrays I fear it might create some overhead.

It’s a bit of a pity I need to come up with a special method just because we need either = or .=, would it be bad if I just defined = to do .= for views?

10

Have you tested it? A vector of vectors will have some overhead, yes, but it may be irrelevant in the context of your actual code. And a vector of SVectors may have no overhead at all. Benchmarking and testing a few implementations will be much more informative than speculating.

Yes, that would definitely be bad, as (if it were possible) it would break all sorts of unrelated code that relied on the old definition of =. It’s also not possible, since = is not a generic function you can overload.

On the other hand, defining your own function (like the my_store! above) would be just as effective and would not introduce any risk of breaking unrelated code.

3 Likes
12

Maybe the whole idea of using a function as a vectorized version is not really good. I’ll just move to loop externally, and if that was really needed than something like you suggested must be the way to go. This is all perhaps a good illustration of how Julia really is different from Numpy and Matlab, and for the good…

13

I’m giving it another try, and I must be doing something wrong. I want to apply this vector function over the columns of a matrix. And I also have an in-place version that I don’t know if it is worth using, but I want to test it. What would be a good way to do this? What is “the Julia way”?

## Original function, just a mildly complicated operation over the input vector.
function myfun(a::Array{T,1}) where T
    out = Array{T}(undef, size(a)...)
    out[1] = 3 * (a[2] - a[1])
    out[2:end-1] = 3 * (a[3:end] - a[1:end-2])
    out[end] = 3 * (a[end] - a[end-1])
    out
end

## An abomination, I don't want to rewrite the code just to support matrices.
function myfunm(a::Array{T,2}) where T
    out = Array{T}(undef, size(a)...)
    for k in 1:size(a,2)
        out[1,k] = 3 * (a[2,k] - a[1,k])
        out[2:end-1,k] = 3 * (a[3:end,k] - a[1:end-2,k])
        out[end,k] = 3 * (a[end,k] - a[end-1,k])
    end
    out
end

## It would be nice to preallocate the output and store the result directly there.
function myfunm2(a::Array{T,2}) where T
    hcat([myfun(a[:,k]) for k in 1:size(a, 2)]...)
end

## Mutable version of the same code
function myfun!(a::Array{T,1}) where T
    N = size(a, 1)
    mem, a[1] = a[1], 3 * (a[2] - a[1])
    for i in 2:N-1
        mem, a[i] = a[i], 3 * (a[i+1] - mem)
    end
    a[N] = 3 * (a[N] - mem)
end

## Another abomination
function myfunm!(a::Array{T,2}) where T
    N = size(a, 1)
    for k in 1:size(a, 2)
        mem, a[1,k] = a[1,k], 3 * (a[2,k] - a[1,k])
        for i in 2:N-1
            mem, a[i,k] = a[i,k], 3 * (a[i+1,k] - mem)
        end
        a[N,k] = 3 * (a[N,k] - mem)
    end
end

## This doesn't work
function myfunm2!(a::Array{T,2}) where T
    for k in 1:size(a, 2)
        myfun!(a[:,k])
    end
end

a = Array(1:5)
myfun(a)

aa = [1:5 11:15]
myfunm(aa) ≈ myfunm2(aa) #true

x = copy(a)
myfun!(x)
x ≈ myfun(a) # true

xx = copy(aa)
myfunm!(xx)
xx ≈ myfunm(aa) # true

xx = copy(aa)
myfunm2!(xx)
xx ≈ myfunm(aa) # false
14 
# This works!!... I'm happy again
function myfunm2!(a::Array{T,2}) where T
    for k in 1:size(a, 2)
        myfun!(@view a[:,k])
    end
end