Getting rid of BlockArrays

[nantonel]

That could simply code substantially.

Open to any idea!

[nantonel]

Using macros:

x,   y,  z = randn(5), randn(5), randn(5) # normal arrays
xb, yb, zb = (randn(5),randn(2,5)), (randn(5),randn(2,5)), (randn(5),randn(2,5)) # block arrays

macro block(ex1,ex2)

    if eval(:(eltype($ex1) <: AbstractArray))
           # if ex1 is Tuple of Arrays simply return the expression
           return esc(ex2)

       else eval(:( eltype($ex1) <: Tuple ))
           # if ex1 is Tuple of Tuple of Arrays construct
           # ( (x, y, z) -> x .= y .+ z ).(xb,yb,zb)
           return  esc( :( ( $ex1 -> $ex2 ).( $ex1... ) ) )
    end

end

using Test
# create new array
xx = @block (y,z) y .+ z
@test xx == y .+ z

# in place
@block (x,y,z) x .= y .+ z
@test x == y .+ z

# create new block array
zz = @block (yb,zb) yb .+ zb .+ 3
@test zz[1] == (yb[1] + zb[1] .+3)
@test zz[2] == (yb[2] + zb[2] .+3)

# in place
@block (xb,yb,zb) xb .= yb .+ 4 .*zb .+ 3
@test xb[1] == (yb[1] .+ 4 .*zb[1] .+ 3)
@test xb[2] == (yb[2] .+ 4 .*zb[2] .+ 3)

f(xb,yb,zb) = @block (xb,yb,zb) xb .= yb .+ 4 .*zb .+ 3
function g(xb,yb,zb)
    xb[1] .= yb[1] .+ 4 .*zb[1] .+ 3
    xb[2] .= yb[2] .+ 4 .*zb[2] .+ 3
end

using BenchmarkTools
println("\n\n with block arrays \n\n")
println("\n\n @block \n\n")
@btime f($xb,$yb,$zb)
println("\n\n non-block \n\n")
@btime g($xb,$yb,$zb)

f(x,y,z) = @block (x,y,z) x .= y .+ 4 .*z .+ 3
function g(x,y,z)
    x .= y .+ 4 .*z .+ 3
end

println("\n\n with arrays \n\n")
println("\n\n @block \n\n")
@btime f($x,$y,$z)
println("\n\n non-block \n\n")
@btime g($x,$y,$z)
return

[nantonel]

ok better macro, which avoids overhead (inspired by the base '@.' macro)

__block__(x, i::Int) = x
__block__(x::Symbol,i::Int) = Base.isoperator(x) || typeof(eval(x)) <: Number  ? x : Meta.parse(string(x)*"[$i]")
function __block__(x::Expr, i::Int)
    dotargs = map(__block__, x.args, i*ones(Int,length(x.args)))
    if x.head == :call && Base.isoperator(x.args[1])
        Expr(:., dotargs[1], Expr(:tuple, dotargs[2:end]...))
    else
        head = string(x.head)
        if last(head) == '=' && first(head) != '.'
            Expr(Symbol('.',head), dotargs...)
        else
            Expr(x.head, dotargs...)
        end
    end
end

macro block!(ex)
    if eval(:(typeof($(ex.args[1])) <: AbstractArray))
        # if ex begins with an Arrays add dots (same as @.)
        return esc(Base.Broadcast.__dot__(ex)) 
    else typeof(eval(ex.args[1])) <: NTuple{N,AbstractArray} where N
        # if ex begins with Tuple of Arrays construct
        # xb[1] .= yb[1] .+ 4 .*zb[1] .+ 3
        # xb[2] .= yb[2] .+ 4 .*zb[2] .+ 3
        new_ex = :()
        for i = 1:length(eval(ex.args[1]))
            blk = __block__(ex,i)
            new_ex = :( $(new_ex); $(blk))
        end
        return  esc( new_ex )
    end

end

