# Assigning an array to a view

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>
``````

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

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?

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.

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.

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]``````

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

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

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?

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

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âŚ

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
``````
``````# 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
``````