text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random

module TestBroadcastInternals

using Base.Broadcast: check_broadcast_axes, check_broadcast_shape, newindex, _bcs
using Base: OneTo
using Test, Random

@test @inferred(_bcs((3,5), (3,5))) == (3,5)
@test @inferred(_bcs((3,1), (3,5))) == (3,5)
@test @inferred(_bcs((3,),  (3,5))) == (3,5)
@test @inferred(_bcs((3,5), (3,)))  == (3,5)
@test_throws DimensionMismatch _bcs((3,5), (4,5))
@test_throws DimensionMismatch _bcs((3,5), (3,4))
@test @inferred(_bcs((-1:1, 2:5), (-1:1, 2:5))) == (-1:1, 2:5)
@test @inferred(_bcs((-1:1, 2:5), (1, 2:5)))    == (-1:1, 2:5)
@test @inferred(_bcs((-1:1, 1),   (1, 2:5)))    == (-1:1, 2:5)
@test @inferred(_bcs((-1:1,),     (-1:1, 2:5))) == (-1:1, 2:5)
@test_throws DimensionMismatch _bcs((-1:1, 2:6), (-1:1, 2:5))
@test_throws DimensionMismatch _bcs((-1:1, 2:5), (2, 2:5))

@test @inferred(Broadcast.combine_axes(zeros(3,4), zeros(3,4))) == (OneTo(3),OneTo(4))
@test @inferred(Broadcast.combine_axes(zeros(3,4), zeros(3)))   == (OneTo(3),OneTo(4))
@test @inferred(Broadcast.combine_axes(zeros(3),   zeros(3,4))) == (OneTo(3),OneTo(4))
@test @inferred(Broadcast.combine_axes(zeros(3), zeros(1,4), zeros(1))) == (OneTo(3),OneTo(4))

check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,5))
check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,1))
check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3))
check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,5), zeros(3))
check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,5), 1)
check_broadcast_axes((OneTo(3),OneTo(5)), 5, 2)
@test_throws DimensionMismatch check_broadcast_axes((OneTo(3),OneTo(5)), zeros(2,5))
@test_throws DimensionMismatch check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,4))
@test_throws DimensionMismatch check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,4,2))
@test_throws DimensionMismatch check_broadcast_axes((OneTo(3),OneTo(5)), zeros(3,5), zeros(2))
check_broadcast_axes((-1:1, 6:9), 1)

check_broadcast_shape((-1:1, 6:9), (-1:1, 6:9))
check_broadcast_shape((-1:1, 6:9), (-1:1, 1))
check_broadcast_shape((-1:1, 6:9), (1, 6:9))
@test_throws DimensionMismatch check_broadcast_shape((-1:1, 6:9), (-1, 6:9))
@test_throws DimensionMismatch check_broadcast_shape((-1:1, 6:9), (-1:1, 6))

ci(x) = CartesianIndex(x)
@test @inferred(newindex(ci((2,2)), (true, true), (-1,-1)))   == ci((2,2))
@test @inferred(newindex(ci((2,2)), (true, false), (-1,-1)))  == ci((2,-1))
@test @inferred(newindex(ci((2,2)), (false, true), (-1,-1)))  == ci((-1,2))
@test @inferred(newindex(ci((2,2)), (false, false), (-1,-1))) == ci((-1,-1))
@test @inferred(newindex(ci((2,2)), (true,), (-1,-1)))   == ci((2,))
@test @inferred(newindex(ci((2,2)), (true,), (-1,)))   == ci((2,))
@test @inferred(newindex(ci((2,2)), (false,), (-1,))) == ci((-1,))
@test @inferred(newindex(ci((2,2)), (), ())) == ci(())

end

function as_sub(x::AbstractVector)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = view(y, 2:2:length(y))
    y[:] = x[:]
    y
end
function as_sub(x::AbstractMatrix)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = view(y, 2:2:size(y,1), 2:2:size(y,2))
    for j=1:size(x,2)
        for i=1:size(x,1)
            y[i,j] = x[i,j]
        end
    end
    y
end
function as_sub(x::AbstractArray{T,3}) where T
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = view(y, 2:2:size(y,1), 2:2:size(y,2), 2:2:size(y,3))
    for k=1:size(x,3)
        for j=1:size(x,2)
            for i=1:size(x,1)
                y[i,j,k] = x[i,j,k]
            end
        end
    end
    y
end

bittest(f::Function, a...) = (@test f.(a...) == BitArray(broadcast(f, a...)))
n1 = 21
n2 = 32
n3 = 17
rb = 1:5

for arr in (identity, as_sub)
    @test broadcast(+, arr([1 0; 0 1]), arr([1, 4])) == [2 1; 4 5]
    @test broadcast(+, arr([1 0; 0 1]), arr([1  4])) == [2 4; 1 5]
    @test broadcast(+, arr([1  0]), arr([1, 4])) == [2 1; 5 4]
    @test broadcast(+, arr([1, 0]), arr([1  4])) == [2 5; 1 4]
    @test broadcast(+, arr([1, 0]), arr([1, 4])) == [2, 4]
    @test broadcast(+, arr([1, 0]), 2) == [3, 2]

    @test @inferred(broadcast(+, arr([1 0; 0 1]), arr([1, 4]))) == arr([2 1; 4 5])
    @test arr([1 0; 0 1]) .+ arr([1  4]) == arr([2 4; 1 5])
    @test arr([1  0]) .+ arr([1, 4]) == arr([2 1; 5 4])
    @test arr([1, 0]) .+ arr([1  4]) == arr([2 5; 1 4])
    @test arr([1, 0]) .+ arr([1, 4]) == arr([2, 4])
    @test arr([1]) .+ arr([]) == arr([])

    A = arr([1 0; 0 1]); @test broadcast!(+, A, A, arr([1, 4])) == arr([2 1; 4 5])
    A = arr([1 0; 0 1]); @test broadcast!(+, A, A, arr([1  4])) == arr([2 4; 1 5])
    A = arr([1  0]); @test_throws DimensionMismatch broadcast!(+, A, A, arr([1, 4]))
    A = arr([1  0]); @test broadcast!(+, A, A, arr([1  4])) == arr([2 4])
    A = arr([1  0]); @test broadcast!(+, A, A, 2) == arr([3 2])

    @test arr([ 1    2])   .* arr([3,   4])   == [ 3 6; 4 8]
    @test arr([24.0 12.0]) ./ arr([2.0, 3.0]) == [12 6; 8 4]
    @test arr([1 2]) ./ arr([3, 4]) == [1/3 2/3; 1/4 2/4]
    @test arr([1 2]) .\ arr([3, 4]) == [3 1.5; 4 2]
    @test arr([3 4]) .^ arr([1, 2]) == [3 4; 9 16]
    @test arr(BitArray([true false])) .* arr(BitArray([true, true])) == [true false; true false]
    @test arr(BitArray([true false])) .^ arr(BitArray([false, true])) == [true true; true false]
    @test arr(BitArray([true false])) .^ arr([0, 3]) == [true true; true false]

    M = arr([11 12; 21 22])
    @test getindex.((M,), [2 1; 1 2], arr([1, 2])) == [21 11; 12 22]
    @test_throws BoundsError getindex.((M,), [2 1; 1 2], arr([1, -1]))
    @test_throws BoundsError getindex.((M,), [2 1; 1 2], arr([1, 2]), [2])
    @test getindex.((M,), [2 1; 1 2],arr([2, 1]), [1]) == [22 12; 11 21]

    A = arr(zeros(2,2))
    setindex!.((A,), arr([21 11; 12 22]), [2 1; 1 2], arr([1, 2]))
    @test A == M
    setindex!.((A,), 5, [1,2], [2 2])
    @test A == [11 5; 21 5]
    setindex!.((A,), 7, [1,2], [1 2])
    @test A == fill(7, 2, 2)
    A = arr(zeros(3,3))
    setindex!.((A,), 10:12, 1:3, 1:3)
    @test A == [10 0 0; 0 11 0; 0 0 12]
    @test_throws BoundsError setindex!.((A,), 7, [1,-1], [1 2])

    for f in ((==), (<) , (!=), (<=))
        bittest(f, arr([1 0; 0 1]), arr([1, 4]))
        bittest(f, arr([1 0; 0 1]), arr([1  4]))
        bittest(f, arr([0, 1]), arr([1  4]))
        bittest(f, arr([0  1]), arr([1, 4]))
        bittest(f, arr([1, 0]), arr([1, 4]))
        bittest(f, arr(rand(rb, n1, n2, n3)), arr(rand(rb, n1, n2, n3)))
        bittest(f, arr(rand(rb,  1, n2, n3)), arr(rand(rb, n1,  1, n3)))
        bittest(f, arr(rand(rb,  1, n2,  1)), arr(rand(rb, n1,  1, n3)))
        bittest(f, arr(bitrand(n1, n2, n3)), arr(bitrand(n1, n2, n3)))
    end
end

let r1 = 1:1,
    r2 = 1:5,
    ratio = [1,1/2,1/3,1/4,1/5],
    m = [1:2;]'
    @test r1.*r2 == [1:5;]
    @test r1./r2 == ratio
    @test m.*r2 == [1:5 2:2:10]
    @test m./r2 ≈ [ratio 2ratio]
    @test m./[r2;] ≈ [ratio 2ratio]
end

@test @inferred(broadcast(+,[0,1.2],reshape([0,-2],1,1,2))) == reshape([0 -2; 1.2 -0.8],2,1,2)
rt = Base.return_types(broadcast, Tuple{typeof(+), Array{Float64, 3}, Array{Int, 1}})
@test length(rt) == 1 && rt[1] == Array{Float64, 3}
rt = Base.return_types(broadcast!, Tuple{Function, Array{Float64, 3}, Array{Float64, 3}, Array{Int, 1}})
@test length(rt) == 1 && rt[1] == Array{Float64, 3}

# f.(args...) syntax (#15032)
let x = [1, 3.2, 4.7],
    y = [3.5, pi, 1e-4],
    α = 0.2342
    @test sin.(x) == broadcast(sin, x)
    @test sin.(α) == broadcast(sin, α)
    @test sin.(3.2) == broadcast(sin, 3.2) == sin(3.2)
    @test factorial.(3) == broadcast(factorial, 3)
    @test atan.(x, y) == broadcast(atan, x, y)
    @test atan.(x, y') == broadcast(atan, x, y')
    @test atan.(x, α) == broadcast(atan, x, α)
    @test atan.(α, y') == broadcast(atan, α, y')
end

