The Julia Language Challenge

using BenchmarkTools

reference(out, a, b, c) = (out .= a .+ b .- sin.(c))

a = rand(1000, 1000);
b = rand(1000);
c = 1.0
out1 = similar(a);
out2 = similar(a);
br = @broadcast a + b - sin(c)

@btime materialize!($out1, $br)
@btime reference($out2, $a, $b, $c)
@assert out1 == out2
nothing

# Library A - I chose GeometryTypes. 
pkg"add GeometryTypes" # use package manager to install the package
# Any library with NVectors specializing to the length will do!
using GeometryTypes
# You need to take my word for it, that this is just a normal Julia struct
# Or you can inspect it with dump(Point3)
const Point3 = Point{3, Float32} 

# function needs to come from different library
module LibraryB
  # no odd stuff, no functors, no special lambda expression! 
  # this function needs to be a normal language function
  # as can be found in the wild
	super_custom_func(a, b) = sqrt(sum(a .* b))
end
# emulate that the function comes from a different library
using .LibraryB: super_custom_func 

using BenchmarkTools

a = rand(Point3, 10^6)
b = rand(Point3, 10^6)
out = fill(0f0, 10^6)

@btime $out .= super_custom_func.($a, $b)
br = @broadcast super_custom_func(a, b)
@btime materialize!($out, $br)
nothing

# Simplified implementation to take the axes of the array with the largest 
# dimensionality (axes -> the index range an array iterates over)
biggest(a, b, c, rest...) = biggest(biggest(a, b), biggest(c, rest...))
biggest(a::NTuple{N1, Any}, b::NTuple{N2, Any}) where {N1, N2} = 
	ifelse(N1 > N2, a, b)
biggest(a) = a
flatten_args(t::LazyBroadcast, rest...) = 
	(flatten_args(t.args...)..., flatten_args(rest...)...)
flatten_args(t::Any, rest...) = (t, flatten_args(rest...)...)
flatten_args() = ()
# the indexing axes of our array
axes(br::LazyBroadcast) = biggest(map(axes, flatten_args(br))...)
# lazy view that can be used to index all elements in br
eachindex(br::LazyBroadcast) = CartesianIndices(axes(br)) 

iterate(br::LazyBroadcast) = iterate(br, (eachindex(br),))
@propagate_inbounds function iterate(bc::LazyBroadcast, s)
    istate = iterate(s...) # iterate indices
    istate === nothing && return nothing # if done
    i, newstate = istate
    return (bc[i], (s[1], newstate))
end

@btime sum($br) 
@btime sum($out .= super_custom_func.($a, $b))

@assert sum(br) ≈ sum(super_custom_func.(a, b))

nothing

using Statistics
function threaded_materialize!(out::AbstractArray, x::LazyBroadcast)
    Threads.@threads for i in CartesianIndices(axes(out))
        @inbounds out[i] = x[i]
    end
    return out
end
a = @benchmark threaded_materialize!($out, $br)
b = @benchmark materialize!($out, $br)
hyper_threads = Sys.CPU_THREADS
# intels hyper threads don't really count in such intense calculations
cores = hyper_threads / 2
# scaling
println("best scaling ", minimum(b).time / minimum(a).time / cores)
println("average scaling ", mean(b).time / mean(a).time / cores)
nothing