using Test
x,   y,  z = randn(5), randn(5), randn(5) # normal arrays
xb = (randn(5),randn(2,5),randn(3)+im*randn(3)) # block arrays
yb = (randn(5),randn(2,5),randn(3)+im*randn(3)) # block arrays
zb = (randn(5),randn(2,5),randn(3)+im*randn(3)) # block arrays

# in place
@block! x = y + z
@test x == y .+ z

a = 4
## in place
@block! xb = yb + a *zb + 3.0
@test xb[1] == (yb[1] .+ a .*zb[1] .+ 3.0)
@test xb[2] == (yb[2] .+ a .*zb[2] .+ 3.0)
@test xb[3] == (yb[3] .+ a .*zb[3] .+ 3.0)
#
f(xb,yb,zb,a) = @block! xb = yb + a*zb + 3.0
function g(xb,yb,zb,a)  
    xb[1] .= yb[1] .+ a .*zb[1] .+ 3.0 
    xb[2] .= yb[2] .+ a .*zb[2] .+ 3.0
    xb[3] .= yb[3] .+ a .*zb[3] .+ 3.0
end

using BenchmarkTools
println("\n\n with block arrays \n\n")
println("\n\n @block! \n\n")
@btime f($xb,$yb,$zb,a) 
println("\n\n non-block \n\n")
@btime g($xb,$yb,$zb,a)

f(x,y,z,a) = @block! x = y + a *z + 3
function g(x,y,z,a)  
    x .= y .+ a .*z .+ 3
end

println("\n\n with arrays \n\n")
println("\n\n @block! \n\n")
@btime f($x,$y,$z,a)
println("\n\n non-block \n\n")
@btime g($x,$y,$z,a)
return

[lostella]

This one looks good! I'll now use it as a case-study to learn macros ;-)

[nantonel]

well probably there are simpler macros to begin with! 😆

[lostella]

Unfortunately, I think this will not work: doing evals and other runtime-related checks breaks the macro, which works in your example by coincidence I believe (because variables are defined in the same namespace of the macro, or something like that).

But not everything is lost! I’m working on fixing this! Merry Christmas!

[nantonel]

yes that macro proposal has surely some problems... but it's a modified version of the base @. macro! Possibly, here there will be potential solutions for this! Merry Xmas! 🎅

[lostella]

I think this has no easy solution in the realm of macros, at least I can see none. The only easy way of easily handling this type of "tuples of arrays" is to equip them with whatever makes AbstractArray easy to use, that is: iteration, indexing, broadcasting, and vector space operations (sum, scaling).

Luckily, someone has done this already, see the ArrayPartition <: AbstractArray type there: this is simply a wrapper type around tuples of arrays, with all necessary methods defined (similar, copy, copyto!, iterate, +, -, *, /, as well as all the indexing and broadcasting boilerplate) and can be used in place of other AbstractArray types. Give it a try, it is very nice and it is mantained, tested and versioned by the JuliaDiffEq people, which is also a plus. I mean, it is a dependency of DifferentialEquations.jl, so it has to be stable.

For example, I've added a test in PANOC.jl that runs the algorithm on a (fake) two-blocks problem using ArrayPartition as type for arrays. Note that in order to do that I did not touch any of the algorithm code, nor L-BFGS code: this simply relies on a couple of compatibility definitions to make prox and gradient work with ArrayPartition objects instead of Tuple. As long as writing algorithms is concerned, simply using ArrayPartition as type for arrays in this way seems to be enough. However, I had some issue with overriding the behaviour of AbstractOperator with this definition, due to some problem with the typing of methods in AbstractOperators.jl. I'll open a separate issue for this.

More importantly: If we feel like we could improve things internally (inside of AbstractOperators.jl and/or ProximalOperators.jl) we could directly switch to ArrayPartition inside of our operators.

I think these were my last Julia lines for 2018, happy new year! 🎉 🎉 🎉

[nantonel]

I’ll check this in detail next year! For the moment happy new year 🥳

[nantonel]

I think I'm quite convinced in adopting RecursiveArrayTools.jl. I'll try to make PR in the following days!