# issue 14725
let a = Number[2, 2.0, 4//2, 2+0im] / 2
    @test eltype(a) == Number
end
let a = Real[2, 2.0, 4//2] / 2
    @test eltype(a) == Real
end
let a = Real[2, 2.0, 4//2] / 2.0
    @test eltype(a) == Float64
end

# issue 16164
let a = broadcast(Float32, [3, 4, 5])
    @test eltype(a) == Float32
end

# broadcasting scalars:
@test sin.(1) === broadcast(sin, 1) === sin(1)
@test (()->1234).() === broadcast(()->1234) === 1234

# issue #4883
@test isa(broadcast(tuple, [1 2 3], ["a", "b", "c"]), Matrix{Tuple{Int,String}})
@test isa(broadcast((x,y)->(x==1 ? 1.0 : x, y), [1 2 3], ["a", "b", "c"]), Matrix{Tuple{Real,String}})
let a = length.(["foo", "bar"])
    @test isa(a, Vector{Int})
    @test a == [3, 3]
end
let a = sin.([1, 2])
    @test isa(a, Vector{Float64})
    @test a ≈ [0.8414709848078965, 0.9092974268256817]
end

# PR #17300: loop fusion
@test (x->x+1).((x->x+2).((x->x+3).(1:10))) == 7:16
let A = [sqrt(i)+j for i = 1:3, j=1:4]
    @test atan.(log.(A), sum(A, dims=1)) == broadcast(atan, broadcast(log, A), sum(A, dims=1))
end
let x = sin.(1:10)
    @test atan.((x->x+1).(x), (x->x+2).(x)) == broadcast(atan, x.+1, x.+2)
    @test sin.(atan.([x.+1,x.+2]...)) == sin.(atan.(x.+1 ,x.+2)) == @. sin(atan(x+1,x+2))
    @test sin.(atan.(x, 3.7)) == broadcast(x -> sin(atan(x,3.7)), x)
    @test atan.(x, 3.7) == broadcast(x -> atan(x,3.7), x) == broadcast(atan, x, 3.7)
end
# Use side effects to check for loop fusion.
let g = Int[]
    f17300(x) = begin; push!(g, x); x+2; end
    f17300.(f17300.(f17300.(1:3)))
    @test g == [1,3,5, 2,4,6, 3,5,7]
    empty!(g)
    @. f17300(f17300(f17300(1:3)))
    @test g == [1,3,5, 2,4,6, 3,5,7]
end
# fusion with splatted args:
let x = sin.(1:10), a = [x]
    @test cos.(x) == cos.(a...)
    @test atan.(x,x) == atan.(a..., a...) == atan.([x, x]...)
    @test atan.(x, cos.(x)) == atan.(a..., cos.(x)) == broadcast(atan, x, cos.(a...)) == broadcast(atan, a..., cos.(a...))
    @test ((args...)->cos(args[1])).(x) == cos.(x) == ((y,args...)->cos(y)).(x)
end
@test atan.(3, 4) == atan(3, 4) == (() -> atan(3, 4)).()
# fusion with keyword args:
let x = [1:4;]
    f17300kw(x; y=0) = x + y
    @test f17300kw.(x) == x
    @test f17300kw.(x, y=1) == f17300kw.(x; y=1) == f17300kw.(x; [(:y,1)]...) == x .+ 1 == [2, 3, 4, 5]
    @test f17300kw.(sin.(x), y=1) == f17300kw.(sin.(x); y=1) == sin.(x) .+ 1
    @test sin.(f17300kw.(x, y=1)) == sin.(f17300kw.(x; y=1)) == sin.(x .+ 1)
end

# issue #23236
let X = [[true,false],[false,true]]
    @test [.!x for x in X] == [[false,true],[true,false]]
end

# splice escaping of @.
let x = [4, -9, 1, -16]
    @test [2, 3, 4, 5] == @.(1 + sqrt($sort(abs(x))))
end

# interaction of @. with let
@test [1,4,9] == @. let x = [1,2,3]; x^2; end

# interaction of @. with for loops
let x = [1,2,3], y = x
    @. for i = 1:3
        y = y^2 # should convert to y .= y.^2
    end
    @test x == [1,256,6561]
end

# interaction of @. with function definitions
let x = [1,2,3]
    @. f(x) = x^2
    @test f(x) == [1,4,9]
end

# Issue #23622: @. with chained comparisons
let x = [1,2,3]
    @test (1 .< x .< 3) == @.(1 < x < 3) == (@. 1 .< x .< 3) == [false, true, false]
    @test (x .=== 1:3 .=== [1,2,3]) == @.(x === 1:3 === [1,2,3]) == [true, true, true]
end

# PR #17510: Fused in-place assignment
let x = [1:4;], y = x
    y .= 2:5
    @test y === x == [2:5;]
    y .= factorial.(x)
    @test y === x == [2,6,24,120]
    y .= 7
    @test y === x == [7,7,7,7]
    y .= factorial.(3)
    @test y === x == [6,6,6,6]
    f17510() = 9
    y .= f17510.()
    @test y === x == [9,9,9,9]
    y .-= 1
    @test y === x == [8,8,8,8]
    @. y -= 1:4          # @. should convert to .-=
    @test y === x == [7,6,5,4]
    x[1:2] .= 1
    @test y === x == [1,1,5,4]
    @. x[1:2] .+= [2,3]  # use .+= to make sure @. works with dotted assignment
    @test y === x == [3,4,5,4]
    @. x[:] .= 0         # use .= to make sure @. works with dotted assignment
    @test y === x == [0,0,0,0]
    @. x[2:end] = 1:3    # @. should convert to .=
    @test y === x == [0,1,2,3]
end
let a = [[4, 5], [6, 7]], b = reshape(a, 1, 2)
    a[1] .= 3
    @test a == [[3, 3], [6, 7]]
    a[CartesianIndex(1)] .= 4
    @test a == [[4, 4], [6, 7]]
    b[1, CartesianIndex(1)] .= 5
    @test a == [[5, 5], [6, 7]]
end
let d = Dict(:foo => [1,3,7], (3,4) => [5,9])
    d[:foo] .+= 2
    @test d[:foo] == [3,5,9]
    d[3,4] .-= 1
    @test d[3,4] == [4,8]
end
let identity = error, x = [1,2,3]
    x .= 1 # make sure it goes to broadcast!(Base.identity, ...), not identity
    @test x == [1,1,1]
end

# make sure scalars are inlined, which causes f.(x,scalar) to lower to a "thunk"
import Base.Meta: isexpr
@test isexpr(Meta.lower(Main, :(f.(x,1))), :thunk)
@test isexpr(Meta.lower(Main, :(f.(x,1.0))), :thunk)
@test isexpr(Meta.lower(Main, :(f.(x,$π))), :thunk)
@test isexpr(Meta.lower(Main, :(f.(x,"hello"))), :thunk)
@test isexpr(Meta.lower(Main, :(f.(x,$("hello")))), :thunk)

# PR #17623: Fused binary operators
@test [true] .* [true] == [true]
@test [1,2,3] .|> (x->x+1) == [2,3,4]
let g = Int[], ⊕ = (a,b) -> let c=a+2b; push!(g, c); c; end
    @test [1,2,3] .⊕ [10,11,12] .⊕ [100,200,300] == [221,424,627]
    @test g == [21,221,24,424,27,627] # test for loop fusion
end

# Fused unary operators
@test .√[3,4,5] == sqrt.([3,4,5])
@test .![true, true, false] == [false, false, true]
@test .-[1,2,3] == -[1,2,3] == .+[-1,-2,-3] == [-1,-2,-3]

# PR 16988
@test Base.promote_op(+, Bool) === Int
@test isa(broadcast(+, [true]), Array{Int,1})

# issue #17304
let foo = [[1,2,3],[4,5,6],[7,8,9]]
    @test max.(foo...) == broadcast(max, foo...) == [7,8,9]
end

# Issue 17314
@test broadcast(x->log(log(log(x))), [1000]) == [log(log(log(1000)))]
let f17314 = x -> x < 0 ? false : x
    @test eltype(broadcast(f17314, 1:3)) === Int
    @test eltype(broadcast(f17314, -1:1)) === Integer
    @test eltype(broadcast(f17314, Int[])) === Integer
end
let io = IOBuffer()
    broadcast(x->print(io,x), 1:5) # broadcast with side effects
    @test take!(io) == [0x31,0x32,0x33,0x34,0x35]
end

# Issue 18176
let f18176(a, b, c) = a + b + c
    @test f18176.(1.0:2, 3, 4) == f18176.(3.0, 1.0:2, 4.0) == broadcast(f18176, 3, 4, 1.0:2)
end

# Issue #17984
let A17984 = []
    @test isa(abs.(A17984), Array{Any,1})
end

# Issue #16966
@test parse.(Int, "1") == 1
@test parse.(Int, ["1", "2"]) == [1, 2]
@test trunc.((Int,), [1.2, 3.4]) == [1, 3]
@test abs.((1, -2)) == (1, 2)
@test broadcast(+, 1.0, (0, -2.0)) == (1.0,-1.0)
@test broadcast(+, 1.0, (0, -2.0), [1]) == [2.0, 0.0]
@test broadcast(*, ["Hello"], ", ", ["World"], "!") == ["Hello, World!"]
let s = "foo"
    @test s .* ["bar", "baz"] == ["foobar", "foobaz"] == "foo" .* ["bar", "baz"]
end

# Ensure that even strange constructors that break `T(x)::T` work with broadcast
struct StrangeType18623 end
StrangeType18623(x) = x
StrangeType18623(x,y) = (x,y)
@test @inferred(broadcast(StrangeType18623, 1:3)) == [1,2,3]
@test @inferred(broadcast(StrangeType18623, 1:3, 4:6)) == [(1,4),(2,5),(3,6)]

@test typeof(Int.(Number[1, 2, 3])) === typeof((x->Int(x)).(Number[1, 2, 3]))

@test @inferred(broadcast(CartesianIndex, 1:2)) == [CartesianIndex(1), CartesianIndex(2)]
@test @inferred(broadcast(CartesianIndex, 1:2, 3:4)) == [CartesianIndex(1,3), CartesianIndex(2,4)]

# Issue 18622
@test @inferred(broadcast(muladd, [1.0], [2.0], [3.0])) == [5.0]
@test @inferred(broadcast(tuple, 1:3, 4:6, 7:9)) == [(1,4,7), (2,5,8), (3,6,9)]

# 19419
@test @inferred(broadcast(round, Int, [1])) == [1]

# https://discourse.julialang.org/t/towards-broadcast-over-combinations-of-sparse-matrices-and-scalars/910
let
    f(A, n) = broadcast(x -> +(x, n), A)
    @test @inferred(f([1.0], 1)) == [2.0]
    g() = (a = 1; Broadcast.combine_eltypes(x -> x + a, (1.0,)))
    @test @inferred(g()) === Float64
end

# Ref as 0-dimensional array for broadcast
@test (-).(C_NULL, C_NULL)::UInt == 0
@test (+).(1, Ref(2)) == 3
@test (+).(Ref(1), Ref(2)) == 3
@test (+).([[0,2], [1,3]], Ref{Vector{Int}}([1,-1])) == [[1,1], [2,2]]

# Check that broadcast!(f, A) populates A via independent calls to f (#12277, #19722),
# and similarly for broadcast!(f, A, numbers...) (#19799).
@test let z = 1; A = broadcast!(() -> z += 1, zeros(2)); A[1] != A[2]; end
@test let z = 1; A = broadcast!(x -> z += x, zeros(2), 1); A[1] != A[2]; end

## broadcasting for custom AbstractArray
abstract type ArrayData{T,N} <: AbstractArray{T,N} end
Base.getindex(A::ArrayData, i::Integer...) = A.data[i...]
Base.setindex!(A::ArrayData, v::Any, i::Integer...) = setindex!(A.data, v, i...)
Base.size(A::ArrayData) = size(A.data)
Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{A}}, ::Type{T}) where {A,T} =
    A(Array{T}(undef, size(bc)))

struct Array19745{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:Array19745} = Broadcast.ArrayStyle{Array19745}()

# Two specialized broadcast rules with no declared precedence
struct AD1{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD1} = Broadcast.ArrayStyle{AD1}()
struct AD2{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD2} = Broadcast.ArrayStyle{AD2}()

# Two specialized broadcast rules with explicit precedence
struct AD1P{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD1P} = Broadcast.ArrayStyle{AD1P}()
struct AD2P{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD2P} = Broadcast.ArrayStyle{AD2P}()

Base.BroadcastStyle(a1::Broadcast.ArrayStyle{AD1P}, ::Broadcast.ArrayStyle{AD2P}) = a1

# Two specialized broadcast rules where users unnecessarily
# define `BroadcastStyle` for both argument orders (but do so consistently)
struct AD1B{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD1B} = Broadcast.ArrayStyle{AD1B}()
struct AD2B{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD2B} = Broadcast.ArrayStyle{AD2B}()

Base.BroadcastStyle(a1::Broadcast.ArrayStyle{AD1B}, a2::Broadcast.ArrayStyle{AD2B}) = a1
Base.BroadcastStyle(a2::Broadcast.ArrayStyle{AD2B}, a1::Broadcast.ArrayStyle{AD1B}) = a1

# Two specialized broadcast rules with conflicting precedence
struct AD1C{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD1C} = Broadcast.ArrayStyle{AD1C}()
struct AD2C{T,N} <: ArrayData{T,N}
    data::Array{T,N}
end
Base.BroadcastStyle(::Type{T}) where {T<:AD2C} = Broadcast.ArrayStyle{AD2C}()

Base.BroadcastStyle(a1::Broadcast.ArrayStyle{AD1C}, a2::Broadcast.ArrayStyle{AD2C}) = a1
Base.BroadcastStyle(a2::Broadcast.ArrayStyle{AD2C}, a1::Broadcast.ArrayStyle{AD1C}) = a2

# A Custom type with specific dimensionality
struct AD2Dim{T} <: ArrayData{T,2}
    data::Array{T,2}
end
struct AD2DimStyle <: Broadcast.AbstractArrayStyle{2}; end
AD2DimStyle(::Val{2}) = AD2DimStyle()
AD2DimStyle(::Val{N}) where {N} = Broadcast.DefaultArrayStyle{N}()
Base.similar(bc::Broadcast.Broadcasted{AD2DimStyle}, ::Type{T}) where {T} =
    AD2Dim(Array{T}(undef, size(bc)))
Base.BroadcastStyle(::Type{T}) where {T<:AD2Dim} = AD2DimStyle()

@testset "broadcasting for custom AbstractArray" begin
    a  = randn(10)
    aa = Array19745(a)
    fadd(aa) = aa .+ 1
    fadd2(aa) = aa .+ 1 .* 2
    fadd3(aa) = aa .+ [missing; 1:9]
    fprod(aa) = aa .* aa'
    @test a .+ 1  == @inferred(fadd(aa))
    @test a .+ 1 .* 2  == @inferred(fadd2(aa))
    @test a .* a' == @inferred(fprod(aa))
    @test isequal(a .+ [missing; 1:9], fadd3(aa))
    @test Core.Compiler.return_type(fadd3, (typeof(aa),)) <: Array19745{<:Union{Float64, Missing}}
    @test isa(aa .+ 1, Array19745)
    @test isa(aa .+ 1 .* 2, Array19745)
    @test isa(aa .* aa', Array19745)
    @test isa(aa .* [missing; 1:9], Array19745)
    a1 = AD1(rand(2,3))
    a2 = AD2(rand(2))
    @test a1 .+ 1 isa AD1
    @test a2 .+ 1 isa AD2
    @test a1 .+ 1 .* 2 isa AD1
    @test a2 .+ 1 .* 2 isa AD2
    @test a1 .+ a2 isa Array
    @test a2 .+ a1 isa Array
    @test a1 .+ a2 .+ a1 isa Array
    @test a1 .+ a2 .+ a2 isa Array
    a1 = AD1P(rand(2,3))
    a2 = AD2P(rand(2))
    @test a1 .+ 1 isa AD1P
    @test a2 .+ 1 isa AD2P
    @test a1 .+ 1 .* 2 isa AD1P
    @test a2 .+ 1 .* 2 isa AD2P
    @test a1 .+ a2 isa AD1P
    @test a2 .+ a1 isa AD1P
    @test a1 .+ a2 .+ a1 isa AD1P
    @test a1 .+ a2 .+ a2 isa AD1P
    a1 = AD1B(rand(2,3))
    a2 = AD2B(rand(2))
    @test a1 .+ 1 isa AD1B
    @test a2 .+ 1 isa AD2B
    @test a1 .+ 1 .* 2 isa AD1B
    @test a2 .+ 1 .* 2 isa AD2B
    @test a1 .+ a2 isa AD1B
    @test a2 .+ a1 isa AD1B
    @test a1 .+ a2 .+ a1 isa AD1B
    @test a1 .+ a2 .+ a2 isa AD1B
    a1 = AD1C(rand(2,3))
    a2 = AD2C(rand(2))
    @test a1 .+ 1 isa AD1C
    @test a2 .+ 1 isa AD2C
    @test a1 .+ 1 .* 2 isa AD1C
    @test a2 .+ 1 .* 2 isa AD2C
    @test_throws ErrorException a1 .+ a2
    a2d = AD2Dim(rand(2, 3))
    a2 = AD2(rand(2))
    @test a2d .+ 1 isa AD2Dim
    @test a2d .+ a2 isa Matrix
    @test a2d .+ (1:2) isa AD2Dim
    @test a2d .+ ones(2, 3) isa AD2Dim
    @test a2d .+ ones(2, 3, 4) isa Array{Float64, 3}
end

# broadcast should only "peel off" one container layer
@test getindex.([Ref(1), Ref(2)]) == [1, 2]
let io = IOBuffer()
    broadcast(x -> print(io, x), [Ref(1.0)])
    @test String(take!(io)) == "Base.RefValue{Float64}(1.0)"
end

# Test that broadcast's promotion mechanism handles closures accepting more than one argument.
# (See issue #19641 and referenced issues and pull requests.)
let f() = (a = 1; Broadcast.combine_eltypes((x, y) -> x + y + a, (1.0, 1.0)))
    @test @inferred(f()) == Float64
end

@testset "broadcast resulting in BitArray" begin
    let f(x) = x ? true : "false"
        ba = f.([true])
        @test ba isa BitArray
        @test ba == [true]
        a = f.([false])
        @test a isa Array{String}
        @test a == ["false"]
        @test f.([true, false]) == [true, "false"]
    end
end

# Test that broadcast treats type arguments as scalars, i.e. containertype yields Any,
# even for subtypes of abstract array. (https://github.com/JuliaStats/DataArrays.jl/issues/229)
@testset "treat type arguments as scalars, DataArrays issue 229" begin
    @test Broadcast.combine_styles(Broadcast.broadcastable(AbstractArray)) == Base.Broadcast.DefaultArrayStyle{0}()
    @test broadcast(==, [1], AbstractArray) == BitArray([false])
    @test broadcast(==, 1, AbstractArray) == false
end

@testset "broadcasting falls back to iteration (issues #26421, #19577, #23746)" begin
    @test_throws ArgumentError broadcast(identity, Dict(1=>2))
    @test_throws ArgumentError broadcast(identity, (a=1, b=2))
    @test_throws ArgumentError length.(Dict(1 => BitSet(1:2), 2 => BitSet(1:3)))
    @test_throws MethodError broadcast(identity, Base)

    @test broadcast(identity, Iterators.filter(iseven, 1:10)) == 2:2:10
    d = Dict([1,2] => 1.1, [3,2] => 0.1)
    @test length.(keys(d)) == [2,2]
    @test Set(exp.(Set([1,2,3]))) == Set(exp.([1,2,3]))
end

# Test that broadcasting identity where the input and output Array shapes do not match
# yields the correct result, not merely a partial copy. See pull request #19895 for discussion.
let N = 5
    @test iszero(fill(1, N, N) .= zeros(N, N))
    @test iszero(fill(1, N, N) .= zeros(N, 1))
    @test iszero(fill(1, N, N) .= zeros(1, N))
    @test iszero(fill(1, N, N) .= zeros(1, 1))
end

@testset "test broadcast for matrix of matrices" begin
    A = fill([0 0; 0 0], 4, 4)
    A[1:3,1:3] .= [[1 1; 1 1]]
    @test all(A[1:3,1:3] .== [[1 1; 1 1]])
end

# Test that broadcast does not confuse eltypes. See also
# https://github.com/JuliaLang/julia/issues/21325
@testset "eltype confusion (#21325)" begin
    foo(x::Char, y::Int) = 0
    foo(x::String, y::Int) = "hello"
    @test broadcast(foo, "x", [1, 2, 3]) == ["hello", "hello", "hello"]

    @test isequal(
        [Set([1]), Set([2])] .∪ Ref(Set([3])),
        [Set([1, 3]), Set([2, 3])])
end

# A bare bones custom type that supports broadcast
struct Foo26601{T}
    data::T
end
Base.axes(f::Foo26601) = axes(f.data)
Base.getindex(f::Foo26601, i...) = getindex(f.data, i...)
Base.ndims(::Type{Foo26601{T}}) where {T} = ndims(T)
Base.Broadcast.broadcastable(f::Foo26601) = f
@testset "barebones custom object broadcasting" begin
    for d in (rand(Float64, ()), rand(5), rand(5,5), rand(5,5,5))
        f = Foo26601(d)
        @test f .* 2 == d .* 2
        @test f .* (1:5) == d .* (1:5)
        @test f .* reshape(1:25,5,5) == d .* reshape(1:25,5,5)
        @test sqrt.(f) == sqrt.(d)
        @test f .* (1,2,3,4,5) == d .* (1,2,3,4,5)
    end
end

@testset "broadcast resulting in tuples" begin
    # Issue #21291
    let t = (0, 1, 2)
        o = 1
        @test @inferred(broadcast(+, t, o)) == (1, 2, 3)
    end

    # Issue #23647
    @test (1, 2, 3) .+ (1,) == (1,) .+ (1, 2, 3) == (2, 3, 4)
    @test (1,) .+ () == () .+ (1,) == () .+ () == ()
    @test (1, 2) .+ (1, 2) == (2, 4)
    @test_throws DimensionMismatch (1, 2) .+ (1, 2, 3)
end

@testset "broadcasted assignment from tuples and tuple styles (#33020)" begin
    a = zeros(3)
    @test_throws DimensionMismatch a .= (1,2)
    @test_throws DimensionMismatch a .= sqrt.((1,2))
    a .= (1,)
    @test all(==(1), a)
    a .= sqrt.((2,))
    @test all(==(√2), a)
    a = zeros(3, 2)
    @test_throws DimensionMismatch a .= (1,2)
    @test_throws DimensionMismatch a .= sqrt.((1,2))
    a .= (1,)
    @test all(==(1), a)
    a .= sqrt.((2,))
    @test all(==(√2), a)
    a .= (1,2,3)
    @test a == [1 1; 2 2; 3 3]
end

@testset "scalar .=" begin
    A = [[1,2,3],4:5,6]
    A[1] .= 0
    @test A[1] == [0,0,0]
    @test_throws ErrorException A[2] .= 0
    @test_throws MethodError A[3] .= 0
    A = [[1,2,3],4:5]
    A[1] .= 0
    @test A[1] == [0,0,0]
    @test_throws ErrorException A[2] .= 0
end

# Issue #22180
@test convert.(Any, [1, 2]) == [1, 2]

# Issue #24944
let n = 1
    @test ceil.(Int, n ./ (1,)) == (1,)
    @test ceil.(Int, 1 ./ (1,)) == (1,)
end

# Issue #29266
@testset "deprecated scalar-fill .=" begin
    a = fill(1, 10)
    @test_throws ArgumentError a[1:5] = 0

    x = randn(10)
    @test_throws ArgumentError x[x .> 0.0] = 0.0
end

# lots of splatting!
let x = [[1, 4], [2, 5], [3, 6]]
    y = .+(x..., .*(x..., x...)..., x[1]..., x[2]..., x[3]...)
    @test y == [14463, 14472]

    z = zeros(2)
    z .= .+(x..., .*(x..., x...)..., x[1]..., x[2]..., x[3]...)
    @test z == Float64[14463, 14472]
end

# Issue #21094
@generated function foo21094(out, x)
    quote
        out .= x .+ x
        out
    end
end
@test foo21094([0.0], [1.0]) == [2.0]

# Issue #22053
struct T22053
    t
end
Broadcast.BroadcastStyle(::Type{T22053}) = Broadcast.Style{T22053}()
Broadcast.axes(::T22053) = ()
Broadcast.broadcastable(t::T22053) = t
function Base.copy(bc::Broadcast.Broadcasted{Broadcast.Style{T22053}})
    all(x->isa(x, T22053), bc.args) && return 1
    return 0
end
Base.:*(::T22053, ::T22053) = 2
let x = T22053(1)
    @test x*x == 2
    @test x.*x == 1
end

# Issue https://github.com/JuliaLang/julia/pull/25377#discussion_r159956996
let X = Any[1,2]
    X .= nothing
    @test X[1] == X[2] == nothing
end

# Ensure that broadcast styles with custom indexing work
let X = zeros(2, 3)
    X .= (1, 2)
    @test X == [1 1 1; 2 2 2]
end

# issue #27988: inference of Broadcast.flatten
using .Broadcast: Broadcasted
let
    bc = Broadcasted(+, (Broadcasted(*, (1, 2)), Broadcasted(*, (Broadcasted(*, (3, 4)), 5))))
    @test @inferred(Broadcast.cat_nested(bc)) == (1,2,3,4,5)
    @test @inferred(Broadcast.materialize(Broadcast.flatten(bc))) == @inferred(Broadcast.materialize(bc)) == 62
    bc = Broadcasted(+, (Broadcasted(*, (1, Broadcasted(/, (2.0, 2.5)))), Broadcasted(*, (Broadcasted(*, (3, 4)), 5))))
    @test @inferred(Broadcast.cat_nested(bc)) == (1,2.0,2.5,3,4,5)
    @test @inferred(Broadcast.materialize(Broadcast.flatten(bc))) == @inferred(Broadcast.materialize(bc)) == 60.8
end

let
  bc = Broadcasted(+, (Broadcasted(*, ([1, 2, 3], 4)), 5))
  @test isbits(Broadcast.flatten(bc).f)
end

# Issue #26127: multiple splats in a fused dot-expression
let f(args...) = *(args...)
    x, y, z = (1,2), 3, (4, 5)
    @test f.(x..., y, z...) == broadcast(f, x..., y, z...) == 120
    @test f.(x..., f.(x..., y, z...), y, z...) == broadcast(f, x..., broadcast(f, x..., y, z...), y, z...) == 120*120
end

@testset "Issue #27911: Broadcasting over collections with big indices" begin
    @test iszero.(Int128(0):Int128(2)) == [true, false, false]
    @test iszero.((Int128(0):Int128(2)) .- 1) == [false, true, false]
    @test iszero.(big(0):big(2)) == [true, false, false]
    @test iszero.((big(0):big(2)) .- 1) == [false, true, false]
end

@testset "Issue #27775: Broadcast!ing over nested scalar operations" begin
    a = zeros(2)
    a .= 1 ./ (1 + 2)
    @test a == [1/3, 1/3]
    a .= 1 ./ (1 .+ 3)
    @test a == [1/4, 1/4]
    a .= sqrt.(1 ./ 2)
    @test a == [sqrt(1/2), sqrt(1/2)]
    rng = MersenneTwister(1234)
    a .= rand.((rng,))
    rng = MersenneTwister(1234)
    @test a == [rand(rng), rand(rng)]
    @test a[1] != a[2]
    rng = MersenneTwister(1234)
    broadcast!(rand, a, (rng,))
    rng = MersenneTwister(1234)
    @test a == [rand(rng), rand(rng)]
    @test a[1] != a[2]
end

# Issue #27446: Broadcasting pair operator
let
    c = ["foo", "bar"]
    d = [1,2]
    @test Dict(c .=> d) == Dict("foo" => 1, "bar" => 2)
end

# Broadcasted iterable/indexable APIs
let
    bc = Broadcast.instantiate(Broadcast.broadcasted(+, zeros(5), 5))
    @test IndexStyle(bc) == IndexLinear()
    @test eachindex(bc) === Base.OneTo(5)
    @test length(bc) === 5
    @test ndims(bc) === 1
    @test ndims(typeof(bc)) === 1
    @test bc[1] === bc[CartesianIndex((1,))] === 5.0
    @test copy(bc) == [v for v in bc] == collect(bc)
    @test eltype(copy(bc)) == eltype([v for v in bc]) == eltype(collect(bc))
    @test ndims(copy(bc)) == ndims([v for v in bc]) == ndims(collect(bc)) == ndims(bc)

    bc = Broadcast.instantiate(Broadcast.broadcasted(+, zeros(5), 5*ones(1, 4)))
    @test IndexStyle(bc) == IndexCartesian()
    @test eachindex(bc) === CartesianIndices((Base.OneTo(5), Base.OneTo(4)))
    @test length(bc) === 20
    @test ndims(bc) === 2
    @test ndims(typeof(bc)) === 2
    @test bc[1,1] == bc[CartesianIndex((1,1))] === 5.0
    @test copy(bc) == [v for v in bc] == collect(bc)
    @test eltype(copy(bc)) == eltype([v for v in bc]) == eltype(collect(bc))
    @test ndims(copy(bc)) == ndims([v for v in bc]) == ndims(collect(bc)) == ndims(bc)
end

# issue #31295
let a = rand(5), b = rand(5), c = copy(a)
    view(identity(a), 1:3) .+= view(b, 1:3)
    @test a == [(c+b)[1:3]; c[4:5]]

    x = [1]
    x[[1,1]] .+= 1
    @test x == [2]
end

@testset "broadcasted mapreduce" begin
    xs = 1:10
    ys = 1:2:20
    bc = Broadcast.instantiate(Broadcast.broadcasted(*, xs, ys))
    @test IndexStyle(bc) == IndexLinear()
    @test sum(bc) == mapreduce(Base.splat(*), +, zip(xs, ys))

    xs2 = reshape(xs, 1, :)
    ys2 = reshape(ys, 1, :)
    bc = Broadcast.instantiate(Broadcast.broadcasted(*, xs2, ys2))
    @test IndexStyle(bc) == IndexCartesian()
    @test sum(bc) == mapreduce(Base.splat(*), +, zip(xs, ys))

    xs = 1:5:3*5
    ys = 1:4:3*4
    bc = Broadcast.instantiate(
        Broadcast.broadcasted(iseven, Broadcast.broadcasted(-, xs, ys)))
    @test count(bc) == count(iseven, map(-, xs, ys))

    xs = reshape(1:6, (2, 3))
    ys = 1:2
    bc = Broadcast.instantiate(Broadcast.broadcasted(*, xs, ys))
    @test reduce(+, bc; dims=1, init=0) == [5 11 17]

    # Let's test that `Broadcasted` actually hits the efficient
    # `mapreduce` method as intended.  We are going to invoke `reduce`
    # with this *NON-ASSOCIATIVE* binary operator to see what
    # associativity is chosen by the implementation:
    paren = (x, y) -> "($x,$y)"
    # Next, we construct data `xs` such that `length(xs)` is greater
    # than short array cutoff of `_mapreduce`:
    alphabets = 'a':'z'
    blksize = Base.pairwise_blocksize(identity, paren) ÷ length(alphabets)
    xs = repeat(alphabets, 2 * blksize)
    @test length(xs) > blksize
    # So far we constructed the data `xs` and reducing function
    # `paren` such that `reduce` and `foldl` results are different.
    # That is to say, this `reduce` does not hit the fall-back `foldl`
    # branch:
    @test foldl(paren, xs) != reduce(paren, xs)

    # Now let's try it with `Broadcasted`:
    bcraw = Broadcast.broadcasted(identity, xs)
    bc = Broadcast.instantiate(bcraw)
    # If `Broadcasted` has `IndexLinear` style, it should hit the
    # `reduce` branch:
    @test IndexStyle(bc) == IndexLinear()
    @test reduce(paren, bc) == reduce(paren, xs)
    # If `Broadcasted` does not have `IndexLinear` style, it should
    # hit the `foldl` branch:
    @test IndexStyle(bcraw) == IndexCartesian()
    @test reduce(paren, bcraw) == foldl(paren, xs)

    # issue #41055
    bc = Broadcast.instantiate(Broadcast.broadcasted(Base.literal_pow, Ref(^), [1,2], Ref(Val(2))))
    @test sum(bc, dims=1, init=0) == [5]
    bc = Broadcast.instantiate(Broadcast.broadcasted(*, ['a','b'], 'c'))
    @test prod(bc, dims=1, init="") == ["acbc"]
end

# treat Pair as scalar:
@test replace.(split("The quick brown fox jumps over the lazy dog"), r"[aeiou]"i => "_") ==
      ["Th_", "q__ck", "br_wn", "f_x", "j_mps", "_v_r", "th_", "l_zy", "d_g"]

# 28680
@test 1 .+ 1 .+  (1, 2) == (3, 4)

# PR #35260 no allocations in simple broadcasts
u = rand(100)
k1 = similar(u)
k2 = similar(u)
k3 = similar(u)
k4 = similar(u)
f(a,b,c,d,e) = @. a = a + 1*(b+c+d+e)
@allocated f(u,k1,k2,k3,k4)
@test (@allocated f(u,k1,k2,k3,k4)) == 0

ret =  @macroexpand @.([Int, Number] <: Real)
@test ret == :([Int, Number] .<: Real)

ret =  @macroexpand @.([Int, Number] >: Real)
@test ret == :([Int, Number] .>: Real)

# Threw mapany not defined
p = rand(4,4); r = rand(2,4);
p0 = copy(p)
@views @. p[1:2, :] += r
@test p[1:2, :] ≈ p0[1:2, :] + r

@test identity(.+) == Broadcast.BroadcastFunction(+)
@test identity.(.*) == Broadcast.BroadcastFunction(*)
@test map(.+, [[1,2], [3,4]], [5, 6]) == [[6,7], [9,10]]
@test repr(.!) == "Base.Broadcast.BroadcastFunction(!)"
@test eval(:(.+)) == Base.BroadcastFunction(+)

@testset "Issue #5187: Broadcasting of short-circuiting ops" begin
    ex = Meta.parse("A .< 1 .|| A .> 2")
    @test ex == :((A .< 1) .|| (A .> 2))
    @test ex.head == :.||
    ex = Meta.parse("A .< 1 .&& A .> 2")
    @test ex == :((A .< 1) .&& (A .> 2))
    @test ex.head == :.&&

    A = -1:4
    @test (A .< 1 .|| A .> 2) == [true, true, false, false, true, true]
    @test (A .>= 1 .&& A .<= 2) == [false, false, true, true, false, false]

    mutable struct F5187; x; end
    (f::F5187)(x) = (f.x += x)
    @test (iseven.(1:4) .&& (F5187(0)).(ones(4))) == [false, 1, false, 2]
    @test (iseven.(1:4) .|| (F5187(0)).(ones(4))) == [1, true, 2, true]
    r = 1:4; o = ones(4); f = F5187(0);
    @test (@. iseven(r) && f(o)) == [false, 1, false, 2]
    @test (@. iseven(r) || f(o)) == [3, true, 4, true]

    @test (iseven.(1:8) .&& iseven.((F5187(0)).(ones(8))) .&& (F5187(0)).(ones(8))) == [false,false,false,1,false,false,false,2]
    @test (iseven.(1:8) .|| iseven.((F5187(0)).(ones(8))) .|| (F5187(0)).(ones(8))) == [1,true,true,true,2,true,true,true]
    r = 1:8; o = ones(8); f1 = F5187(0); f2 = F5187(0)
    @test (@. iseven(r) && iseven(f1(o)) && f2(o)) == [false,false,false,1,false,false,false,2]
    @test (@. iseven(r) || iseven(f1(o)) || f2(o)) == [3,true,true,true,4,true,true,true]
    @test (iseven.(1:8) .&& iseven.((F5187(0)).(ones(8))) .&& (F5187(0)).(ones(8))) == [false,false,false,1,false,false,false,2]
    @test (iseven.(1:8) .|| iseven.((F5187(0)).(ones(8))) .|| (F5187(0)).(ones(8))) == [1,true,true,true,2,true,true,true]
end

@testset "Issue #28382: inferrability of broadcast with Union eltype" begin
    @test isequal([1, 2] .+ [3.0, missing], [4.0, missing])
    @test Core.Compiler.return_type(broadcast, Tuple{typeof(+), Vector{Int},
                                                     Vector{Union{Float64, Missing}}}) ==
        Union{Vector{Missing}, Vector{Union{Missing, Float64}}, Vector{Float64}}
    @test Core.Compiler.return_type(+, Tuple{Vector{Int},
                                             Vector{Union{Float64, Missing}}}) ==
        Union{Vector{Missing}, Vector{Union{Missing, Float64}}, Vector{Float64}}
    @test isequal(tuple.([1, 2], [3.0, missing]), [(1, 3.0), (2, missing)])
    @test Core.Compiler.return_type(broadcast, Tuple{typeof(tuple), Vector{Int},
                                                     Vector{Union{Float64, Missing}}}) ==
        Union{Vector{Tuple{Int, Missing}}, Vector{Tuple{Int, Any}}, Vector{Tuple{Int, Float64}}}
    # Check that corner cases do not throw an error
    @test isequal(broadcast(x -> x === 1 ? nothing : x, [1, 2, missing]),
                  [nothing, 2, missing])
    @test isequal(broadcast(x -> x === 1 ? nothing : x, Any[1, 2, 3.0, missing]),
                  [nothing, 2, 3, missing])
    @test broadcast((x,y)->(x==1 ? 1.0 : x, y), [1 2 3], ["a", "b", "c"]) ==
        [(1.0, "a") (2, "a") (3, "a")
         (1.0, "b") (2, "b") (3, "b")
         (1.0, "c") (2, "c") (3, "c")]
    @test typeof.([iszero, isdigit]) == [typeof(iszero), typeof(isdigit)]
    @test typeof.([iszero, iszero]) == [typeof(iszero), typeof(iszero)]
    @test isequal(identity.(Vector{<:Union{Int, Missing}}[[1, 2],[missing, 1]]),
                  [[1, 2],[missing, 1]])
end

@testset "Issue #28382: eltype inconsistent with getindex" begin
    struct Cyclotomic <: Number
    end

    Base.eltype(::Type{<:Cyclotomic}) = Tuple{Int,Int}

    Base.:*(c::T, x::Cyclotomic) where {T<:Real} = [1, 2]
    Base.:*(x::Cyclotomic, c::T) where {T<:Real} = [1, 2]

    @test Cyclotomic() .* [2, 3] == [[1, 2], [1, 2]]
end

@testset "inplace broadcast with trailing singleton dims" begin
    for (a, b, c) in (([1, 2], reshape([3 4], :, 1), reshape([5, 6], :, 1, 1)),
            ([1 2; 3 4], reshape([5 6; 7 8], 2, 2, 1), reshape([9 10; 11 12], 2, 2, 1, 1)))

        a_ = copy(a)
        a_ .= b
        @test a_ == dropdims(b, dims=(findall(==(1), size(b))...,))

        a_ = copy(a)
        a_ .= b
        @test a_ == dropdims(b, dims=(findall(==(1), size(b))...,))

        a_ = copy(a)
        a_ .= b .+ c
        @test a_ == dropdims(b .+ c, dims=(findall(==(1), size(c))...,))

        a_ = copy(a)
        a_ .*= c
        @test a_ == dropdims(a .* c, dims=(findall(==(1), size(c))...,))
    end
end

@testset "Issue #40309: still gives a range after #40320" begin
    @test Base.broadcasted_kwsyntax(+, [1], [2]) isa Broadcast.Broadcasted{<:Any, <:Any, typeof(+)}
    @test Broadcast.BroadcastFunction(+)(2:3, 2:3) == 4:2:6
    @test Broadcast.BroadcastFunction(+)(2:3, 2:3) isa AbstractRange
end

@testset "#42063" begin
    buf = IOBuffer()
    @test println.(buf, [1,2,3]) == [nothing, nothing, nothing]
    @test String(take!(buf)) == "1\n2\n3\n"
end

@testset "Memory allocation inconsistency in broadcasting #41565" begin
    function test(y)
        y .= 0 .- y ./ (y.^2) # extra allocation
        return y
    end
    arr = rand(1000)
    @allocated test(arr)
    @test (@allocated test(arr)) == 0
end

max_stars_repo_name: pepijndevos/julia

max_issues_repo_name: pepijndevos/julia

max_forks_repo_name: pepijndevos/julia


text:
using LinearAlgebra, DataFrames, Findpeaks, Smoothers, Interpolations, LsqFit

"""
    expected_energy(θ, Eγ)

Compute expected compton energy from scattering off electrons.

# Arguments
- `θ`: scattering angle in degrees
- `Eγ`: photon source energy
"""
function expected_energy(θ, Eγ)
    α = Eγ / (511)
    E = Eγ * (1 / (1 + α*(1 - cos(θ * (π / 180)))))
    return E
end

"""
    read_Txt(path_to_file::String)

Create DataFrame from .Txt file generated by Maestro.

To get .Txt file from Maestro: `file > print > print to text`
and save as a .Txt file in a location you can easily locate.
"""
function read_Txt(path_to_file::String)

    file = open(path_to_file, "r")

    ## collect data. Note that this is 1-indexed while data is 0-indexedx
    data = Vector{Int64}()

    for ln in eachline(file)
        if !isempty(ln) && ln[1] == ' '
            tempNum::String = "";
            dataPass = false
            for i in 1:length(ln)
                if dataPass && ln[i] != ' '
                    tempNum *= ln[i]
                elseif dataPass && ln[i - 1] != ' ' && ln[i] == ' ' && !isempty(tempNum)
                    append!(data, parse(Int64, tempNum)) # add number to vector
                    tempNum = "";
                end
                if dataPass && i == length(ln) # if add the end of a row, append number to vector
                    append!(data, parse(Int64, tempNum)) # add number to vector
                    tempNum = "";
                end
                if ln[i] == ':' dataPass = true end
            end
        end
    end
    close(file)

    # add indices of each channel starting with 1
    lines = [i for i in 1:length(data)]

    return DataFrame(Lines = lines, Counts = data)
end

"""
    get_xy_data(data::DataFrame, lowerLim::Int64, f::Function)

Scale x-axis to computed energy calibration via the linear function `f`

# Arguments
- `data::DataFrame`: data to be scaled
- `lowerLim::Int64`: cut off data below this index
- `f::Function`: computed calibration
"""
function get_xy_data(data::DataFrame, lowerLim::Int64, f::Function)
    xCh = data[!, 1]
    y = data[!, 2]

    #exclude low energy noise
    xCh = xCh[lowerLim:length(xCh)]
    y = y[lowerLim:length(y)]

    return f.(xCh), y
end

"""
    select_part(x0::Vector, y0::Vecotor, lower, upper)

Selects data between `lower` and `upper` limits of the x-axis
"""
function select_part(x0, y0, lower, upper)

    lims = Vector{Int64}()

    #find limits for filtering ydata
    for i in 1:length(y0)
        if x0[i] > lower
            push!(lims, i)
            break
        end
    end
    for i in 1:length(y0)
        if x0[length(x0) - i] < upper
            push!(lims, length(x0) - i)
            break
        end
    end

    return x0[lims[1]:lims[2]], y0[lims[1]:lims[2]]
end

"""
    weighted_mean(x::Vector, w::Vector)

Compute weighted mean of `x` with weights `w`
"""
function weighted_mean(x, w)
    #find weighted mean
    meanVal = 0.0
    w = normalize(w, 1)
    for i in 1:length(x)
        meanVal += x[i] * w[i]
    end
    #find uncertainty in mean
    uncertaintyVal = 0.0
    for i in 1:length(x)
        uncertaintyVal += w[i] * (x[i] - meanVal)^2 * (1/(1 - w[i]^2))
    end
    return meanVal, sqrt(uncertaintyVal);
end

"""
    remove_missing(x::Vector)

Remove `missing` values from a vector
"""
function remove_missing(x::Vector)
    xNew = Vector{Float64}()
    for i in 1:length(x)
        if ismissing(x[i])
            push!(xNew, 0.0)
        else
            push!(xNew, x[i])
        end
    end
    return xNew
end

"""
    find_peak_ends(x::Vector, y::Vector, peak::Float64, tol::Float64)

Compute endings of peak up to a give tolerance `tol`.
The value of `tol` represents the value of the derivative at which to
truncate the peak.

See also [`peak_parameters`](@ref), [`get_peak_value`](@ref)
"""
function find_peak_ends(x::Vector, y::Vector, peak, tol)

    #create interpolation
    itp = interpolate((x,), y, Gridded(Linear()));
    dx(x) = only(Interpolations.gradient(itp, x))

    # vector of outputs
    turningPoints = Vector{Int64}()

    #going left from the peak
    for i in 1:(peak + 20)
        if abs(dx(x[peak - i - 20])) < tol
            push!(turningPoints, peak - i - 20)
            break
        end
    end
    push!(turningPoints, peak) #add peak
    #going to the right from the peak
    for i in 1:(length(x) - peak - 15)
        if abs(dx(x[peak + i + 20])) < tol
            push!(turningPoints, peak + i + 20)
            break
        end
    end

    return turningPoints
end

"""
    peak_parameters(x, y0, SmoothParam, promParam, tol)

Compute location of peaks and endings of peaks as well as
a smoothing of the data.

# Arguments
- `x::Vector`: x-values of graph
- `y0::Vector`: y-values of graph
- `SmoothParam::Int`: size of weighted average region for data smoothing
- `promParam::Float64`: prominance of peaks to search for
- `tol::Float64`: tolance of location of peak ending
```
See also [`find_peak_ends`](@ref), [`get_peak_value`](@ref)
"""
function peak_parameters(x::Vector, y0::Vector, SmoothParam::Int, promParam::Float64, tol::Float64)

    #smooth data
    y = remove_missing(sma(y0, SmoothParam, true))

    #find peaks
    peaks = findpeaks(y, x, min_prom=promParam)

    #get peaks and ends of peaks
    peakData = Vector{Vector{Int64}}()
    for i in 1:length(peaks)
        append!(peakData, [find_peak_ends(x, y, peaks[i], tol)])
    end

    return peakData, y
end

"""
    get_peak_value(x::Vector, y::Vector, peak::Vector)

Compute weighted mean about a peak

# Arguments
- `x::Vector`: x-values
- `y::Vector`: y-values
- `peak::Vector`: one of the peaks obtained from the function `peak_parameters()`

See also [`peak_parameters`](@ref), [`find_peak_ends`](@ref), [`weighted_mean`](@ref)
"""
function get_peak_value(x::Vector, y::Vector, peak::Vector)

    #find value of the background
    background = (y[peak[1]] + y[peak[3]]) / 2
    filterPeak = 0.2 * (y[peak[2]] - background) + background

    #find indices of where 20% of background intersects with plot
    scaledVals = Vector{Int64}()
    #from left edge to right
    for i in 1:length(y)
        if y[peak[1] + i] > filterPeak
            push!(scaledVals, peak[1] + i)
            break
        end
    end
    #from right edge to left
    for i in 1:length(y)
        if y[peak[3] - i] > filterPeak
            push!(scaledVals, peak[3] - i)
            break
        end
    end

    #compute weighted mean over this new section
    xScaled, yScaled = x[scaledVals[1]:scaledVals[2]], y[scaledVals[1]:scaledVals[2]]
    peakVal = weighted_mean(xScaled, yScaled)

    return peakVal;

end

"""
    gauss(x::Vector, p::Vector)

Compute standard gaussion from `x` data with fit parameters `p`

# Fit parameters
- `p[1]`: scaling constant
- `p[2]`: σ
- `p[3]`: μ
"""
@. gauss(x, p) = p[1] * (1/(p[2] * sqrt(2*π))) * exp((-1/2)*((x - p[3])^2) / p[2]^2)

"""
    fit_to_gauss(x::Vector, y::Vector, pk::Vector)

Fit data to gaussian around a peak
"""
function fit_to_gauss(x::Vector, y::Vector, pk::Vector)

    #find value of the background
    background = (y[pk[1]] + y[pk[3]]) / 2
    filterPeak = 0.2 * (y[pk[2]] - background) + background

    #find indices of where 20% of background intersects with plot
    scaledVals = Vector{Int64}()
    #from left edge to right
    for i in 1:length(y)
        if y[pk[1] + i] > filterPeak
            push!(scaledVals, pk[1] + i)
            break
        end
    end
    #from right edge to left
    for i in 1:length(y)
        if y[pk[3] - i] > filterPeak
            push!(scaledVals, pk[3] - i)
            break
        end
    end

    # starting conditions
    p0 = [500.0, 500.0, 500.0]
    fit = curve_fit(gauss, x[scaledVals[1]:scaledVals[2]], y[scaledVals[1]:scaledVals[2]], p0) #perform fit
    Params = fit.param
    Errors = standard_errors(fit)

    return Params, Errors

end

ext: jl


text:
using BenchmarkTools, ConvolutionalOperatorLearning
using Random

const SUITE = BenchmarkGroup()

## CAOL 2D benchmarks: 62 images (128 x 128)
SUITE["CAOL 2D"] = BenchmarkGroup()

rng = MersenneTwister(1)
x = [randn(rng,128,128) for l in 1:62]

R = (3,3)
H0 = generatefilters(:DCT,R,form=:matrix)
λ = 0.0001

maxiters = 30

SUITE["CAOL 2D"]["notrace,30"]  = @benchmarkable CAOL($x,$λ,($H0,$R),maxiters=30)
SUITE["CAOL 2D"]["trace,30"]    = @benchmarkable CAOL($x,$λ,($H0,$R),maxiters=30,trace=true)

SUITE["CAOL 2D"]["notrace,200"] = @benchmarkable CAOL($x,$λ,($H0,$R),maxiters=200)
SUITE["CAOL 2D"]["trace,200"]   = @benchmarkable CAOL($x,$λ,($H0,$R),maxiters=200,trace=true)

for p in 2:2:size(H0,2)-1
    SUITE["CAOL 2D"]["p=$p"] = @benchmarkable CAOL($x,$λ,($H0,$R),p=$p,maxiters=30)
end

# todo: add more benchmarks, incorporate in package ci appropriately?


text:
module Account

import HTTP, JSON3, Dates

export listAccounts, getAccount, getAccountSummary, getAccountInstruments, setAccountConfig

# TODO: Remove this?
include("Position.jl")
include("Trade.jl")
include("Order.jl")

# TODO: Add strict types

"""
The account struct given by Oanda

# Fields
- 'NAV': The Net Asset Value of an account
- 'alias': User defined alias if one exists
- 'balance': Current Account Balance
- 'createdByUserID': The User ID of the account creator
- 'createdTime': The time the account was created
- 'currency': The primary currency of the account
- 'hedgingEnabled': If the account is allowed to hedge
- 'id': The account ID
- 'lastTransactionID': The last transaction ID
- 'marginAvailable': The margin still available on the account
- 'marginCloseoutMarginUsed': The closeout margin used
- 'marginCloseoutNAV': Margins closeout NAV
- 'marginCloseoutPercent': Margin closeout percent
- 'marginCloseoutPositionValue': Margin closeout position value
- 'marginCloseoutUnrealizedPL': Margin closeout unrealised profit/loss
- 'marginRate': The margin rate
- 'marginUsed': Amount of margin used
- 'openPositionCount': Number of open positions
- 'openTradeCount': Number of open trades
- 'orders': Orders of the account
- 'pendingOrderCount': Number of pending orders
- 'pl': The profit or loss over the lifetime of the account
- 'positionValue': Value of an accounts open positions
- 'positions': Positions of the account
- 'resettablePL': The resetable profit/loss since last reset
- 'trades': Trades of the account
- 'unrealizedPL': The unrealised profit/loss of the account
- 'withdrawalLimit': The withdrawal limit of the account
"""
mutable struct account
    NAV # The Net Asset Value of an account
    alias # User defined alias if one exists
    balance # Current Account Balance
    createdByUserID # The User ID of the account creator
    createdTime # The time the account was created
    currency # The primary currency of the account
    hedgingEnabled # If the account is allowed to hedge
    id # The account ID
    lastTransactionID # The last transaction ID
    marginAvailable # The margin still available on the account
    marginCloseoutMarginUsed # The closeout margin used
    marginCloseoutNAV # Margins closeout NAV
    marginCloseoutPercent # Margin closeout percent
    marginCloseoutPositionValue # Margin closeout position value
    marginCloseoutUnrealizedPL # Margin closeout unrealised profit/loss
    marginRate # The margin rate
    marginUsed # Amount of margin used
    openPositionCount # Number of open positions
    openTradeCount # Number of open trades
    orders::Vector{Order.orderRequest} # Orders of the account
    pendingOrderCount # Number of pending orders
    pl # The profit or loss over the lifetime of the account
    positionValue # Value of an accounts open positions
    positions::Vector{Position.position} # Positions of the account
    resettablePL # The resetable profit/loss since last reset
    trades::Vector{Trade.trade} # Trades of the account
    unrealizedPL # The unrealised profit/loss of the account
    withdrawalLimit # The withdrawal limit of the account

    account() = new()
end

"Nessecary for automatic JSON parsing, not for regular use"
mutable struct topLayer
    account::account

    topLayer() = new()
end

"The ID and tag of each account"
mutable struct accountListed
    id # Account id
    tags # Account Tags

    accountListed() = new()
end

"The list of accounts returned by Oanda"
mutable struct accountsList
    accounts::Vector{accountListed} # Array of accounts

    accountsList() = new()
end

"""
Tradeable Instrument data

# Fields
- 'displayName': Instrument name
- 'displayPrecision': Decimal precision of the instrument
- 'marginRate': Margin rate on the instrument
- 'maximumOrderUnits': Max units that can be ordered
- 'maximumPositionSize': max position size of the instrument
- 'maximumTrailingStopDistance': max trailing stop distance
- 'minimumTrailingStopDistance': min trailing stop distance
- 'name': Request usable instrument name
- 'pipLocation': current pip location
- 'tradeUnitsPrecision': Decimal precision of trade units
- 'type': Type of instrument
"""
mutable struct instrumentDetail
    displayName # Instrument name
    displayPrecision # Decimal precision of the instrument
    marginRate # Margin rate on the instrument
    maximumOrderUnits # Max units that can be ordered
    maximumPositionSize # max position size of the instrument
    maximumTrailingStopDistance # max trailing stop distance
    #minimumPositionSize # min position size of the instrument
    minimumTrailingStopDistance # min trailing stop distance
    name # Request usable instrument name
    pipLocation # current pip location
    tradeUnitsPrecision # Decimal precision of trade units
    type # Type of instrument

    instrumentDetail() = new()
end

"Nessecary for automatic JSON parsing, not for regular use"
mutable struct instrumentTopLayer
    instruments::Vector{instrumentDetail} # Tradable instruments

    instrumentTopLayer() = new()
end

"For configuring Accounts"
struct accountConfig
    alias::String
    marginRate::String
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{account}) = JSON3.Mutable()
JSON3.StructType(::Type{topLayer}) = JSON3.Mutable()
JSON3.StructType(::Type{accountsList}) = JSON3.Mutable()
JSON3.StructType(::Type{accountListed}) = JSON3.Mutable()
JSON3.StructType(::Type{instrumentDetail}) = JSON3.Mutable()
JSON3.StructType(::Type{instrumentTopLayer}) = JSON3.Mutable()
JSON3.StructType(::Type{accountConfig}) = JSON3.Struct()

"Coerce a given Account Summary into its proper types (Used for getAccount Method)"
function coerceAccountSummary(acc::account)
    acc.NAV = parse(Float32, acc.NAV)
    acc.balance = parse(Float32, acc.balance)
    acc.marginAvailable = parse(Float32, acc.marginAvailable)
    acc.marginCloseoutMarginUsed = parse(Float32, acc.marginCloseoutMarginUsed)
    acc.marginCloseoutNAV = parse(Float32, acc.marginCloseoutNAV)
    acc.marginCloseoutPercent = parse(Float32, acc.marginCloseoutPercent)
    acc.marginCloseoutPositionValue = parse(
        Float32,
        acc.marginCloseoutPositionValue,
    )
    acc.marginCloseoutUnrealizedPL = parse(
        Float32,
        acc.marginCloseoutUnrealizedPL,
    )
    acc.marginRate = parse(Float32, acc.marginRate)
    acc.marginUsed = parse(Float32, acc.marginUsed)
    #acc.openPositionCount = parse(Int32, acc.openPositionCount)
    #acc.openTradeCount = parse(Int32, acc.openTradeCount)
    #acc.pendingOrderCount = parse(Int32, acc.pendingOrderCount)
    acc.pl = parse(Float32, acc.pl)
    acc.positionValue = parse(Float32, acc.positionValue)
    acc.resettablePL = parse(Float32, acc.resettablePL)
    acc.unrealizedPL = parse(Float32, acc.unrealizedPL)
    acc.withdrawalLimit = parse(Float32, acc.withdrawalLimit)

    return acc
end

"Coerce a given Account into its proper types (Used for getAccount Method)"
function coerceAccount(acc::account)
    temp = Vector{Order.orderRequest}()
    for order in acc.orders
        order = Order.coerceOrder(order)
        push!(temp, order)
    end
    acc.orders = temp

    temp = Vector{Trade.trade}()
    for trade in acc.trades
        trade = Trade.coerceTrade(trade)
        push!(temp, trade)
    end
    acc.trades = temp

    temp = Vector{Position.position}()
    for position in acc.positions
        position = Position.coercePos(position)
        push!(temp, position)
    end
    acc.positions = temp

    acc = coerceAccountSummary(acc)

    return acc
end

"Coerce a given instrument detail into its proper types"
function coerceInstrumentDetail(inst::instrumentDetail)
    inst.marginRate = parse(Float32, inst.marginRate)
    inst.maximumOrderUnits = parse(Float64, inst.maximumOrderUnits)
    inst.maximumPositionSize = parse(Float32, inst.maximumPositionSize)
    inst.maximumTrailingStopDistance = parse(
        Float32,
        inst.maximumTrailingStopDistance,
    )
    inst.minimumTrailingStopDistance = parse(
        Float32,
        inst.minimumTrailingStopDistance,
    )

    return inst
end

"""
    listAccounts(config::config)

Returns a list of all account IDs and tags authorized for the given Token

Returns object of type 'Vector{accountListed}'

# Arguments
- 'config::config': a valid config object

"""
function listAccounts(config)
    r = HTTP.request(
        "GET",
        string("https://", config.hostname, "/v3/accounts"),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )

    data = JSON3.read(r.body, accountsList)

    return data.accounts

end

"""
    getAccount(config::config)

Returns an Oanda account struct when given a valid config

Returns object of type 'account'

# Arguments
- 'config::config': a valid config object

"""
function getAccount(config)
    r = HTTP.request(
        "GET",
        string("https://", config.hostname, "/v3/accounts/", config.account),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )

    temp = JSON3.read(r.body, topLayer)

    temp = temp.account

    # Type Coersions
    temp = coerceAccount(temp)

    return temp

end

"""
    getAccountSummary(config::config)

Similar to getAccount but doesnt return the order & trade & positions lists, however
it still returns a full account struct, just with these fields left undefined

Returns object of type 'account'

# Arguments
- 'config::config': a valid config object

"""
function getAccountSummary(config)
    r = HTTP.request(
        "GET",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/summary",
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )

    temp = JSON3.read(r.body, topLayer)

    temp = temp.account

    # Type Coersions
    temp = coerceAccountSummary(temp)

    return temp

end

"""
    getAccountInstruments(config::config, inst=nothing)

Returns a list of tradeable instruments details for the account

Returns object of type 'Vector{instrumentDetail}'

# Arguments
- 'config::config': a valid config object
- 'inst::String': Can be left blank to return all tradeable instruments, or as a string csv of instruments to return their details

"""
function getAccountInstruments(config, inst::String = "")
    request = string(
        "https://",
        config.hostname,
        "/v3/accounts/",
        config.account,
        "/instruments",
    )
    if inst != ""
        request = string(request, "?instruments=", inst)
    end
    r = HTTP.request(
        "GET",
        request,
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
        ],
    )

    data = JSON3.read(r.body, instrumentTopLayer)
    data = data.instruments
    instruments = Vector{instrumentDetail}()
    for inst in data
        inst = coerceInstrumentDetail(inst)
        push!(instruments, inst)
    end

    return instruments
end

"""
    setAccountConfig(config::config, alias::String, marginRate::String)

Set client configurable configuration settings

Returns true on success

# Arguments
= 'config::config': a valid config object
- 'alias::String': The account alias
- 'marginRate::String': The desired decimal margin rate formatted as a string

"""
function setAccountConfig(config, alias::String, marginRate::String)
    data = accountConfig(alias, marginRate)
    r = HTTP.request(
        "PATCH",
        string(
            "https://",
            config.hostname,
            "/v3/accounts/",
            config.account,
            "/configuration",
        ),
        [
         "Authorization" => string("Bearer ", config.token),
         "Accept-Datetime-Format" => config.datetime,
         "Content-Type" => "application/json",
        ],
        JSON3.write(data),
    )

    return true
end

end

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
cd(@__DIR__)
length(workers()) == 1 && @show addprocs(4)

# @everywhere using Revise
using ParallelDataTransfer
@everywhere include(Pkg.dir("DynamicMovementPrimitives","src","two_link.jl"))
@everywhere using TwoLink, Parameters, JacProp, OrdinaryDiffEq, LTVModels, LTVModelsBase
@everywhere using Flux: params, jacobian
@everywhere using Flux.Optimise: Param, optimiser, expdecay
@everywhere begin
    BLAS.set_num_threads(1)
    @with_kw struct TwoLinkSys
        N  = 1000
        nu  = 2
        nx = 4
        h = 0.02
        σ0 = 0
        sind = 1:nx
        uind = nx+1:(nx+nu)
        s1ind = (nx+nu+1):(nx+nu+nx)
    end

    function generate_data(sys::TwoLinkSys, seed, validation=false; ufun=u->filt(ones(200),[200], 20u')')
        @unpack N, nu, h, σ0 = sys
        srand(seed)
        done = false
        local x,u
        while !done
            fs = logspace(-3,-1.5,50)[randperm(50)[1:5]]
            # u    = ufun(randn(nu,N+2))
            u = sum(f->sin.(2π*f.* (1:N+2)), fs)'
            u = [u;reverse(u')']
            t    = 0:h:N*h
            x0   = [-0.4,0,0,0]
            prob = OrdinaryDiffEq.ODEProblem((x,p,t)->time_derivative(x, u[:,floor(Int,t/h)+1]),x0,(t[[1,end]]...))
            sol  = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
            x    = hcat(sol(t)...)
            done = all(abs.(x[1:2,:]) .< 0.9π)
        end
        u = u[:,1:N]
        validation || (x .+= σ0 * randn(size(x)))
        @assert all(isfinite, x)
        @assert all(isfinite, u)
        x,u
    end

    function true_jacobian(sys::TwoLinkSys, x::AbstractVector, u::AbstractVector)
        Jtrue = ForwardDiff.jacobian(x->time_derivative(x, u), x) # TODO: why the indexing
        Jtrue = [Jtrue ForwardDiff.jacobian(u->time_derivative(x, u), u)]
        Jtrue = expm([sys.h*Jtrue;zeros(2,6)])[1:4,:] # Discretize
    end

    function callbacker(epoch, loss,d,trace,model,mt)
        i = length(trace) + epoch - 1
        function ()
            l = sum(d->Flux.data(loss(d...)),d)
            increment!(trace,epoch,l)
            i % 500 == 0 && println(@sprintf("Loss: %.4f", l))
        end
    end

    num_params = 40
    wdecay     = 0
    stepsize   = 0.01
    const sys  = TwoLinkSys(N=200, h=0.02, σ0 = 0.0)
    true_jacobian(x,u) = true_jacobian(sys,x,u)
    nu         = sys.nu
    nx         = sys.nx

end

## Generate validation data
function valdata()
    vx,vu,vy = Vector{Float64}[],Vector{Float64}[],Vector{Float64}[]
    for i = 20:60
        x,u = generate_data(sys,i, true)
        for j in 10:5:(sys.N-1)
            push!(vx, x[:,j])
            push!(vy, x[:,j+1])
            push!(vu, u[:,j])
        end
    end
    hcat(vx...),hcat(vu...),hcat(vy...)
end
vx,vu,vy = valdata()
vt = Trajectory(vx,vu,vy)

trajs = [Trajectory(generate_data(sys, i)...) for i = 1:3]

sendto(workers(), trajs=trajs, vt=vt)
#
# ## Without jacprop
# f2 = @spawnat 2 begin
#     srand(1)
#     models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
#     opts       = ADAM.(params.(models), stepsize, decay=0.0005)#; [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
#     trainer  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker, P = 10, R2 = 10000I, σdivider = 20)
#     for i = 1:3
#         trainer(trajs[i], epochs=20000, jacprop=0, useprior=false)
#         # traceplot(trainer)
#     end
#     trainer
# end
#
#
# ## With jacprop and prior
# f3 = @spawnat 3 begin
#     srand(1)
#     models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
#     opts       = ADAM.(params.(models), stepsize, decay=0.0005)#; [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
#     trainerj  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker, P = 10, R2 = 1000I, σdivider = 20)
#     for i = 1:3
#         trainerj(trajs[i], epochs=10000, jacprop=1, useprior=true)
#         # traceplot(trainerj)
#     end
#     trainerj
# end
#
# ## With jacprop no prior
# f4 = @spawnat 4 begin
#     srand(1)
#     models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
#     opts       = ADAM.(params.(models), stepsize, decay=0.0005)#; [expdecay(Param(p), wdecay) for p in params(models[i]) if p isa AbstractMatrix]] for i = 1:length(models)]
#     trainerjn  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker, P = 10, R2 = 1000I, σdivider = 200)
#     for i = 1:3
#         trainerjn(trajs[i], epochs=100, jacprop=1, useprior=false)
#         # traceplot(trainerjn)
#     end
#     trainerjn(epochs=10000)
#     trainerjn
# end
# # predsimplot(trainerjn.models, trajs[1], title=["\$q_1\$" "\$\\dotq_1\$" "\$q_2\$" "\$\\dotq_2\$"])
# # predsimplot(trainerjn.models, Trajectory(generate_data(sys, 6)...), title=["\$q_1\$" "\$\\dotq_1\$" "\$q_2\$" "\$\\dotq_2\$"])
#
# pyplot(reuse=false)
# trainer,trainerj,trainerjn = fetch(f2), fetch(f3),fetch(f4)
#
# mutregplot(trainer, vt, true_jacobian, title="Witout jacprop", subplot=1, layout=(2,3), reuse=false, useprior=false)
# mutregplot!(trainerj, vt, true_jacobian, title="With jacprop and prior", subplot=2, link=:y, useprior=true, show=false)
# mutregplot!(trainerjn, vt, true_jacobian, title="With jacprop, no prior", subplot=3, link=:y, useprior=false)
# traceplot!(trainer, subplot=4)
# traceplot!(trainerj, subplot=5, show=false)
# traceplot!(trainerjn, subplot=6)
# gui()
# ##
#
#
# jacplot(trainer.models, vt, true_jacobian, label="Without", c=:red, reuse=false, fillalpha=0.2)
# # jacplot!(KalmanModel(trainer, vt), vt, label="Without", c=:pink)
# jacplot!(trainerj.models, vt, label="With", c=:blue, fillalpha=0.2)
# # jacplot!(KalmanModel(trainerj, vt), vt, label="With", c=:cyan)
# jacplot!(trainerjn.models, vt, label="With no", c=:green, fillalpha=0.2)
# # jacplot!(KalmanModel(trainerjn, vt), vt, label="With no", c=:green)
# gui()
#
#
# ui = display_modeltrainer(trainer, size=(800,600))
# jacplot(trainer.models, trainer.trajs[3], true_jacobian, ds=20)
# @gif for i = 1:length(t)
#     eigvalplot(trainer.models, trainer.trajs[1], true_jacobian, ds=20, onlyat=i)
# end
#
# # TODO: Validate vt in callback

## Monte-Carlo evaluation
num_montecarlo = 2

res = map(1:num_montecarlo) do it
    r2 = @spawn begin
        srand(it)
        model     = JacProp.ADDiffSystem(nx,nu,num_params,tanh) # TODO: tanh has no effect
        opt       = LTVModels.ADAMOptimizer(model.w, α = stepsize)
        trainerad = ADModelTrainer(;model=model, opt=opt, λ=0.01, testdata = vt)
        for i = 1:3
            trainerad(trajs[i], epochs=0)
        end
        trainerad(epochs=2000)
        trainerad
    end
    r1 = @spawn begin
        srand(it)
        # models     = [DiffSystem(nx,nu,num_params, a) for a in default_activations]
        models     = [DiffSystem(nx,nu,num_params, tanh)]
        opts       = ADAM.(params.(models), stepsize, decay=0.0005)#;
        trainer  = ModelTrainer(models = models, opts = opts, losses = JacProp.loss.(models), cb=callbacker, P = 10, R2 = 10000I, σdivider = 20)
        for i = 1:3
            trainer(trajs[i], epochs=0, jacprop=0, useprior=false)
            # traceplot(trainer)
        end
        trainer(epochs=2000)
        trainer
    end

    println("Done with montecarlo run $it")
    r1,r2
end
res = [(fetch.(rs)...) for rs in res]

# serialize("results", res)
# res = deserialize("results")
eigvalplot(res[1][1].models, vt, true_jacobian;layout=2,subplot=1,cont=false,title="Standard", ylims=[-0.1,0.1], xlims=[0.5,1.1])
eigvalplot!(res[1][2].model, vt, true_jacobian;subplot=2,cont=false,title="AD Jacprop", ylims=[-0.1,0.1], xlims=[0.5,1.1]);gui()
plot(res[1][2].trace.iterations,[res[2].trace.values for res in res], c=:blue)
plot!(res[1][2].trace.iterations,[res[2].tracev.values for res in res], c=:orange)
plot!(res[1][1].trace.iterations,[res[1].trace.values for res in res]./4, c=:red, xscale=:log10, yscale=:log10, legend=false)

resdiff = getindex.(res,1)
resad = getindex.(res,2)

nr = length(res[1])÷2
labelvec = ["f" "g"]
infostring = @sprintf("Num hidden: %d, sigma: %2.2f, Montecarlo: %d", num_params, sys.σ0, num_montecarlo)
simvals = [Trajectory(generate_data(sys, i)...) for i = 4:6]
pred  = [JacProp.eval_pred.(resdiff, vt) JacProp.eval_pred.(resad, vt)]
sim   = vcat([[JacProp.eval_sim.(resdiff, simval) JacProp.eval_sim.(resad, simval)] for simval in simvals]...)
jac   = [JacProp.eval_jac.(resdiff, vt, true_jacobian,3) JacProp.eval_jac.(resad, vt, true_jacobian,3)]

##
using StatPlots
vio1 = boxplot(pred, lab=["Standard" "Jacobian propagation"], ylabel="Prediction RMS", reuse=false, yscale=:log10)
vio2 = boxplot(min.((sim),2), lab=["Standard" "Jacobian propagation"], ylabel="Simulation RMS", yscale=:log10)
vio3 = boxplot(jac, lab=["Standard" "Jacobian propagation"], ylabel="Jacobian Error", yscale=:identity)
plot(vio1,vio2,vio3,title=infostring); gui()
##
# savefig2("/local/home/fredrikb/papers/nn_prior/figs/valerr.tex")

plot(trajs[1])
predsimplot(resdiff[1].models, trajs[3]); gui()
predsimplot(resad[1].model, trajs[3], reuse=false); gui()
# simulate(resad[1].model, trajs[1])' |> plot
##
i = 1
j = 3
plot(simvals[i], control=false, layout=4, lab="Traj")
simplot!(resdiff[j].models, simvals[i], lab="Standard")
simplot!(resad[j].model, simvals[i], lab="AD JAcprop")
##

jacplot(resdiff[1].models, trajs[1], true_jacobian, ds=2, reuse=false)
jacplot!(resad[1].model, trajs[1], ds=2); gui()


text:
# construct solver for fixed-point iterations
function build_fpsolver(alg, fpalg::Union{NLFunctional,NLAnderson}, u, uEltypeNoUnits,
                        uBottomEltypeNoUnits, ::Val{true})
  # no fixed-point iterations if the algorithm is constrained
  isconstrained(alg) && return

  # define unitless type
  uTolType = real(uBottomEltypeNoUnits)

  # could use integrator.cache.atmp if guaranteed to exist
  atmp = similar(u, uEltypeNoUnits)
  dz = similar(u)

  # build cache
  if fpalg isa NLFunctional
    fpcache = FPFunctionalCache(atmp, dz)
  elseif fpalg isa NLAnderson
    max_history = min(fpalg.max_history, fpalg.max_iter, length(u))
    Δz₊s = [zero(u) for i in 1:max_history]
    Q = Matrix{uEltypeNoUnits}(undef, length(u), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)

   dzold = zero(u)
   z₊old = zero(u)

   fpcache = FPAndersonCache(atmp,dz,dzold,z₊old,Δz₊s,Q,R,γs,0,fpalg.aa_start,fpalg.droptol)
  end

  # build solver
  ηold = one(uTolType)

  FPSolver{typeof(fpalg),true,uTolType,typeof(fpcache)}(
    fpalg, uTolType(fpalg.κ), uTolType(fpalg.fast_convergence_cutoff), ηold, 10000,
    fpalg.max_iter, SlowConvergence, fpcache, 0)
end

function build_fpsolver(alg, fpalg::Union{NLFunctional,NLAnderson}, u, uEltypeNoUnits,
                        uBottomEltypeNoUnits, ::Val{false})
  # no fixed-point iterations if the algorithm is constrained
  isconstrained(alg) && return

  # define unitless type
  uTolType = real(uBottomEltypeNoUnits)

  # build cache
  if fpalg isa NLFunctional
    fpcache = FPFunctionalConstantCache()
  elseif fpalg isa NLAnderson
    max_history = min(fpalg.max_history, fpalg.max_iter, length(u))
    Δz₊s = Vector{typeof(u)}(undef, max_history)
    Q = Matrix{uEltypeNoUnits}(undef, length(u), max_history)
    R = Matrix{uEltypeNoUnits}(undef, max_history, max_history)
    γs = Vector{uEltypeNoUnits}(undef, max_history)

    dz = u
    dzold = u
    z₊old = u

    fpcache = FPAndersonConstantCache(dz,dzold,z₊old,Δz₊s,Q,R,γs,0,fpalg.aa_start,
                                      fpalg.droptol)
  end

  # build solver
  ηold = one(uTolType)

  FPSolver{typeof(fpalg),false,uTolType,typeof(fpcache)}(
    fpalg, uTolType(fpalg.κ), uTolType(fpalg.fast_convergence_cutoff), ηold, 10_000,
    fpalg.max_iter, SlowConvergence, fpcache, 0)
end

## resize

function resize_fpsolver!(integrator::DDEIntegrator, i::Int)
  @unpack fpsolver = integrator

  if fpsolver !== nothing
    resize!(fpsolver, integrator, i)
  end

  nothing
end

function Base.resize!(fpsolver::FPSolver, integrator::DDEIntegrator, i::Int)
  resize!(fpsolver.cache, fpsolver, integrator, i)
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
using FiniteDiff, LinearAlgebra, SparseArrays, Test, StaticArrays

function f(x)
  xm1 = [0;x[1:end-1]]
  xp1 = [x[2:end];0]
  xm1 - 2x + xp1
end

function second_derivative_stencil(N)
  A = zeros(N,N)
  for i in 1:N, j in 1:N
      (j-i==-1 || j-i==1) && (A[i,j]=1)
      j-i==0 && (A[i,j]=-2)
  end
  A
end

x = @SVector ones(30)
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:forward}, eltype(x))
@test J ≈ second_derivative_stencil(30)

J = FiniteDiff.finite_difference_jacobian(f, x, Val{:central}, eltype(x))
@test J ≈ second_derivative_stencil(30)

J = FiniteDiff.finite_difference_jacobian(f, x, Val{:complex}, eltype(x))
@test J ≈ second_derivative_stencil(30)

spJ = sparse(second_derivative_stencil(30))
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:forward}, eltype(x), jac_prototype=spJ)
@test J ≈ second_derivative_stencil(30)
@test typeof(J) == typeof(spJ)
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:forward}, eltype(x),
    colorvec=SVector{30}(repeat(1:3, 10)), sparsity=spJ, jac_prototype=spJ)
@test J ≈ second_derivative_stencil(30)
@test typeof(J) == typeof(spJ)

f(x) = x
@testset "1x1 test of $x" for
  (x, y) in ((SVector{1}([1.]), SMatrix{1,1}), ([1.0], Matrix)),
    difftype in (:forward, :central, :complex)

  J = FiniteDiff.finite_difference_jacobian(f, x, Val{difftype}, eltype(x))
  @test J[1, 1] ≈ 1.0
  @test J isa y
end

x = SVector{1}([1.])
f(x) = vcat(x, x)
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:forward}, eltype(x))
@test J ≈ fill(1.0, 2, 1)
@test J isa SMatrix{2,1}
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:central}, eltype(x))
@test J ≈ fill(1.0, 2, 1)
@test J isa SMatrix{2,1}
J = FiniteDiff.finite_difference_jacobian(f, x, Val{:complex}, eltype(x))
@test J ≈ fill(1.0, 2, 1)
@test J isa SMatrix{2,1}