text:
using Jets, LinearAlgebra, Random, Test

JopFoo_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
function JopFoo(diag)
    spc = JetSpace(Float64, length(diag))
    JopLn(;df! = JopFoo_df!, dom = spc, rng = spc, s = (diagonal=diag,))
end
Jets.perfstat(J::Jet{D,R,typeof(JopFoo_df!)}) where {D,R} = Float64(π)

JopClose_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
function JopClose(diag)
    spc = JetSpace(Float64, length(diag))
    file = tempname()
    touch(file)
    JopLn(;df! = JopClose_df!, dom = spc, rng = spc, s = (diagonal=diag, file=file))
end
Base.close(J::Jet{D,R,typeof(JopClose_df!)}) where {D,R} = rm(state(J).file)

JopBar_f!(d,m;kwargs...) = d .= m.^2
JopBar_df!(δd,δm;mₒ,kwargs...) = δd .= 2 .* mₒ .* δm
function JopBar(n)
    spc = JetSpace(Float64, n)
    JopNl(f! = JopBar_f!, df! = JopBar_df!, dom = spc, rng = spc)
end

JopBaz_df!(d,m;A,kwargs...) =  d .= A * m
JopBaz_df′!(m,d;A,kwargs...) = m .= A' * d
function JopBaz(A)
    dom = JetSpace(eltype(A), size(A,2))
    rng = JetSpace(eltype(A), size(A,1))
    JopLn(;df! = JopBaz_df!, df′! = JopBaz_df′!, dom = dom, rng = rng, s = (A=A,))
end

JopFooBar_df!(d,m;A,kwargs...) =  d .= A.*m
function JopFooBar(n)
    spc = JetSpace(Float64,n,n)
    JopLn(df! = JopFooBar_df!, dom = spc, rng = spc, s= (A=rand(n,n),))
end

JopRosenbrock_f!(d,m;kwargs...) = d.= [1-m[1], 10*(m[2]-m[1]^2)]
JopRosenbrock_df!(d,m;J,kwargs...) = d .= J*m
JopRosenbrock_df′!(m,d;J,kwargs...) = m .= J'*d
JopRosenbrock_upstate!(m,s) = s.J[2,1] = -20.0*m[1]
function JopRosenbrock()
    dom = JetSpace(Float64, 2)
    rng = JetSpace(Float64, 2)
    JopNl(f! = JopRosenbrock_f!, df! = JopRosenbrock_df!, df′! = JopRosenbrock_df′!,
        upstate! = JopRosenbrock_upstate!, dom = dom, rng = rng, s = (J=[-1.0 0.0;0.0 10.0],))
end

JopFoo2_df!(d,m;diagonal,kwargs...) = d .= diagonal .* m
function JopFoo2(diag)
    spc = JetSpace(Float64, size(diag))
    JopLn(;df! = JopFoo2_df!, dom = spc, rng = spc, s = (diagonal=diag,))
end

@testset "JetSpace, construction, n=$n, T=$T" for n in ((2,),(2,3),(2,3,4)), T in (Float32,Float64,Complex{Float32},Complex{Float64})
    N = length(n)
    R = JetSpace(T, n...)
    @test size(R) == n
    @test eltype(R) == T
    @test eltype(typeof(R)) == T
    @test ndims(R) == N
end

@testset "JetSpace, operations, n=$n, T=$T" for n in ((2,),(2,3),(2,3,4)), T in (Float32,Float64,Complex{Float32},Complex{Float64})
    R = JetSpace(T, n...)
    N = length(n)
    x = rand(R)
    @test eltype(x) == T
    @test size(x) == n
    @test ndims(R) == N
    y = rand(T, n)
    z = copy(y)[:]
    @test y ≈ reshape(z,  R)
    @test ones(R) ≈ ones(T, n)
    @test size(rand(R)) == size(R)
    @test zeros(R) ≈ zeros(T, n)
    @test size(Array(R)) == size(R)
    @test length(randperm(R, 2)) == 2
    @test typeof(randperm(R, 2)) == Vector{Int}
    @test R == space(rand(R))
end

@testset "Jet, construction" begin
    f!(d,m;a) = d .= a .* m.^2
    df!(δd,δm;a,mₒ) = δd .= 2 .* a .* mₒ .* δm
    a = rand(20)
    ✈ = Jet(;
        dom = JetSpace(Float64,20),
        rng = JetSpace(Float64,10,2),
        f! = f!,
        df! = df!,
        df′! = df!,
        s = (a=a,))
    @test domain(✈) == JetSpace(Float64,20)
    @test range(✈) == JetSpace(Float64,10,2)
    @test point(✈) ≈ zeros(eltype(domain(✈)),ntuple(i->0,ndims(domain(✈))))
    mₒ = rand(domain(✈))
    Jets.point!(✈, mₒ)
    @test point(✈) ≈ mₒ
    a .= rand(20)
    state!(✈, (a=a,))
    s = state(✈)
    @test s.a ≈ a
    @test state(✈, :a) ≈ a
    @test shape(✈, 1) == (10,2)
    @test shape(✈, 2) == (20,)
    @test shape(✈) == ((10,2),(20,))
    @test size(✈, 1) == 20
    @test size(✈, 2) == 20
    @test close(✈) == false
end

@testset "linear operator" begin
    diag = rand(10)
    A = JopFoo(diag)
    m = rand(10)
    d = A*m
    @test d ≈ diag .* m
    a = A'*d
    @test a ≈ diag .* d

    d .= 0
    mul!(d, A, m)
    @test d ≈ diag .* m
    a .= 0
    mul!(a, A', d)
    @test a ≈ diag .* d

    @test size(A) == (10,10)
    @test shape(A) == ((10,), (10,))
    @test size(A,1) == 10
    @test size(A,2) == 10
    @test shape(A,1) == (10,)
    @test shape(A,2) == (10,)
    @test domain(A) == JetSpace(Float64,10)
    @test range(A) == JetSpace(Float64,10)
    @test eltype(A) == Float64
    @test convert(Array, A) ≈ diagm(0=>diag)

    @test state(A).diagonal ≈ diag
    @test state(A, :diagonal) ≈ diag

    B = JopFooBar(5)
    @test convert(Array, B) ≈ diagm(0=>vec(state(B).A))
end

@testset "nonlinear operator" begin
    n = 10
    F = JopBar(n)
    m = rand(domain(F))
    d = F*m
    @test d ≈ m.^2
    d .= 0
    mul!(d, F, m)
    @test d ≈ m.^2
    J = jacobian!(F, m)
    @test point(J) ≈ m
    d = J*m
    @test d ≈ 2 .* point(J) .* m
    a = J'*d
    @test a ≈ 2 .* point(J) .* d

    @test size(F) == (10,10)
    @test shape(F) == ((10,), (10,))
    @test size(F,1) == 10
    @test size(F,2) == 10
    @test shape(F,1) == (10,)
    @test shape(F,2) == (10,)
    @test domain(F) == JetSpace(Float64,10)
    @test range(F) == JetSpace(Float64,10)
end

@testset "upstate" begin
    F = JopRosenbrock()
    m = rand(2)
    J = jacobian!(F, m)
    @test state(J).J ≈ [-1.0 0.0;-20*m[1] 10.0]
end

@testset "multiple simultaneous linearizations" begin
    F = JopBar(2)
    J₁ = jacobian(F, [1.0,2.0])
    J₂ = jacobian(F, [3.0,4.0])

    δm = [1.0, 2.0]
    @test J₁ * δm ≈ 2 .* [1.0,2.0] .* δm
    @test J₂ * δm ≈ 2 .* [3.0,4.0] .* δm

    J₁ = jacobian!(F, [1.0,2.0])
    J₂ = jacobian!(F, [3.0,4.0])

    @test J₁ * δm ≈ 2 .* [3.0,4.0] .* δm
    @test J₂ * δm ≈ J₁ * δm
end

function indexmap(I)
    if I[1] < 5
        return CartesianIndex(I)
    else
        return CartesianIndex((I[1]-4,I[2]))
    end
end

@testset "Symmetric space" begin
    R = Jets.JetSSpace(Complex{Float64}, (8,4), (4,4), indexmap)
    @test size(R) == (8,4)
    @test eltype(R) == Complex{Float64}
    real.(ones(R)) ≈ real.(ones(Complex{Float64},8,4))
    imag.(ones(R)) ≈ imag.(ones(Complex{Float64},8,4))
    @test real.(zeros(R)) ≈ zeros(8,4)
    @test imag.(zeros(R)) ≈ zeros(8,4)
    @test size(rand(R)) == (8,4)
    @test size(Array(R)) == (8,4)
    @test eltype(Array(R)) == Complex{Float64}
    x = rand(R)
    z = similar(x)
    @test typeof(z) == Jets.SymmetricArray{Complex{Float64},2,typeof(indexmap)}
    y = x.A
    @test norm(x) ≈ sqrt(2*(norm(y)^2))
    @test norm(x,2) ≈ norm(x)
    @test norm(x,1) ≈ 2*norm(y,1)
    @test norm(x,Inf) ≈ norm(y,Inf)
    x[1,1] = x[6,1] = 0
    @test norm(x,0) ≈ 2*norm(y,0)
    @test space(rand(R)) == R
end

@testset "Symmetric spaces, broadcast" begin
    R = Jets.JetSSpace(Complex{Float64}, (8,4), (4,4), indexmap)
    u = rand(R)
    v = rand(R)
    w = rand(R)
    a = rand(Float64)
    b = rand(Float64)
    c = rand(Float64)
    x = a*u .+ b*v .+ c*w
    @test typeof(x) == Jets.SymmetricArray{Complex{Float64},2,typeof(indexmap)}

    for i2 = 1:4, i1 = 1:4
        @test x.A[i1,i2] ≈ a*u.A[i1,i2] + b*v.A[i1,i2] + c*w.A[i1,i2]
    end

    y = zeros(R)
    y .= x
    @test typeof(x) == Jets.SymmetricArray{Complex{Float64},2,typeof(indexmap)}

    for i2 = 1:4, i1 = 1:4
        @test y.A[i1,i2] ≈ a*u.A[i1,i2] + b*v.A[i1,i2] + c*w.A[i1,i2]
    end

    y .= x
    @test y ≈ x
    y .*= -1
    for i = 1:length(parent(y))
        parent(y)[i] ≈ -parent(x)[i]
    end

    z = abs.(x)
    @test typeof(z) == Array{Float64,2}
    for i = 1:length(z)
        @test z[i] ≈ abs(x[i])
    end
end

@testset "composition, linear" begin
    A₁,A₂,A₃,A₄ = map(d->JopBaz(rand(10,10)), 1:4)
    B₁,B₂,B₃,B₄ = map(A->state(A).A, (A₁,A₂,A₃,A₄))
    A₂₁ = A₂ ∘ A₁
    A₃₂₁ = A₃ ∘ A₂ ∘ A₁
    A₄₃₂₁ = A₄ ∘ A₃ ∘ A₂ ∘ A₁
    m = rand(domain(A₁))
    d = A₂₁*m
    @test d ≈ B₂ * ( B₁ * m )
    d = A₃₂₁*m
    @test d ≈ B₃ * (B₂ * ( B₁ * m))
    d = A₄₃₂₁*m
    @test d ≈ B₄ * (B₃ * ( B₂ * ( B₁ * m)))
    @test length(state(A₄₃₂₁).ops) == 4

    a = A₂₁'*d
    @test a ≈ (B₁' * ( B₂' * d))
    a = A₃₂₁'*d
    @test a ≈ (B₁' * ( B₂' * ( B₃' * d )))
    a = A₄₃₂₁'*d
    @test a ≈ (B₁' * ( B₂' * ( B₃' * ( B₄' * d))))

    @test domain(A₄₃₂₁) == JetSpace(Float64, 10)
    @test eltype(A₄₃₂₁) == Float64

    B₄₃₂₁ = convert(Array, A₄₃₂₁)
    @test B₄₃₂₁*m ≈ A₄₃₂₁*m

    C₄₃₂₁ = A₄ ∘ A₃ ∘ A₂₁'
    @test C₄₃₂₁*m ≈ A₄ * (A₃ * (A₁' * (A₂' * m)))
end

@testset "composition, linear with matrices" begin
    A₁,A₂,A₄ = map(d->JopBaz(rand(10,10)), 1:4)
    A₃ = rand(10,10)
    B₁,B₂,B₄ = map(A->state(A).A, (A₁,A₂,A₄))
    B₃ = A₃
    A₂₁ = A₂ ∘ A₁
    A₃₂₁ = A₃ ∘ A₂ ∘ A₁
    A₄₃₂₁ = A₄ ∘ A₃ ∘ A₂ ∘ A₁
    m = rand(domain(A₁))
    d = A₂₁*m
    @test d ≈ B₂ * ( B₁ * m )
    d = A₃₂₁*m
    @test d ≈ B₃ * (B₂ * ( B₁ * m))
    d = A₄₃₂₁*m
    @test d ≈ B₄ * (B₃ * ( B₂ * ( B₁ * m)))

    a = A₂₁'*d
    @test a ≈ (B₁' * ( B₂' * d))
    a = A₃₂₁'*d
    @test a ≈ (B₁' * ( B₂' * ( B₃' * d )))
    a = A₄₃₂₁'*d
    @test a ≈ (B₁' * ( B₂' * ( B₃' * ( B₄' * d))))

    @test domain(A₄₃₂₁) == JetSpace(Float64, 10)
    @test eltype(A₄₃₂₁) == Float64

    B₄₃₂₁ = convert(Array, A₄₃₂₁)
    @test B₄₃₂₁*m ≈ A₄₃₂₁*m
end

@testset "composition, nonlinear" begin
    F₁,F₂,F₃,F₄ = map(i->JopBar(10), 1:4)
    F₂₁ = F₂ ∘ F₁
    F₃₂₁ = F₃ ∘ F₂ ∘ F₁
    F₄₃₂₁ = F₄ ∘ F₃ ∘ F₂ ∘ F₁
    m = rand(domain(F₁))
    d = F₂₁ * m
    @test d ≈ F₂*(F₁*m)
    d = F₃₂₁*m
    @test d ≈ F₃*( F₂ * ( F₁ * m))
    d = F₄₃₂₁ * m
    @test d ≈ F₄ * ( F₃ * ( F₂ * ( F₁ * m)))

    m = ones(10)
    J₁ = jacobian!(F₁, m)
    J₂₁ = jacobian!(F₂, F₁*m) ∘ J₁
    J₃₂₁ = jacobian!(F₃, (F₂ ∘ F₁) * m) ∘ jacobian!(F₂, F₁ * m) ∘ J₁
    J₄₃₂₁ = jacobian!(F₄, (F₃ ∘ F₂ ∘ F₁) * m) ∘ jacobian!(F₃, (F₂ ∘ F₁) * m) ∘ jacobian!(F₂, F₁ * m) ∘ J₁

    L₁ = jacobian!(F₁, m)
    L₂₁ = jacobian!(F₂₁, m)
    L₃₂₁ = jacobian!(F₃₂₁, m)
    L₄₃₂₁ = jacobian!(F₄₃₂₁, m)

    δm = ones(10)
    @test J₁ * δm ≈ L₁ * δm
    @test J₂₁ * δm ≈ L₂₁ * δm
    @test J₃₂₁ * δm ≈ L₃₂₁ * δm
    @test J₄₃₂₁ * δm ≈ L₄₃₂₁ * δm

    δd = J₄₃₂₁*δm
    @test J₄₃₂₁'*δd ≈ L₄₃₂₁'*δd
end

@testset "composition, linear+nonlinear, adjoints" begin
    A₂ = JopBaz(rand(10,10))
    A₄ = JopFoo(rand(10))
    F₁,F₃ = map(i->JopBar(10), 1:2)
    F₂₁ = A₂' ∘ F₁
    F₃₂₁ = F₃ ∘ A₂' ∘ F₁
    F₄₃₂₁ = A₄ ∘ F₃ ∘ A₂' ∘ F₁
    m = rand(domain(F₁))
    d = F₂₁ * m
    @test d ≈ A₂'*(F₁*m)
    d = F₃₂₁*m
    @test d ≈ F₃*( A₂' * ( F₁ * m))
    d = F₄₃₂₁ * m
    @test d ≈ A₄ * ( F₃ * ( A₂' * ( F₁ * m)))

    m = rand(10)
    J₁ = jacobian!(F₁, m)
    J₂₁ = A₂' ∘ jacobian!(F₁, m)
    J₃₂₁ = jacobian!(F₃, A₂'*(F₁*m)) ∘ A₂' ∘ jacobian!(F₁, m)
    J₄₃₂₁ = A₄ ∘ jacobian!(F₃, A₂'*(F₁*m)) ∘ A₂' ∘ jacobian!(F₁, m)

    L₁ = jacobian!(F₁, m)
    L₂₁ = jacobian!(F₂₁, m)
    L₃₂₁ = jacobian!(F₃₂₁, m)
    L₄₃₂₁ = jacobian!(F₄₃₂₁, m)

    δm = rand(10)
    @test J₁ * δm ≈ L₁ * δm
    @test J₂₁ * δm ≈ L₂₁ * δm
    @test J₃₂₁ * δm ≈ L₃₂₁ * δm
    @test J₄₃₂₁ * δm ≈ L₄₃₂₁ * δm
end

@testset "composition operator * operator" begin
    A₁ = JopFoo(rand(2))
    A₂ = @blockop [JopBar(2), JopBar(2)]
    A = A₂ ∘ A₁

    A₁₁ = getblock(A, 1, 1)
    A₂₁ = getblock(A, 2, 1)

    m = rand(domain(A))
    @test getblock(A * m, 1) ≈ A₁₁ * m
    @test getblock(A * m, 2) ≈ A₂₁ * m
end

@testset "composition, state" begin
    diag = rand(2,2)
    A = JopFoo(diag)
    F = JopBar(2)
    G = A ∘ F
    @test state(G, :diagonal) ≈ state(A, :diagonal)
    @test state(G, :diagonal) ≈ state(A).diagonal

    @test_throws Exception state(G, :foo)

    B = JopFoo(diag)
    C = A ∘ B
    @test_throws Exception state(C, :digaonal)
end

@testset "sum, linear" begin
    A₁,A₂,A₃ = map(i->JopBaz(rand(10,10)), 1:3)
    A₁₂ = A₁ + A₂
    A₁₂₃ = A₁ + A₂ - A₃

    m = rand(domain(A₁))
    @test A₁₂ * m ≈ A₁*m + A₂*m
    @test A₁₂₃ * m ≈ A₁*m + A₂*m - A₃*m

    A₁₂₃ = A₁₂ + A₃
    @test A₁₂₃ * m ≈ A₁*m + A₂*m + A₃*m
    A₁₂₃₁₂ = A₁₂₃ - A₁₂
    @test A₁₂₃₁₂*m ≈ A₁*m + A₂*m + A₃*m - A₁*m - A₂*m

    d = rand(range(A₁))
    @test A₁₂₃' * d ≈ A₁'*d + A₂'*d + A₃'*d
end

@testset "sum, linear with compositions" begin
    A₁,A₂,A₃ = map(i->JopBaz(rand(10,10)), 1:3)
    a₁,a₂,a₃ = map(i->rand(), 1:3)
    A₁₂ = a₁*A₁ + a₂*A₂
    A₁₂₃ = a₁*A₁ + a₂*A₂ - a₃*A₃

    m = rand(domain(A₁))
    @test A₁₂ * m ≈ a₁*(A₁*m) + a₂*(A₂*m)
    @test A₁₂₃ * m ≈ a₁*(A₁*m) + a₂*(A₂*m) - a₃*(A₃*m)

    A₁₂₃ = A₁₂ + a₃*A₃
    @test A₁₂₃ * m ≈ a₁*(A₁*m) + a₂*(A₂*m) + a₃*(A₃*m)
    A₁₂₃₁₂ = A₁₂₃ - A₁₂
    @test A₁₂₃₁₂*m ≈ a₁*(A₁*m) + a₂*(A₂*m) + a₃*(A₃*m) - a₁*(A₁*m) - a₂*(A₂*m)

    d = rand(range(A₁))
    @test A₁₂₃' * d ≈ a₁*(A₁'*d) + a₂*(A₂'*d) + a₃*(A₃'*d)
end

@testset "sum, linear with matrices" begin
    A₁,A₃ = map(i->JopBaz(rand(10,10)), 1:3)
    A₂ = rand(10,10)
    A₁₂ = A₁ + A₂
    A₁₂₃ = A₁ + A₂ - A₃
    m = rand(domain(A₁))
    @test A₁₂ * m ≈ A₁*m + A₂*m
    @test A₁₂₃ * m ≈ A₁*m + A₂*m - A₃*m
end

@testset "sum, linear+nonlinar" begin
    A₁ = JopBaz(rand(10,10))
    F₂ = JopBar(10)

    F₁₂ = A₁ + F₂
    m = rand(domain(F₁₂))
    @test F₁₂*m ≈ A₁*m + F₂*m

    J₁₂ = jacobian(F₁₂,m)
    @test J₁₂*m ≈ A₁*m + jacobian(F₂,m)*m
end

@testset "block arrays" begin
    R = Jets.JetBSpace([JetSpace(Float64,2),JetSpace(Float64,2,2),JetSpace(Float64,2,3)])
    x = ones(R)
    @test getblock(x,1) ≈ ones(2)
    @test getblock(x,2) ≈ ones(2,2)
    @test getblock(x,3) ≈ ones(2,3)
    setblock!(x,1,π)
    setblock!(x,2,2π)
    setblock!(x,3,3*pi*ones(2,3))
    @test getblock!(x, 2, Array{Float64}(undef,2,2)) ≈ 2π*ones(2,2)
    _x = convert(Array, x)
    @test _x ≈ x
    @test norm(x) ≈ norm(_x)
    @test norm(x,0) ≈ norm(_x,0)
    @test norm(x,Inf) ≈ norm(_x,Inf)

    x = rand(R)
    _x = convert(Array, x)
    mn,mx = extrema(x)
    _mn,_mx = extrema(_x)
    @test mn ≈ _mn
    @test mx ≈ _mx

    R = Jets.JetBSpace([JetSpace(ComplexF64,2),JetSpace(ComplexF64,2,2),JetSpace(ComplexF64,2,3)])
    x = rand(R)
    _x = convert(Array, x)
    y = abs.(x)
    @test eltype(y) == Float64
    @test y ≈ abs.(_x)
    @test R == space(rand(R))
end

@testset "block arrays, broadcasting" begin
    R = Jets.JetBSpace([JetSpace(Float64,2),JetSpace(Float64,2,2),JetSpace(Float64,2,3)])
    u = rand(R)
    v = rand(R)
    w = rand(R)
    a = rand(Float64)
    b = rand(Float64)
    c = rand(Float64)
    x = a*u .+ b*v .+ c*w
    @test typeof(x) == Jets.BlockArray{Float64,Array{Float64}}

    for i = 1:length(x.arrays)
        for j = 1:length(x.arrays[i])
            @test x.arrays[i][j] ≈ a*u.arrays[i][j] + b*v.arrays[i][j] + c*w.arrays[i][j]
        end
    end

    y = zeros(R)
    y .= x
    @test typeof(y) == Jets.BlockArray{Float64,Array{Float64}}

    for i = 1:length(x.arrays)
        for j = 1:length(x.arrays[i])
            @test y.arrays[i][j] ≈ a*u.arrays[i][j] + b*v.arrays[i][j] + c*w.arrays[i][j]
        end
    end

    z = zeros(12)
    z .= x
    @test typeof(z) == Array{Float64,1}
    @test z ≈ x
    x .= 3.14
    for i = 1:length(x)
        @test x[i] ≈ 3.14
    end

    x = rand(Int32, size(R))
    y = rand(R)
    z = x .* y
    @test isa(z, Jets.BlockArray)
    @test convert(Array, z) ≈ x .* convert(Array, y)

    x = similar(y, Int)
    _x = similar(y, Int, length(y))
    __x = similar(y, Int, (length(y),))
    @test isa(_x, Jets.BlockArray)
    @test isa(__x, Jets.BlockArray)
end

@testset "block array, reshaped from array" begin
    x = rand(5,10)
    R = JetBSpace([JetSpace(Float64,5) for i=1:10])
    _x = reshape(x, R)
    for i = 1:length(_x)
        @test _x[i] ≈ x[i]
        _x[i] = i
        @test x[i] ≈ i
    end
    A = [rand(10,10) for i = 1:5]
    _A = @blockop [JopBaz(A[i]) for i = 1:5]

    m = rand(domain(_A))
    _y = _A*m
    for i = 1:5
        @test getblock(_y, i) ≈ A[i]*m
    end
end

@testset "block operator" begin
    B₁₁,B₁₃,B₁₄,B₂₁,B₂₃,B₂₄,B₃₂,B₃₃ = map(i->rand(10,10), 1:8)
    A₁₁,A₁₃,A₁₄,A₂₁,A₂₃,A₃₂,A₃₃ = map(B->JopBaz(B), (B₁₁,B₁₃,B₁₄,B₂₁,B₂₃,B₃₂,B₃₃))
    A₂₄ = (JopBaz(B₂₄))'
    F₁₂,F₂₃,F₃₁ = map(i->JopBar(10), 1:3)
    Z₂₂,Z₃₄ = map(i->JopZeroBlock(JetSpace(Float64,10), JetSpace(Float64,10)), 1:2)

    @test iszero(Z₂₂)
    @test !iszero(A₁₁)
    @test !iszero(F₁₂)

    C₂₄ = A₂₄ ∘ JopBar(10)

    G = [A₁₁ F₁₂ A₁₃ A₁₄;
         A₂₁ Z₂₂ F₂₃ C₂₄;
         F₃₁ A₃₂ A₃₃ Z₃₄]

    @test isa(@blockop(G), JopNl{<:Jet{<:Jets.JetBSpace,<:Jets.JetBSpace,typeof(Jets.JetBlock_f!)}})

    F = @blockop [A₁₁ F₁₂ A₁₃ A₁₄;
                  A₂₁ Z₂₂ F₂₃ C₂₄;
                  F₃₁ A₃₂ A₃₃ Z₃₄]

    @test nblocks(F) == (3,4)
    @test nblocks(F,1) == 3
    @test nblocks(F,2) == 4
    @test nblocks(A₁₁) == (1,1)

    _F₁₂ = getblock(JopNl,F,1,2)
    m = rand(domain(F₁₂))
    @test F₁₂*m ≈ _F₁₂*m

    _A₁₁ = getblock(JopLn,F,1,1)
    m = rand(domain(A₁₁))
    @test A₁₁*m ≈ _A₁₁*m

    @test ones(domain(F)) ≈ ones(40)
    @test zeros(range(F)) ≈ zeros(30)
    @test size(rand(range(F))) == (30,)

    m = rand(domain(F))
    d = F*m
    @test d[1:10]  ≈ B₁₁ * m[1:10] + F₁₂ * m[11:20] + B₁₃ * m[21:30] + B₁₄ * m[31:40]
    @test d[11:20] ≈ B₂₁ * m[1:10]                  + F₂₃ * m[21:30] + C₂₄ * m[31:40]
    @test d[21:30] ≈ F₃₁ * m[1:10] + B₃₂ * m[11:20] + B₃₃ * m[21:30]

    J = jacobian!(F,m)
    δm = rand(domain(J))
    δd = J * δm

    J₁₂ = jacobian!(F₁₂,m[11:20])
    J₂₃ = jacobian!(F₂₃,m[21:30])
    J₂₄ = jacobian!(C₂₄,m[31:40])
    J₃₁ = jacobian!(F₃₁,m[1:10])

    L = @blockop [A₁₁ J₁₂ A₁₃ A₁₄;
                  A₂₁ Z₂₂ J₂₃ J₂₄;
                  J₃₁ A₃₂ A₃₃ Z₃₄]

    @test δd ≈ L*δm
    @test L'*δd ≈ J'*δd

    @test mul!(rand(domain(L)), L', δd) ≈ J'*δd

    @test eltype(L) == Float64
    K = convert(Array, L)

    @test L*δm ≈ K*δm

    _J₁₂ = getblock(J, 1, 2)
    @test isa(J₁₂, JopLn)
    δm = rand(domain(J₁₂))
    @test J₁₂*δm ≈ _J₁₂*δm
end

@testset begin "block operator with keyword arguments"
    x = @blockop [JopBaz(rand(2,2)) ; JopBaz(rand(2,2))] foo = 3
    @test isa(x, Jop)
    x = @blockop [JopBaz(rand(2,2)) ; JopBaz(rand(2,2))] foo = 3 bar = 4
    @test isa(x, Jop)
end

@testset "block operator, singleton" begin
    B = JopBaz(rand(5,5))
    A = @blockop [B for i=1:1, j=1:1]
    m = rand(domain(A))
    @test A*m ≈ B*m
    d = rand(range(A))
    @test A'*d ≈ B'*d

    F = JopBar(5)
    G = @blockop [F for i=1:1, j=1:1]
    @test F*m ≈ G*m
    J = jacobian!(F, m)
    @test J*m ≈ jacobian!(G,m)*m
    @test J'*d ≈ (jacobian!(G,m)'*d)
end

@testset "block operator, tall-and-skinny" begin
    B = [JopBaz(rand(5,5)) for i=1:3, j=1:1]
    A = @blockop B
    m = rand(domain(A))
    @test A*m ≈ [B[1,1]*m;B[2,1]*m;B[3,1]*m]
    d = rand(range(A))
    @test A'*d ≈ B[1,1]'*d[1:5] + B[2,1]'*d[6:10] + B[3,1]'*d[11:15]

    G = [JopBar(5) for i=1:3, j=1:1]
    F = @blockop G
    @test F*m ≈ [G[1,1]*m;G[2,1]*m;G[3,1]*m]
    J = jacobian!(F,m)
    @test J*m ≈ [jacobian!(G[1,1],m)*m;jacobian!(G[2,1],m)*m;jacobian!(G[3,1],m)*m]
    @test J'*d ≈ jacobian!(G[1,1],m)'*d[1:5] + jacobian!(G[2,1],m)'*d[6:10] + jacobian!(G[3,1],m)'*d[11:15]

    B₁₁ = getblock(A,1,1)
    @test isa(B₁₁,JopLn)
    @test B₁₁*m ≈ B[1,1]*m

    @test nblocks(A) == (3,1)
    @test nblocks(A, 1) == 3
    @test nblocks(A, 2) == 1
end

@testset "block operator, short-and-fat" begin
    B = [JopBaz(rand(5,5)) for i=1:1, j=1:3]
    A = @blockop B
    m = rand(domain(A))
    @test A*m ≈ B[1,1]*m[1:5] + B[1,2]*m[6:10] + B[1,3]*m[11:15]
    d = rand(range(A))
    @test A'*d ≈ [B[1,1]'*d;B[1,2]'*d;B[1,3]'*d]

    G = [JopBar(5) for i=1:1, j=1:3]
    F = @blockop G
    @test F*m ≈ G[1,1]*m[1:5] + G[1,2]*m[6:10] + G[1,3]*m[11:15]
    J = jacobian!(F,m)
    @test J*m ≈ jacobian!(G[1,1],m[1:5])*m[1:5] + jacobian!(G[1,2],m[6:10])*m[6:10] + jacobian!(G[1,3],m[11:15])*m[11:15]
    @test J'*d ≈ [jacobian!(G[1,1],m[1:5])'*d;jacobian!(G[1,2],m[6:10])'*d;jacobian!(G[1,3],m[11:15])'*d]
end

@testset "getblock, adjoint operator" begin
    B = [JopBaz(rand(5,5)) for i=1:2, j=1:3]
    A = @blockop B
    C = A'

    for i = 1:2, j = 1:3
        Cji = getblock(C, j, i)
        @test isa(Cji, JopAdjoint)
        @test state(Cji).A ≈ state(B[i,j]).A
        m = rand(domain(Cji))
        @test Cji*m ≈ B[i,j]'*m
    end

    B = [JopBaz(rand(5,5))' for i=1:2, j=1:3]
    A = @blockop B
    C = A'

    C11 = getblock(C,1,1)
    isa(C11,JopLn)

    for i = 1:2, j=1:3
        Cji = getblock(C, j, i)
        @test isa(Cji, JopLn)
        @test state(Cji).A ≈ state(B[i,j]).A
        m = rand(domain(Cji))
        @test Cji*m ≈ B[i,j]'*m
    end
end

@testset "scalar multplying an operator" begin
    A = JopBaz(rand(10,10))
    a = 3.14
    m = rand(domain(A))
    B = a*A
    @test B*m ≈ a*(A*m)
end

@testset "vectorized operator" begin
    A = JopFoo2(rand(10,11))
    x = rand(domain(A))

    @test size(vec(domain(A))) == (10*11,)

    B = vec(A)
    @test isa(jet(B), Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(Jets.JetVec_f!)})

    d = A * x
    _d = vec(A) * vec(x)
    __d = reshape(_d, range(A))
    @test d[:] ≈ _d
    @test __d ≈ d
    @test size(d) == (10,11)
    @test size(_d) == (10*11,)
    @test size(__d) == (10,11)

    C = JopFoo(rand(10,11))
    @test isa(jet(C), Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(JopFoo_df!)})
end

@testset "vectorized block operator" begin
    A = @blockop [JopFoo2(rand(10,11)), JopFoo2(rand(10,11))]
    x = rand(domain(A))

    @test size(vec(domain(A))) == (10*11,)
    @test size(vec(range(A))) == (2*10*11,)
    @test size(range(A)) == size(vec(range(A)))

    B = vec(A)
    @test isa(jet(B), Jet{<:JetAbstractSpace,<:JetAbstractSpace,typeof(Jets.JetVec_f!)})

    d = A * x
    _d = vec(A) * vec(x)
    @test d ≈ _d

    a = A' * d
    _a = vec(A') * d
    __a = reshape(a, domain(A))
    @test a ≈ __a
end

@testset "close op" begin
    A = JopClose(rand(2))
    @test isfile(state(A).file)
    close(A)
    @test !isfile(state(A).file)
end

@testset "close block op" begin
    A = @blockop [JopClose(rand(2)) JopClose(rand(2)) ; JopClose(rand(2)) JopClose(rand(2))]
    A₁₁ = getblock(A,1,1)
    A₁₂ = getblock(A,1,2)
    A₂₁ = getblock(A,2,1)
    A₂₂ = getblock(A,2,2)
    for B in (A₁₁, A₁₂, A₂₁, A₂₂)
        @test isfile(state(B).file)
    end
    close(A)
    for B in (A₁₁, A₁₂, A₂₁, A₂₂)
        @test !isfile(state(B).file)
    end
end

@testset "close composite op" begin
    A₁ = JopClose(rand(2))
    A₂ = JopClose(rand(2))
    A = A₂ ∘ A₁
    for B in (A₁, A₂)
        @test isfile(state(B).file)
    end
    close(A)
    for B in (A₁, A₂)
        @test !isfile(state(B).file)
    end
end

@testset "perfstat" begin
    A₁ = JopFoo(rand(2))
    A₂ = JopBar(2)
    A = A₂ ∘ A₁

    @test perfstat(A₁) ≈ Float64(π)
    @test perfstat(A₂) == nothing
    @test perfstat(A) ≈ Float64(π)

    B = A₂ + A₁
    @test perfstat(B) ≈ Float64(π)
end


text:
using Dates

"""
    filter_by_ids( selected_ids::Array{<:Integer}, ids::Array{<:Integer})

Checks for matching IDs.
"""
function filter_by_ids( selected_ids::Array{<:Integer}, ids::Array{<:Integer})

    # look for the array positions where the IDs match
    matched = sort( get_index_list( selected_ids, ids ) )

    return matched
end

"""
    read_filtered(snap_file, blockname, parttype, block_position, matched)

Helper function to return the correct filtered arrays.
"""
function read_filtered(snap_file, blockname, parttype, block_position, matched)
    data = read_block(snap_file, blockname, parttype=parttype, block_position=block_position)
    if size(data,2) > 1
        return data[:,matched]
    else
        return data[matched]
    end
end

"""
    read_particles_by_id_single_file(snap_file::String, halo_ids::Array{<:Integer},
                                          blocks::Array{String}, parttype::Integer)

Helper function to read the matching IDs of particles in one whole file.
"""
function read_particles_by_id_single_file(snap_file::String, halo_ids::Array{<:Integer},
                                          blocks::Array{String}, parttype::Integer;
                                          verbose::Bool=true)

    # read ids in file
    ids = read_snap(snap_file, "ID", parttype)

    if verbose
        @info "Matching IDs..."
        t1 = Dates.now()
    end

    # check which IDs match
    matched = filter_by_ids( halo_ids, ids)

    if verbose
        t2 = Dates.now()
        @info "Found $(size(matched,1)) matches. Took: $(t2 - t1)"
    end

    # store block positions for faster read-in
    block_positions = get_block_positions(snap_file)

    # read the data blocks whole, but only store relevant entries
    return Dict(blocks[i] => read_filtered(snap_file, blocks[i], parttype, block_positions[blocks[i]], matched)
                for i = 1:size(blocks,1))

end

"""
    construct_matched_dict(data_in::Dict{String, Union{Vector{T}, Array{T,2}}},
                                blocks::Array{String,1}, matched::Array{<:Integer,1}) where T

Write all matching particles to a new dictionary.
"""
function construct_matched_dict(data_in::Dict{String, VecOrMat{T} where T},
                                blocks::Array{String,1}, matched::Array{<:Integer,1})

    dict_out = Dict{String, VecOrMat{T} where T}()

    for block ∈ blocks
        dim   = size(data_in[block], 2)
        if dim == 1
            dict_out[block] = data_in[block][matched]
        else
            dict_out[block] = data_in[block][:, matched]
        end
    end
    return dict_out
end

"""
    read_particles_by_id(snap_base::String, ids::Array{<:Integer},
                         blocks::Array{String};
                         parttype::Integer=0, verbose::Bool=true,
                         pos0::Array{<:Real}=[-1.234, -1.234, -1.234],
                         r0::Real=0.0)

Reads particles filtered by the provided IDs. Returns all requested blocks as entries in a `Dict`.
"""
function read_particles_by_id(snap_base::String, selected_ids::Array{<:Integer},
                              blocks::Array{String};
                              parttype::Integer=0, verbose::Bool=true,
                              pos0::Union{Array{<:Real},Nothing}=nothing,
                              r0::Real=0.0,
                              use_keys::Bool=true)

    # try reading the first of the distributed snapshots
    filename = select_file(snap_base, 0)

    # check if blocks are present
    blocks, no_mass_block = check_blocks(filename, blocks)

    # extend the list of blocks to read by ID block
    blocks = [ blocks ; "ID" ]

    # remove duplicate blocks
    unique!(blocks)

    # for a given halo position and search radius we can use `read_particles_in_volume`
    if ( !isnothing(pos0) && r0 != 0.0)

        # check if .key files are present
        if isfile(filename * ".key")

            # read all particles in the defined volume
            data = read_particles_in_volume(snap_base, blocks, pos0, r0,
                                            parttype=parttype, verbose=verbose,
                                            use_keys=use_keys)

            if verbose
                println()
                @info "Matching IDs..."
                t1 = Dates.now()
            end

            # find matching entries
            matched = filter_by_ids( selected_ids, data["ID"])

            if verbose
                t2 = Dates.now()
                @info "Found $(size(matched,1)) matches. Took: $(t2 - t1)"
            end

            # construct a new dict only containing the matching particles
            return construct_matched_dict(data, blocks, matched)

        else # if there are no .key files we need to read the while snapshot
            return read_particles_by_id_single_file(filename, selected_ids, blocks, parttype, verbose=verbose)
        end

    else # otherwise we need to search whole files

        # single file read
        if isfile(snap_base)

            return read_particles_by_id_single_file(snap_base, selected_ids, blocks, parttype, verbose=verbose)

        else # multi-file brute-force read -> slow!

            # total number of particles to read
            N_to_read = size(selected_ids,1)

            # number of particles read so far
            N_read = 0

            # read the info block
            snap_info = read_info(filename)

            # prepare data dict for storage
            data = allocate_data_dict(blocks, N_to_read, snap_info, no_mass_block)

            filename = snap_base * ".0"
            h = head_to_obj(filename)

            @warn "Brute-force reading $(h.num_files) files! This may take a while!"

            # loop over all the files in the snap directory
            for i = 0:(h.num_files-1)

                # get current file name
                filename = select_file(snap_base, i)

                # read data from file
                data_file = read_particles_by_id_single_file(filename, selected_ids, blocks, parttype, verbose=verbose)

                N_this_file = size(data_file["ID"],1)

                # write into master dict
                for block ∈ blocks
                    dim   = snap_info[getfield.(snap_info, :block_name) .== block][1].n_dim
                    if dim == 1
                        data[block][N_read+1:N_read+N_this_file]    = data_file[block]
                    else
                        data[block][:, N_read+1:N_read+N_this_file] = data_file[block]
                    end
                end # blocks

                # update number of read particles
                N_read += N_this_file

            end

            # reduce array size
            for block ∈ blocks
                # eltype needed for type stable reshape function call
                dim   = eltype(N_read)(snap_info[getfield.(snap_info, :block_name) .== block][1].n_dim)
                if dim == 1
                    resize!(data[block], N_read)
                else
                    data[block] = reshape(resize!(vec(data[block]), dim*N_read), dim, N_read)
                end
            end # blocks

            return data

        end
    end
end

"""
    read_particles_by_id(snap_base::String, ids::Array{<:Integer},
                         block::String;
                         parttype::Integer=0, verbose::Bool=true,
                         pos0::Array{<:Real}=[-1.234, -1.234, -1.234],
                         r0::Real=0.0)

Reads particles filtered by the provided IDs. Returns the requested block as an `Array`.
"""
function read_particles_by_id(snap_base::String, selected_ids::Array{<:Integer},
                              block::String;
                              parttype::Integer=0, verbose::Bool=true,
                              pos0::Union{Array{<:Real},Nothing}=nothing,
                              r0::Real=0.0,
                              use_keys::Bool=true)

   data = read_particles_by_id(snap_base, selected_ids,
                              [block],
                              parttype=parttype, verbose=verbose,
                              pos0=pos0, r0=r0, use_keys=use_keys)

    return data[block]
end

ext: jl


text:
# ---
# title: 994. Rotting Oranges
# id: problem994
# author: Indigo
# date: 2021-02-18
# difficulty: Medium
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/rotting-oranges/description/>
# hidden: true
# ---
#
# In a given grid, each cell can have one of three values:
#
#   * the value `0` representing an empty cell;
#   * the value `1` representing a fresh orange;
#   * the value `2` representing a rotten orange.
#
# Every minute, any fresh orange that is adjacent (4-directionally) to a rotten
# orange becomes rotten.
#
# Return the minimum number of minutes that must elapse until no cell has a
# fresh orange.  If this is impossible, return `-1` instead.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/16/oranges.png)**
#
#
#
#     Input: [[2,1,1],[1,1,0],[0,1,1]]
#     Output: 4
#
#
# **Example 2:**
#
#
#
#     Input: [[2,1,1],[0,1,1],[1,0,1]]
#     Output: -1
#     Explanation: The orange in the bottom left corner (row 2, column 0) is never rotten, because rotting only happens 4-directionally.
#
#
# **Example 3:**
#
#
#
#     Input: [[0,2]]
#     Output: 0
#     Explanation: Since there are already no fresh oranges at minute 0, the answer is just 0.
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length <= 10`
#   2. `1 <= grid[0].length <= 10`
#   3. `grid[i][j]` is only `0`, `1`, or `2`.
#
#
## @lc code=start
using LeetCode, DataStructures

function oranges_rotting(grid::Vector{Vector{Int}})
    nr, nc = length(grid), length(grid[1])
    visited = fill(false, nr, nc)
    q = Queue{Tuple{Int,Int,Int}}()
    left, res = 0, 0
    for r in 1:nr, c in 1:nc
        if grid[r][c] == 2
            enqueue!(q, (r, c, 0))
        elseif grid[r][c] == 1
            left += 1
        end
    end
    dirs = [(0, 1), (1, 0), (0, -1), (-1, 0)]
    while !isempty(q)
        tp = dequeue!(q)
        for dir in dirs
            nb = (tp[1] + dir[1], tp[2] + dir[2], tp[3] + 1)
            if 0 < nb[1] ≤ nr &&
               0 < nb[2] ≤ nc &&
               !visited[nb[1], nb[2]] &&
               grid[nb[1]][nb[2]] == 1
                left -= 1
                visited[nb[1], nb[2]] = true
                enqueue!(q, nb)
            end
        end
        res = tp[3]
    end
    left == 0 ? res : -1
end
## @lc code=end


text:
using Unitful: @u_str
using UnitfulAtomic: austrip, auconvert

using .Calculators: AbstractCalculator, Calculator
using NQCModels: Model

abstract type AbstractSimulation{M} end

Base.broadcastable(sim::AbstractSimulation) = Ref(sim)

struct Simulation{M,Calc,T,Ttype,C} <: AbstractSimulation{M}
    temperature::Ttype
    cell::C
    atoms::Atoms{T}
    calculator::Calc
    method::M
end

"""
    Simulation(atoms::Atoms{T}, model::Model, method::M;
        temperature=0u"K", cell::AbstractCell=InfiniteCell()) where {M,S,T}

Simulation parameters that controls the types of atoms, interactions,
dynamics method, temperature and simulation cell.
"""
function Simulation(atoms::Atoms{T}, model::Model, method::M;
        temperature=0u"K", cell::AbstractCell=InfiniteCell()) where {M,T}
    calc = Calculator(model, length(atoms), T)
    Simulation(temperature, cell, atoms, calc, method)
end

struct RingPolymerSimulation{M,Calc,T,Ttype,C,B} <: AbstractSimulation{M}
    temperature::Ttype
    cell::C
    atoms::Atoms{T}
    calculator::Calc
    method::M
    beads::B
end

function RingPolymerSimulation(temperature, cell::AbstractCell,
        atoms::Atoms{T}, model::Model, method::M,
        n_beads::Integer, quantum_nuclei::Vector{Symbol}) where {M,T}

    if isempty(quantum_nuclei)
        beads = RingPolymers.RingPolymerParameters{T}(n_beads, get_temperature(temperature), length(atoms))
    else
        beads = RingPolymers.RingPolymerParameters{T}(n_beads, get_temperature(temperature), atoms.types, quantum_nuclei)
    end

    calc = Calculator(model, length(atoms), n_beads, T)
    RingPolymerSimulation(temperature, cell, atoms, calc, method, beads)
end

function RingPolymerSimulation(atoms::Atoms, model::Model, method::M, n_beads::Integer;
        temperature=0,
        cell::AbstractCell=InfiniteCell(), quantum_nuclei::Vector{Symbol}=Symbol[]) where {M}
    RingPolymerSimulation(temperature, cell, atoms, model, method, n_beads, quantum_nuclei)
end

NQCModels.nstates(sim::AbstractSimulation) = NQCModels.nstates(sim.calculator)
NQCModels.eachstate(sim::AbstractSimulation) = NQCModels.eachstate(sim.calculator)
NQCModels.nelectrons(sim::AbstractSimulation) = NQCModels.nelectrons(sim.calculator)
NQCModels.eachelectron(sim::AbstractSimulation) = NQCModels.eachelectron(sim.calculator)
NQCModels.mobileatoms(sim::AbstractSimulation) = NQCModels.mobileatoms(sim.calculator)
NQCModels.dofs(sim::AbstractSimulation) = NQCModels.dofs(sim.calculator)

NQCModels.ndofs(sim::AbstractSimulation) = NQCModels.ndofs(sim.calculator.model)
natoms(sim::AbstractSimulation) = length(sim.atoms)
RingPolymers.nbeads(sim::RingPolymerSimulation) = RingPolymers.nbeads(sim.beads)
masses(sim::AbstractSimulation) = sim.atoms.masses
masses(sim::AbstractSimulation, i::Int) = sim.atoms.masses[i]

function masses(sim::AbstractSimulation, I::CartesianIndex)
    masses(sim, I[2])
end

Base.size(sim::Simulation) = (ndofs(sim), natoms(sim))
Base.size(sim::RingPolymerSimulation) = (ndofs(sim), natoms(sim), RingPolymers.nbeads(sim))

function get_temperature(sim::Simulation, t::Real=0)
    t = auconvert(u"fs", t)
    get_temperature(sim.temperature, t)
end
function get_temperature(sim::RingPolymerSimulation, t::Real=0)
    t = auconvert(u"fs", t)
    get_temperature(sim.temperature, t) * length(sim.beads)
end
get_temperature(temperature::Number, t=0) = austrip(temperature)
get_temperature(temperature::Function, t=0) = austrip(temperature(t))

function Base.show(io::IO, sim::Simulation{M}) where {M}
    print(io, "Simulation{$M}:\n  ", sim.atoms, "\n  ", sim.calculator.model)
end

function Base.show(io::IO, sim::RingPolymerSimulation{M}) where {M}
    print(io, "RingPolymerSimulation{$M}:\n\n  ", sim.atoms, "\n\n  ", sim.calculator.model,
          "\n  with ", length(sim.beads), " beads.")
end


text:
"""
Example

```
using BayesOpt
config = ConfigParameters()         # calls initialize_parameters_to_default of the C API
set_kernel!(config, "kMaternARD5")  # calls set_kernel of the C API
config.sc_type = SC_MAP
f(x) = sum(x .^ 2)
lowerbound = [-2., -2.]; upperbound = [2., 2.]
optimizer, optimum = bayes_optimization(f, lowerbound, upperbound, config)
```

Exports: `KernelParameters`, `MeanParameters`, `ConfigParameters`, `bayes_optimization`,
`bayes_optimization_disc`, `bayes_optimization_categorical`, `set_kernel!`, `set_mean!`,
`set_criteria!`, `set_surrogate!`, `set_log_file!`, `set_load_file!`, `set_save_file!`,
`set_learning!`, `set_score!`, `L_FIXED`, `L_EMPIRICAL`, `L_DISCRETE`, `L_MCMC`, `L_ERROR`,
`SC_MTL`, `SC_ML`, `SC_MAP`, `SC_LOOCV`, `SC_ERROR`

See also https://rmcantin.bitbucket.io/html/usemanual.html
"""
module BayesOpt

const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("BayesOpt not installed properly, run Pkg.build(\"BayesOpt\"), restart Julia and try again")
end
include(depsjl_path)

# Module initialization function
function __init__()
    # Always check your dependencies from `deps.jl`
    check_deps()
end

export KernelParameters, MeanParameters, ConfigParameters, bayes_optimization,
bayes_optimization_disc, bayes_optimization_categorical

import Base: show
struct KernelParameters
    name::Cstring
    hp_mean::NTuple{128, Cdouble}
    hp_std::NTuple{128, Cdouble}
    n_hp::Csize_t
end
function show(io::IO, mime::MIME"text/plain", o::KernelParameters)
    println(io, "$(unsafe_string(o.name)) Kernel")
    if o.n_hp > 0
        println(io, "  hyperparameters mean $(o.hp_mean[1:o.n_hp])")
        println(io, "  hyperparameters std $(o.hp_std[1:o.n_hp])")
    end
end

@enum LearningType::Cint begin
    L_FIXED
    L_EMPIRICAL
    L_DISCRETE
    L_MCMC
    L_ERROR = -1
end
export L_FIXED, L_EMPIRICAL, L_DISCRETE, L_MCMC, L_ERROR

@enum ScoreType::Cint begin
    SC_MTL
    SC_ML
    SC_MAP
    SC_LOOCV
    SC_ERROR = -1
end
export SC_MTL, SC_ML, SC_MAP, SC_LOOCV, SC_ERROR

struct MeanParameters
    name::Cstring
    coef_mean::NTuple{128, Cdouble}
    coef_std::NTuple{128, Cdouble}
    n_coef::Csize_t
end
function show(io::IO, mime::MIME"text/plain", o::MeanParameters)
    println(io, "$(unsafe_string(o.name)) Mean")
    if o.n_coef > 0
        println(io, "  coefficients mean $(o.coef_mean[1:o.n_coef])")
        println(io, "  coefficients std $(o.coef_std[1:o.n_coef])")
    end
end

mutable struct ConfigParameters
    n_iterations::Csize_t
    n_inner_iterations::Csize_t
    n_init_samples::Csize_t
    n_iter_relearn::Csize_t
    init_method::Csize_t
    random_seed::Cint
    verbose_level::Cint
    log_filename::Cstring
    load_save_flag::Csize_t
    load_filename::Cstring
    save_filename::Cstring
    surr_name::Cstring
    sigma_s::Cdouble
    noise::Cdouble
    alpha::Cdouble
    beta::Cdouble
    sc_type::ScoreType
    l_type::LearningType
    l_all::Cint
    epsilon::Cdouble
    force_jump::Csize_t
    kernel::KernelParameters
    mean::MeanParameters
    crit_name::Cstring
    crit_params::NTuple{128, Cdouble}
    n_crit_params::Csize_t
end
"""
    ConfigParameters()

Returns default parameters of BayesOpt (see initialize_parameters_to_default in the C API).
"""
ConfigParameters() = ccall((:initialize_parameters_to_default, libbayesopt), ConfigParameters, ())

function show(io::IO, mime::MIME"text/plain", o::ConfigParameters)
    println(io, "ConfigParameters")
    for field in fieldnames(ConfigParameters)
        val = getfield(o, field)
        if field == :crit_params
            o.n_crit_params == 0 && continue
            println(io, "$field = $(val[1:o.n_crit_params])")
        end
        valtoshow = typeof(val) == Cstring ? unsafe_string(val) : val
        if typeof(val) == KernelParameters || typeof(val) == MeanParameters
            show(io, mime, val)
        else
            println(io, "$field = $valtoshow")
        end
    end
end

for func in [:set_kernel, :set_mean, :set_criteria, :set_surrogate, :set_log_file,
             :set_load_file, :set_save_file, :set_learning, :set_score]
    jfuncname = Symbol(func, "!")
    jfuncstring = split(string(jfuncname), ".")[end]
    @eval begin
        $jfuncname(config, name) = ccall(($(string(func)), libbayesopt), Cvoid, (Ptr{ConfigParameters}, Cstring), Ref(config), name)
        @doc "`$($jfuncstring)(config, name)` $(replace($(string(func)), "_" => " ")) in `config` to `name`." $jfuncname
        export $jfuncname
    end
end

@inline function prepare_cargs(func, n)
    optimizer = zeros(n); optimum  = Ref{Cdouble}(0)
    bofunc = (n, x, g, d) -> func(unsafe_wrap(Array, x, n))
    cfunc = @cfunction $bofunc Cdouble (Cuint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid})
    cfunc, optimizer, optimum
end

"""
    bayes_optimization(func, lb, ub, config)

Runs continuous Bayesian optimization on `func` that takes vectors of length `d`
as argument and returns a real number, within the box defined by the lowerbounds
`lb` (a vector of length `d` with lowerbounds for each dimension) and upperbounds
`ub`, using `config` (see `ConfigParameters`).
"""
function bayes_optimization(func, lb, ub, config)
    n = length(lb)
    length(ub) == n || @error("lowerbounds and upperbounds have different length.")
    cfunc, optimizer, optimum = prepare_cargs(func, n)
    ccall((:bayes_optimization, libbayesopt), Cint,
          (Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble}, ConfigParameters),
          n, cfunc, Ptr{Nothing}(0),
          pointer(lb), pointer(ub), pointer(optimizer), optimum, config)
    (optimizer = optimizer, optimum = optimum.x)
end

"""
    bayes_optimization_disc(func, valid_x, config)

Runs  Bayesian optimization on `func` that takes vectors of length `d`
as argument and returns a real number, on `valid_x` an array of valid points
(vectors of length `d`) using `config` (see `ConfigParameters`).
"""
function bayes_optimization_disc(func, xs::Array{<:Array{<:Number, 1}, 1}, config)
    n_points = length(xs)
    n = length(xs[1])
    valid_x = pointer(vcat(xs...))
    cfunc, optimizer, optimum = prepare_cargs(func, n)
    ccall((:bayes_optimization_disc, libbayesopt), Cint,
          (Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cdouble}, Csize_t,
           Ptr{Cdouble}, Ptr{Cdouble}, ConfigParameters),
          n, cfunc, Ptr{Nothing}(0),
          valid_x, Csize_t(n_points), pointer(optimizer), optimum, config)
    (optimizer = optimizer, optimum = optimum.x)
end

"""
    bayes_optimization_categorical(func, categories, config)

Runs Bayesian optimization on `func` that takes vectors of length `d` as
argument and returns a real number, with `categories` array of size `d` with the
number of categories per dimension, using `config` (see `ConfigParameters`).
"""
function bayes_optimization_categorical(func, categories::Array{Integer, 1}, config)
    n = length(categories)
    cfunc, optimizer, optimum = prepare_cargs(func, n)
    ccall((:bayes_optimization_categorical, libbayesopt), Cint,
          (Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cint},
           Ptr{Cdouble}, Ptr{Cdouble}, ConfigParameters),
          n, cfunc, Ptr{Nothing}(0),
          pointer(categories), pointer(optimizer), optimum, config)
    (optimizer = optimizer, optimum = optimum.x)
end
end


text:
# code for Video 4.4:
fname = "iris.csv"
using DataFrames
data = readtable(fname, separator = ',')

###Distributions:
Pkg.add("Distributions")
using Distributions

typeof(Normal) # DataType
super(Normal) # Distributions.Distribution{Distributions.Univariate,Distributions.Continuous}
super(Poisson) # Distributions.Distribution{Distributions.Univariate,Distributions.Discrete}

n1 = Normal(10, 1.5) # Distributions.Normal(μ=10.0, σ=1.5)
params(n1) # (10.0,1.5)
fieldnames(n1) # 2-element Array{Symbol,1}:  :μ :σ

rand(n1, 100)
# 100-element Array{Float64,1}:
#   8.48529
#   7.90919
#  13.2968
#   6.55532
#   ...
 # 13.1982
 # 10.7262
 # 11.1575
 # 11.5623

randn(50)
50-element Array{Float64,1}:
  1.34669
 -1.73781
  0.801373
 -1.11795
 ...

b = Binomial(10, 0.7) # binomial distribution with 10 trials and succes rate 0.7
# Distributions.Binomial(n=10, p=0.7)
rand(b, 1000)
1000-element Array{Int64,1}:
 6
 9
 7
 7
 6
 6
 9

### Kernel density:
using Winston, KernelDensity
x = convert(Array, data[1:50, :sepal_width])
y = convert(Array, data[1:50, :sepal_length])
kx = kde(x) # kernel density function
ky = kde(y)
plot(kx.x, kx.density, "r--", ky.x, ky.density, "b;")

###
using Gadfly
set_default_plot_size(20cm, 12cm)
Gadfly.plot(data, x = "petal_length", y = "petal_width", color = "species", Geom.point)

### K-means Clustering
Pkg.add("Clustering")
using Clustering
features = convert(Array, data[:, 1:4])' # 4x150 Array{Float64,2}:
# group the data onto 3 clusters
initseeds(:rand, features, 3) # choose 3 random starting points (centroids)
# 3-element Array{Int64,1}:
#  36
#  93
#  41
result = kmeans( features, 3) # make 3 clusters
# Clustering.KmeansResult{Float64}(4x3 Array{Float64,2}:
#  5.006  5.90161  6.85
#  3.418  2.74839  3.07368
#  1.464  4.39355  5.74211
#  0.244  1.43387  2.07105,[1,1,1,1,1,1,1,1,1,1  …  3,3,
Gadfly.plot(data, x = "petal_length", y = "petal_width", color = result.assignments, Geom.point)

### HypothesisTests
Pkg.add("HypothesisTests")
using HypothesisTests
# test sepal_length I. setosa versus I. versicolor:
xset = convert(Array, data[1:50, :sepal_length])
xver = convert(Array, data[51:100, :sepal_length])

uneqtt = UnequalVarianceTTest(xset, xver)
# Two sample t-test (unequal variance)
# ------------------------------------
# Population details:
#     parameter of interest:   Mean difference
#     value under h_0:         0
#     point estimate:          -0.9299999999999997
#     95% confidence interval: (-1.1057073727210804,-0.7542926272789189)

# Test summary:
#     outcome with 95% confidence: reject h_0
#     two-sided p-value:           3.74674261398387e-17 (extremely significant)

# Details:
#     number of observations:   [50,50]
#     t-statistic:              -10.52098626754911
#     degrees of freedom:       86.53800179765493
#     empirical standard error: 0.08839475466938762

pvalue(uneqtt) # 3.74674261398387e-17
ci(uneqtt) # (-1.1057073727210804,-0.7542926272789189)

max_stars_repo_path: Section 4/4_4.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 4/4_4.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 4/4_4.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
using Distributions: Distributions
using UnitfulAtomic: austrip
using Random: Random

struct BoltzmannVelocityDistribution{T,S} <: VelocityDistribution
    dist::Distributions.MvNormal{T}
    size::NTuple{S,Int}
end

"""
    BoltzmannVelocityDistribution(temperature, masses)

Obtain velocities sampled from the Boltzmann distribution for given temperature and masses.

This is effectively just a multivariate normal distribution but when combined with the
[`DynamicalDistribution`](@ref) it ensure samples are generated with the correct size.
If using for ring polymer dynamics, ensure the temperature is given as nbeads*temperature.
"""
function BoltzmannVelocityDistribution(temperature, masses, size)
    dist = Distributions.MvNormal(sqrt.(austrip(temperature) ./ masses))
    BoltzmannVelocityDistribution(dist, size)
end

function Base.rand(rng::AbstractRNG, s::BoltzmannVelocityDistribution{T,S}) where {T,S}
    out = Array{T,S}(undef, s.size...)
    if S == 2
        @views for i=1:size(out, 1)
            Random.rand!(rng, s.dist, out[i,:])
        end
    elseif S == 3
        @views for i=1:size(out, 3)
            for j=1:size(out, 1)
                Random.rand!(rng, s.dist, out[j,:,i])
            end
        end
    end
    return out
end

select_item(x::BoltzmannVelocityDistribution, ::Integer, ::NTuple) = rand(x)


text:
include("../neat_test.jl")

println("--> Dataset 5")
config = "test_dataset5_config.txt"
p5 = run_NEAT_epochs(config, 2000, false)
print_best_solution_perfomance(p5)

#=
plot_neat_fitness(p5, imagename="fitness-dt5.png")
plot_neat_IP(p5, imagename="IP-dt5.png")
plot_neat_species(p5, imagename="species-dt5.png")
plot_output_distribution(p5, imagename="distribution-dt5.png")
plot_complexity(p5, imagename="solution-complexity-dt5.png")
plot_hit_rate(p5, imagename="hit-rate-dt5.png")
=#
#pop = getPopulation("checkpoint-dt5.jld2")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
using ColPack
using LinearAlgebra
using MatrixMarket
using Random
using SparseArrays
using SparseDiffTools
using Test

orderings = Vector{AbstractOrdering}()

push!(orderings, natural_ordering())
push!(orderings, largest_first_ordering())
push!(orderings, dynamic_largest_first_ordering())
push!(orderings, distance_two_largest_first_ordering())
push!(orderings, smallest_last_ordering())
push!(orderings, distance_two_smallest_last_ordering())
push!(orderings, incidence_degree_ordering())
push!(orderings, distance_two_incidence_degree_ordering())
# push!(orderings, "RANDOM")

ncolors = Vector{Int}()

push!(ncolors, 7)
push!(ncolors, 7)
push!(ncolors, 6)
push!(ncolors, 6)
push!(ncolors, 6)
push!(ncolors, 6)
push!(ncolors, 7)
push!(ncolors, 6)
# push!(ncolors, 10)

Random.seed!(2713)
# Create adjacency graph
A = convert(SparseMatrixCSC, Symmetric(sprand(100,100,0.1)))

const filename = joinpath(@__DIR__, "A.mtx")
MatrixMarket.mmwrite(filename, A)
verbose = false

@testset "Test ColPack" begin
    @testset "Test Matrix Market API" begin
        for (i,ordering) in enumerate(orderings)
            coloring = ColPackColoring(filename, d1_coloring(), ordering; verbose=verbose)
            @test length(unique(get_colors(coloring))) == ncolors[i]
        end
    end
    @testset "Test ADOL-C Compressed Row Storage" begin
        for (i,ordering) in enumerate(orderings)
            coloring = ColPackColoring(A, d1_coloring(), ordering; verbose=verbose)
            @test length(unique(get_colors(coloring))) == ncolors[i]
        end
    end
    @testset "Test ColPack Columns Coloring" begin
        A = [
        [1.0 1.0 0.0 0.0 0.0];
        [0.0 0.0 1.0 0.0 0.0];
        [0.0 1.0 0.0 1.0 0.0];
        [0.0 0.0 0.0 1.0 1.0];
        ]

        A = sparse(A)
        adjA = ColPack.matrix2adjmatrix(A; partition_by_rows=false)
        @test issymmetric(adjA)
        coloring = ColPackColoring(adjA, d1_coloring(), natural_ordering(); verbose=true)
        @test get_colors(coloring) == SparseDiffTools.matrix_colors(A)
    end
end


text:
# Julia version of cmu_arctic_demo.sh.

JULIABIN_PATH=get(ENV, "JULIABIN_PATH", "julia")

BIN_DIR = joinpath(Pkg.dir("VoiceConversion"), "bin")
JULIABIN = [JULIABIN_PATH, "--depwarn=no", "--color=yes"]

## Experimental conditions

src_id = "clb" # Source speaker identifier
tgt_id = "slt" # Target speaker identifier

diff = true # enable direct waveform modification based on spectrum differencial
max = 100   # maximum number of training data
order = 40  # order of mel-cepstrum (expect for 0th)
mix = 8     # number of mixtures of GMMs
power_threshold = -14.0 # power threshold to select frames used in training
niter = 30  # number of iteration in EM algorithm

skip_feature_extraction = false
skip_feature_alignment = false
skip_model_training = false
skip_voice_conversion = false

## Directory settings

cmu_arctic_root = joinpath(expanduser("~"), "data", "cmu_arctic")

feature_save_dir = "./features"
model_save_dir = "./models"
vc_save_dir_top = "./converted_wav"

if diff
    vc_save_dir = joinpath(vc_save_dir_top,
                           "$(src_id)_to_$(tgt_id)_order$(order)_gmm$(mix)_diff")
else
    vc_save_dir = joinpath(vc_save_dir_top,
                           "$(src_id)_to_$(tgt_id)_order$(order)_gmm$(mix)")
end

run(`mkdir -p $(feature_save_dir)`)
run(`mkdir -p $(model_save_dir)`)
run(`mkdir -p $(vc_save_dir)`)

## Feature Extraction

if !skip_feature_extraction
    for s in [src_id, tgt_id]
        tgt_wav_path = joinpath(cmu_arctic_root, "cmu_us_$(s)_arctic/wav")
        options = ["--max", max, "--order", order]
        cmd = `$(JULIABIN) $(BIN_DIR)/mcep.jl
        $(tgt_wav_path) $(feature_save_dir)/speakers/$(s) $options`
        run(cmd)
    end
end

## Alignment

if !skip_feature_alignment
    options = ["--max", max, "--threshold", power_threshold]
    run(`$(JULIABIN) $(BIN_DIR)/align.jl
        $(feature_save_dir)/speakers/$(src_id)
        $(feature_save_dir)/speakers/$(tgt_id)
        $(feature_save_dir)/parallel/$(src_id)_and_$(tgt_id) $options`)
end

if diff
    model_path = joinpath(model_save_dir,
                          "$(src_id)_to_$(tgt_id)_gmm$(mix)_order$(order)_diff.jld")
else
    model_path = joinpath(model_save_dir,
                          "$(src_id)_to_$(tgt_id)_gmm$(mix)_order$(order).jld")
end

if !skip_model_training
    options = ["--max", max, "--n_components", mix, "--n_iter", niter,
               "--n_init", 1]
    diff && push!(options, "--diff")
    cmd = `$(JULIABIN) $(BIN_DIR)/train_gmm.jl
    $(feature_save_dir)/parallel/$(src_id)_and_$(tgt_id)
    $(model_path) $options`
    run(cmd)
end

## Voice Conversion

if !skip_voice_conversion
    synthesis_script = diff ? ["diffvc.jl"] : ["vc.jl"]
    for n=100:105
        tgt_wav_path = joinpath(cmu_arctic_root,
                                "cmu_us_$(src_id)_arctic/wav/arctic_a0$(n).wav")
        cmd = `$(JULIABIN) $(BIN_DIR)/$(synthesis_script)
        $(tgt_wav_path)
        $model_path
        $(vc_save_dir)/arctic_a0$(n)_$(src_id)_to_$(tgt_id).wav --order $order`
        run(cmd)
    end
end


text:
### A Pluto.jl notebook ###
# v0.18.0

using Markdown
using InteractiveUtils

# ╔═╡ 456df07b-1be8-4adb-9302-9290314311e3
begin
	@enum  CartesianDirection XDirection = 1 YDirection ZDirection
	@enum  CartesianAxis XAxis = 1 YAxis ZAxis
	@enum  FieldComponent XComponent = 1 YComponent ZComponent
end

# ╔═╡ 4f7a7193-e75e-4756-9423-7d81e2f7e4db
begin

mutable struct StructuredGrid{nD} <: AbstractArray{Real,nD}
	nodes::NTuple{nD,Vector{Float64}}
end

function StructuredGrid(nodes...)
        StructuredGrid{length(nodes)}(nodes)
end

Base.size(grid::StructuredGrid)	= length.(grid.nodes)

nodes(sgrid::StructuredGrid) = sgrid.nodes
xnodes(grid::StructuredGrid)	= grid.nodes[1]
ynodes(grid::StructuredGrid)	= grid.nodes[2]
znodes(grid::StructuredGrid)	= grid.nodes[3]

xnodes(grid::StructuredGrid,i::Int)	= getindex(xnodes(grid),i)
ynodes(grid::StructuredGrid,j::Int)	=  getindex(ynodes(grid),j)
znodes(grid::StructuredGrid,k::Int)	= getindex(znodes(grid),k)

Base.getindex(grid::StructuredGrid{1},i::Int) = (xnodes(grid,i))
Base.getindex(grid::StructuredGrid{2},i::Int,j::Int) = (xnodes(grid,i),ynodes(grid,j))
Base.getindex(grid::StructuredGrid{3},i::Int,j::Int,k::Int) =  (xnodes(grid,i),ynodes(grid,j), znodes(grid,k))

cell_sizes(sgrid::StructuredGrid) = StructuredGrid(diff.(nodes(sgrid)))
cell_volumes(sgrid::StructuredGrid) = map(prod,cell_sizes(sgrid))
_avg(nodes::AbstractVector) = 0.5(nodes[2:end] + nodes[1:end-1])
cell_centroids(sgrid::StructuredGrid) = StructuredGrid(_avg.(nodes(sgrid)))

abstract type AbstractNodeType end
struct 	Primal <: AbstractNodeType  end
struct 	Dual <: AbstractNodeType  end

function primal_nodes(ax::FieldComponent, reference_grid::StructuredGrid{Dim}) where {Dim}
    StructuredGrid(ntuple(i -> i == Int(ax) ? _avg(nodes(reference_grid)[i]) : nodes(reference_grid)[i], Dim))
end
function dual_nodes(ax::FieldComponent, reference_grid::StructuredGrid{Dim}) where {Dim}
    StructuredGrid(ntuple(i -> i == Int(ax) ? nodes(reference_grid)[i] :  _avg(nodes(reference_grid)[i])  , Dim))
end

function nodes(::Type{Primal},ax::FieldComponent, reference_grid::StructuredGrid)
primal_nodes(ax,reference_grid)
end
function nodes(::Type{Dual},ax::FieldComponent, reference_grid::StructuredGrid)
dual_nodes(ax,reference_grid)
end

end

# ╔═╡ 1bc97871-d861-47d4-bde9-e2373b1e903a
begin
struct YeeGrid{T,nD}
	x::StructuredGrid{nD}
	y::StructuredGrid{nD}
	z::StructuredGrid{nD}
end

function YeeGrid{T}(refgrid::StructuredGrid{nD}) where {nD,T}
	YeeGrid{T,nD}(
		nodes(T,XComponent,refgrid),
		nodes(T,YComponent,refgrid),
		nodes(T,ZComponent,refgrid))
end
function YeeGrid{T}(nodes...) where T
	refgrid = StructuredGrid(nodes...)
	YeeGrid{T}(refgrid)
end

end

# ╔═╡ 4b3b120f-52b0-47ae-bec2-207ea12ddd6a
begin
	grid_3d_primal = YeeGrid{Primal}(LinRange(-2,2,10),LinRange(-5,5,100),LinRange(-3,3,40))
	grid_3d_dual = YeeGrid{Dual}(LinRange(-2,2,10),LinRange(-5,5,100),LinRange(-3,3,40))
	# nonuniform grid
	grid_2d_primal = YeeGrid{Primal}([1,4,7,8,9,10,15,20],LinRange(-1,1,15))
	grid_2d_dual = YeeGrid{Primal}([1,4,7,8,9,10,15,20],LinRange(-1,1,15))

	grid_1d_dual = YeeGrid{Primal}([1,4,7,8,9,10,15,20])
end

# ╔═╡ Cell order:
# ╠═456df07b-1be8-4adb-9302-9290314311e3
# ╠═4f7a7193-e75e-4756-9423-7d81e2f7e4db
# ╠═1bc97871-d861-47d4-bde9-e2373b1e903a
# ╠═4b3b120f-52b0-47ae-bec2-207ea12ddd6a


text:
"""
Struct for user supbroblem callbacks.

**Fields**

* `name::String` : User's name for the callback. Defaults to `string(h)`.
* `h::Function` : Callback function. See notes below for calling signature.
* `ext::Dict{Symbol,Any}` : Dictionary to store data between callback calls or needed parameters.

The callback function `h` must have the signature
    `h(ext::Dict{Symbol,Any}, sp::T, niter::Int, scenario_id::ScenarioID)  where T <: AbstractSubproblem`
where `ext` is the same dictionary given to the `Callback` constructor, `sp` is a concrete type of `AbstractSubproblem` (see `AbstractSubproblem`), `niter` is the current iteration and `scenario_id` is a scenario identifier (see `ScenarioID`).
"""
struct SubproblemCallback
    name::String
    h::Function
    ext::Dict{Symbol,Any}
end

function Base.show(io::IO, spcb::SubproblemCallback)
    print(io, spcb.name)
    return
end

"""
    SubproblemCallback(f::Function)

Creates a `SubproblemCallback` structure for function `f`.
"""
function SubproblemCallback(f::Function)
    return SubproblemCallback(string(f), f, Dict{Symbol,Any}())
end

"""
    SubproblemCallback(f::Function, ext::Dict{Symbol,Any})

Creates a `SubproblemCallback` structure for function `f` with the external data dictionary `ext`.
"""
function SubproblemCallback(f::Function, ext::Dict{Symbol,Any})
    return SubproblemCallback(string(f), f, ext)
end

"""
    spcb(f::Function)
    spcb(f::Function, ext::Dict{Symbol,Any})
    spcb(name::String, f::Function, ext::Dict{Symbol,Any})

Shorthand for `SubproblemCallback`.
"""
spcb(f::Function) = SubproblemCallback(f)
spcb(f::Function, ext::Dict{Symbol,Any}) = SubproblemCallback(f, ext)
spcb(name::String, f::Function, ext::Dict{Symbol,Any}) = SubproblemCallback(name, f, ext)


text:
using Distributions, UncertainData

# Test all combinations of different types of uncertain values
M = MixtureModel([Normal(3, 0.2), Normal(2, 1)])

r1 = UncertainValue(Normal, rand(), rand())
r2 = UncertainValue(rand(M, 10000))
r3 = UncertainValue(Normal, rand(Normal(4, 3.2), 10000))

uvals = [r1; r2; r3]

n = 5
for uval in uvals
    # Sine and inverse sine
    @test sin(uval) isa Vector{Float64}
    @test sin(uval, n) isa Vector{Float64}
    @test sind(uval) isa Vector{Float64}
    @test sind(uval, n) isa Vector{Float64}
    @test sinh(uval) isa Vector{Float64}
    @test sinh(uval, n) isa Vector{Float64}

    #@test asin(uval) isa Vector{Float64}
    #@test asin(uval, n) isa Vector{Float64}
    #@test asind(uval) isa Vector{Float64}
    #@test asind(uval, n) isa Vector{Float64}
    #@test asinh(uval) isa Vector{Float64}
    #@test asinh(uval, n) isa Vector{Float64}

    # Cosine and inverse cosine
    @test cos(uval) isa Vector{Float64}
    @test cos(uval, n) isa Vector{Float64}
    @test cosd(uval) isa Vector{Float64}
    @test cosd(uval, n) isa Vector{Float64}
    @test cosh(uval) isa Vector{Float64}
    @test cosh(uval, n) isa Vector{Float64}

    # @test acos(uval) isa Vector{Float64}
    # @test acos(uval, n) isa Vector{Float64}
    # @test acosd(uval) isa Vector{Float64}
    # @test acosd(uval, n) isa Vector{Float64}
    # @test acosh(uval) isa Vector{Float64}
    # @test acosh(uval, n) isa Vector{Float64}

    # Tangent and inverse tangent
    @test tan(uval) isa Vector{Float64}
    @test tan(uval, n) isa Vector{Float64}
    @test tand(uval) isa Vector{Float64}
    @test tand(uval, n) isa Vector{Float64}
    @test tanh(uval) isa Vector{Float64}
    @test tanh(uval, n) isa Vector{Float64}

    # @test atan(uval) isa Vector{Float64}
    # @test atan(uval, n) isa Vector{Float64}
    # @test atand(uval) isa Vector{Float64}
    # @test atand(uval, n) isa Vector{Float64}
    # @test atanh(uval) isa Vector{Float64}
    # @test atanh(uval, n) isa Vector{Float64}

    # Cosecant and inverse cosecant
    @test csc(uval) isa Vector{Float64}
    @test csc(uval, n) isa Vector{Float64}
    @test cscd(uval) isa Vector{Float64}
    @test cscd(uval, n) isa Vector{Float64}
    @test csch(uval) isa Vector{Float64}
    @test csch(uval, n) isa Vector{Float64}

    # @test acsc(uval) isa Vector{Float64}
    # @test acsc(uval, n) isa Vector{Float64}
    # @test acscd(uval) isa Vector{Float64}
    # @test acscd(uval, n) isa Vector{Float64}
    # @test acsch(uval) isa Vector{Float64}
    # @test acsch(uval, n) isa Vector{Float64}

    # Secant and inverse secant
    @test sec(uval) isa Vector{Float64}
    @test sec(uval, n) isa Vector{Float64}
    @test secd(uval) isa Vector{Float64}
    @test secd(uval, n) isa Vector{Float64}
    @test sech(uval) isa Vector{Float64}
    @test sech(uval, n) isa Vector{Float64}

    # @test asec(uval) isa Vector{Float64}
    # @test asec(uval, n) isa Vector{Float64}
    # @test asecd(uval) isa Vector{Float64}
    # @test asecd(uval, n) isa Vector{Float64}
    # @test asech(uval) isa Vector{Float64}
    # @test asech(uval, n) isa Vector{Float64}

    # Cotangent and inverse cotangent
    @test cot(uval) isa Vector{Float64}
    @test cot(uval, n) isa Vector{Float64}
    @test cotd(uval) isa Vector{Float64}
    @test cotd(uval, n) isa Vector{Float64}
    @test coth(uval) isa Vector{Float64}
    @test coth(uval, n) isa Vector{Float64}

    # @test acot(uval) isa Vector{Float64}
    # @test acot(uval, n) isa Vector{Float64}
    # @test acotd(uval) isa Vector{Float64}
    # @test acotd(uval, n) isa Vector{Float64}
    # @test acoth(uval) isa Vector{Float64}
    # @test acoth(uval, n) isa Vector{Float64}

    # Related trig functions
    @test sincos(uval) isa Vector{Tuple{Float64, Float64}}
    @test sincos(uval, n) isa Vector{Tuple{Float64, Float64}}

    @test sinc(uval) isa Vector{Float64}
    @test sinc(uval, n) isa Vector{Float64}

    @test sinpi(uval) isa Vector{Float64}
    @test sinpi(uval, n) isa Vector{Float64}

    @test cosc(uval) isa Vector{Float64}
    @test cosc(uval, n) isa Vector{Float64}

    @test cospi(uval) isa Vector{Float64}
    @test cospi(uval, n) isa Vector{Float64}
end


text:
#############################
##  Numerical experiments associated with section 5.2.2
#############################
## This file contains the numerical experiments used to compare our method with
## Chebfun's Volterra integral equation implementation when high polynomial orders are required.
## We only include our own side of the implementation here. See the references in the paper for Chebfun.
##

using ApproxFun, MultivariateOrthogonalPolynomials, BandedMatrices, BlockBandedMatrices, SpecialFunctions, Plots, SparseArrays
using SparseVolterraExamples

#####
## First, we can check that the stated arctan function indeed approximates a step-like function.
## This is in Figure 6(a).
v(x,K) = atan(K*x)
    plot(x->v(x,10),0,1, legend=:true , xlabel = "x", label="k=10" ,ylabel = "u_2(x,K)",  legendfontsize=12, tickfontsize=10, thickness_scaling = 1.2 , grid=:none)
    plot!(x->v(x,50),0,1, legend=:true , xlabel = "x", label="k=50" ,ylabel = "u_2(x,K)",  legendfontsize=12, tickfontsize=10, thickness_scaling = 1.2 , grid=:none)
    plot!(x->v(x,100),0,1, legend=:true , xlabel = "x", label="k=100" , ylabel = "u_2(x,k)", legendfontsize=12, tickfontsize=10, thickness_scaling = 1.2 , grid=:none)
    plot!(x->v(x,200),0,1, legend=:bottomright , xlabel = "x", label="k=200" , ylabel = "u_2(x,k)", legendfontsize=12, tickfontsize=10, thickness_scaling = 1.2 , grid=:none)

#####
## First we define the Kernel and g(x,k) as in section 5.2.2.
gf(x,k) = k/(k^2*x^2+1)-(exp(x^2)*atan(k*x))/(2*k^2)+(exp(x^2)*x)/(2*k)-1/2*exp(x^2)*x^2*atan(k*x)
Kfun(x,y) = y*exp(x^2)

#####
## Solver function with step-by-step explanation for given k and with polynomial degree 'n'
function solveSec532(k,n,gf,Kfun)
    gF = Fun(x->gf(x,k),Jacobi(1,2, 0..1),n)             # Approximate g in the appropriate basis
    V = triVolterraFullKernelOpP01(Kfun,n,true,155)      # The following steps generate the appropriate Volterra operator
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V
        V = Conversion(Jacobi(0,0,0..1),Jacobi(1,2,0..1))[1:n,1:n]*V[1:n,1:n]
        V = Derivative(Jacobi(0,1,0..1),1)[1:n,1:n]-V
    coeff = [DirectEvalLHSP10at0(n);V[1:n-1,1:n]] \ [0;pad(gF.coefficients,n-1)] # Append evaluation and initial conditions and then solve the equation
    return Fun(Jacobi(0,1, 0..1),coeff) # returns solution Fun
end

#####
## Now we can compute and plot a specific example for k=100, plotting both the numerical approximation and the analytic solution.
## Computing errors for various k is a bit more work since automatic convergence testing is not currently implemented but a straightforward for loop can do it from here.
u = solveSec532(100,50,gf,Kfun)
plot(u, label="sparse method")
    plot!(x->v(x,100),0,1, legend=:bottomright , xlabel = "x", label="analytic" ,ylabel = "u_2(x,100)",  legendfontsize=12, tickfontsize=10, thickness_scaling = 1.2 , grid=:none)

#####
## We can also use BenchmarkTools to get a somewhat robust time estimation.
## Note however that this obviously strongly depends on the hardware you are using.
using BenchmarkTools
@benchmark solveSec532(100,50,gf,Kfun)
@benchmark solveSec532(200,300,gf,Kfun)

#####
## Now we plot the operator bandedness for the k=100 example, this is basically Figure 6(b).
function OperatorSec532(k,n,Kfun)
    V = triVolterraFullKernelOpP01(Kfun,n,true,155)    # The following steps generate the appropriate Volterra operator
        V = reflectPabtoPba(n)*WLoweringP01P00(n)*V
        V = Conversion(Jacobi(0,0,0..1),Jacobi(1,2,0..1))[1:n,1:n]*V[1:n,1:n]
        V = Derivative(Jacobi(0,1,0..1),1)[1:n,1:n]-V
        V = [DirectEvalLHSP10at0(n);V[1:n-1,1:n]] # Append evaluation and initial conditions
    return V # returns operator
end
#####
## Plots.jl's spy plot is not currently compatible with these types, so we instead convert to generic sparse for visualization.
## The exact bandedness properties depend on the chosen degrees and parameters but here is a standard example
spy(sparse(OperatorSec532(100,300,Kfun)),markersize=2.8,marker=:rect)

max_stars_repo_path: examples/Section 5-2-2 - High order solutions with step-like function.jl


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

"""
LEFT INVERSE / INVOLUTION

Models are pairs of involutions inv:B -> B
AND
monomorphisms f:A->B (with right inverse g: B->A)

The main purpose of this is a simple example which shows off
diagrams.
"""

fgschema = @acset LabeledGraph_{Symbol} begin
  V = 2
  E = 3
  vlabel = [:A,:B]
  elabel = [:f,:g, :inv]
  src    = [1,  2,  2]
  tgt    = [2,  1,  2]
end

fg = FLS(:fg, fgschema,[], [
  (:fginv,      [:f, :g],    []),
  (:involution, [:inv,:inv], [])])


text:
######
##FEM Heat dt Convergence Tests
######
using FiniteElementDiffEq, DiffEqDevTools, DiffEqProblemLibrary
#Convergences estimate has not converged in this range
#Should decrease dx/dt for better estimate

cs = cs_femheat_moving_dt
alg=FEMDiffEqHeatEuler(); println(alg) #Unstable due to μ
sim = test_convergence(cs,alg)

alg=FEMDiffEqHeatImplicitEuler(); println(alg)
sim2 = test_convergence(cs,alg)

alg=FEMDiffEqHeatCrankNicholson(); println(alg) #Bound by spatial discretization error at low dt, decrease dx for full convergence
cs = cs_femheat_moving_faster_dt
sim3 = test_convergence(cs,alg)

#plot(plot(sim),plot(sim2),plot(sim3),layout=@layout([a b c]),size=(1200,400))
#Note: Stabilizes in H1 due to high dx-error, reduce dx and it converges further.

#Returns true if ImplicitEuler converges like dt and
#CN convergeces like >dt^2 (approaches dt^2 as dt and dx is smaller
minimum([abs(sim2.𝒪est[:L2]-1)<.3 abs(sim3.𝒪est[:L2]-2)<.1])


text:
import Base: @__doc__
using PowerDynamics: AbstractLine
import PowerDynamics: construct_vertex, dimension, symbolsof, showdefinition, AbstractNode
using LinearAlgebra: Diagonal, UniformScaling
using NetworkDynamics: ODEVertex

##

function total_current_real(e_s, e_d)
    # Keeping with the convention of negative sign for outging current
    real_current = 0.0
    imag_current = 0.0
    for e in e_s
        real_current -= e[1]
        imag_current -= e[2]
    end
    for e in e_d
        real_current += e[3]
        imag_current += e[4]
    end
    [real_current; imag_current]
end

##
begin
    Base.@__doc__ struct SlackReal <: AbstractNode
            U
            Y_n
        end
    SlackReal(; U, Y_n = 0) = SlackReal(U, Y_n)
    function construct_vertex(par::SlackReal)
        Ur = real(par.U)
        Ui = imag(par.U)
        Y_n = par.Y_n
        function rhs!(dx, x, e_s, e_d, p, t)
            try
                dx[1] = x[1] - Ur
                dx[2] = x[2] - Ui
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim = 2, mass_matrix = [0, 0], sym = [:u_r, :u_i])
    end
    symbolsof(::SlackReal) = begin
            [:u_r, :u_i]
        end
    dimension(::SlackReal) = begin
            2
        end
end

##

import PowerDynamics: construct_edge

begin
    @__doc__ struct LineReal <: AbstractLine
            from
            to
            Y
            LineReal(; from, to, Y) = new(from, to, Y)
        end
    function construct_edge(par::LineReal)
        from = par.from
        to = par.to
        G = real(par.Y)
        B = imag(par.Y)
        function rhs!(e, v_s, v_d, p, t)
            # source_voltage = v_s[1] + v_s[2] * im
            # destination_voltage = v_d[1] + v_d[2] * im
            # complex_current = Y * (destination_voltage - source_voltage)
            # current_vector = [complex_current, complex_current]
            e[1] = G * (v_d[1] - v_s[1]) - B * (v_d[2] - v_s[2])#real
            e[2] = G * (v_d[2] - v_s[2]) + B * (v_d[1] - v_s[1])#imag
            e[3] = G * (v_d[1] - v_s[1]) - B * (v_d[2] - v_s[2])#real
            e[4] = G * (v_d[2] - v_s[2]) + B * (v_d[1] - v_s[1])#imag
        end
        return StaticEdge(f! = rhs!, dim = 4)
    end
end

##

begin
    Base.@__doc__ struct PQReal <: AbstractNode
            P
            Q
            Y_n
        end
    PQReal(; P, Q, Y_n = 0) = PQReal(P, Q, Y_n)
    function construct_vertex(par::PQReal)
        P = par.P
        Q = par.Q
        G_n = real(par.Y_n)
        B_n = imag(par.Y_n)
        function rhs!(dx, x, e_s, e_d, p, t)
            i = total_current_real(e_s, e_d)
            i[1] += G_n * x[1]
            i[2] += B_n * x[2]
            p = x[1] * i[1] + x[2] * i[2]
            q = x[2] * i[1] - x[1] * i[2]
            #s = u * conj(i)
            #du = (P + (1im) * Q) - s
            try
                dx[1] = P - p #real(du)
                dx[2] = Q - q #imag(du)
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim = 2, mass_matrix = [0, 0], sym = [:u_r, :u_i])
    end
    symbolsof(::PQReal) = begin
            [:u_r, :u_i]
        end
    dimension(::PQReal) = begin
            2
        end
end

##

begin
    Base.@__doc__ struct ThirdOrderReal <: AbstractNode
            H
            P
            D
            Ω
            E_f
            T_d_dash
            X_q_dash
            X_d_dash
            X_d
            Y_n
        end
    ThirdOrderReal(; H, P, D, Ω, E_f, T_d_dash, X_q_dash, X_d_dash, X_d, Y_n = 0) = ThirdOrderReal(H, P, D, Ω, E_f, T_d_dash, X_q_dash, X_d_dash, X_d, Y_n)
    function construct_vertex(par::ThirdOrderReal)
        H = par.H
        P = par.P
        D = par.D
        Ω = par.Ω
        E_f = par.E_f
        T_d_dash = par.T_d_dash
        X_q_dash = par.X_q_dash
        X_d_dash = par.X_d_dash
        X_d = par.X_d
        G_n = real(par.Y_n)
        B_n = imag(par.Y_n)
        @assert H > 0 "inertia (H) should be >0"
        @assert D >= 0 "damping (D) should be >=0"
        @assert T_d_dash > 0 "time constant of d-axis (T_d_dash) should be >0"
        @assert X_q_dash >= 0
        @assert X_d_dash >= 0 "transient reactance of d-axis (X_d_dash) should be >=0"
        @assert X_d >= 0 "reactance of d-axis (X_d_dash) should be >=0"
        Ω_H = Ω / (2H)
        function rhs!(dx, x, e_s, e_d, p, t)
            i = total_current_real(e_s, e_d)
            i[1] += G_n * x[1]
            i[2] += B_n * x[2]
            θ = x[3]
            ω = x[4]
            #i_c = (1im) * i * exp((-1im) * θ)
            #e_c = (1im) * u * exp((-1im) * θ)
            #p = real(u * conj(i))
            p = x[1] * i[1] + x[2] * i[2]
            e_d = 0
            e_q = x[2] * sin(θ) + x[1] * cos(θ)
            i_d = i[1] * sin(θ) - i[2] * cos(θ)
            i_q = i[2] * sin(θ) + i[1] * cos(θ)
            dθ = ω
            dω = (((P - D * ω) - p) - (X_q_dash - X_d_dash) * i_d * i_q) * Ω_H
            de_q = (1 / T_d_dash) * ((E_f - e_q) + i_d * (X_d - X_d_dash))
            #du = de_q * exp((1im) * θ) + u * (1im) * ω
            try
                dx[1] = de_q * cos(θ) - x[2] * ω #real(du)
                dx[2] = de_q * sin(θ) + x[1] * ω #imag(du)
                dx[3] = dθ
                dx[4] = dω
                return nothing
            catch e
                if typeof(e) === UndefVarError
                    throw(NodeDynamicsError("you need to provide $(e.var)"))
                else
                    throw(e)
                end
            end
        end
        ODEVertex(f! = rhs!, dim = 4, mass_matrix = UniformScaling{Bool}(true), sym = [:u_r, :u_i, :θ, :ω])
    end
    symbolsof(::ThirdOrderReal) = begin
            [:u_r, :u_i, :θ, :ω]
        end
    dimension(::ThirdOrderReal) = begin
            4
        end
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
using Combinatorics

function parse_line(l)
    m = match(r"^(?<from>\w+) to (?<to>\w+) = (?<dist>\d+)$", l)
    (from=m[:from], to=m[:to], dist=parse(Int, m[:dist]))
end

data = readlines() .|> parse_line

function tsp(g, f)
    n = length(g)
    cities = keys(g) |> collect

    function cost(order)
        c = 0
        pos = order[1]
        for next_pos in order[2:end]
            c += g[cities[pos]][cities[next_pos]]
            pos = next_pos
        end
        c
    end

    f(cost.(permutations(1:n)))
end

function solve_1()
    g = Dict{String, Dict{String, Int}}()
    for d in data
        from, to, dist = d
        dift = get!(g, from, Dict())
        dift[to] = dist
        ditf = get!(g, to, Dict())
        ditf[from] = dist
    end

    tsp(g, minimum)
end

function solve_2()
    g = Dict{String, Dict{String, Int}}()
    for d in data
        from, to, dist = d
        dift = get!(g, from, Dict())
        dift[to] = dist
        ditf = get!(g, to, Dict())
        ditf[from] = dist
    end

    tsp(g, maximum)
end

solve_1() |> println
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
module TwoPlayerTest

using  Test

import Cribbage: CribbageGame, GameState, UnexpectedPlayerException
import Cribbage.RandomPlay: RandomPlayer
using  Cribbage.TwoPlayer

@testset "Test TwoPlayerGame Constructor" begin
    p₁ = RandomPlayer("player one")
    @test_throws AssertionError TwoPlayerGame(p₁, p₁)
    p₂ = RandomPlayer("player two")
    @test_throws AssertionError TwoPlayerGame(p₁, p₂, 42)
    @test isa(TwoPlayerGame(p₁, p₂), CribbageGame)
    @test isa(TwoPlayerGame(p₁, p₂, 61), CribbageGame)
    @test isa(TwoPlayerGame(p₁, p₂, 121), CribbageGame)

end

examplegame = TwoPlayerGame(
    RandomPlayer("1"),
    RandomPlayer("2"),
    121
)
badplayer = RandomPlayer("-1")

@testset "Test TwoPlayerGameState Constructor" begin
    @test_throws AssertionError TwoPlayerGameState(examplegame, badplayer)
    @test_throws AssertionError TwoPlayerGameState(examplegame, examplegame.p₁, -42)
    @test_throws AssertionError TwoPlayerGameState(examplegame, examplegame.p₁, 0, -53)
    @test isa(TwoPlayerGameState(examplegame, examplegame.p₁, 0, 0), GameState)
end

@testset "Test players" begin
    @test length(players(examplegame)) == 2
    for player in players(examplegame)
        @test player ≡ examplegame.p₁ || player ≡ examplegame.p₂
    end
end

@testset "Test incrementscore" begin
    initialstate = TwoPlayerGameState(
        examplegame, examplegame.p₁, 0, 0
    )
    @test_throws UnexpectedPlayerException incrementscore(initialstate, badplayer, 0)
    @test_throws AssertionError incrementscore(initialstate, examplegame.p₁, -5)

    @test incrementscore(initialstate, examplegame.p₁, 0) == initialstate

    newstate = incrementscore(initialstate, examplegame.p₁, 6)
    @test newstate.dealer ≡ initialstate.dealer
    @test newstate.s₁ == 6
    @test newstate.s₂ == 0

    initialstate = TwoPlayerGameState(
        examplegame, examplegame.p₁, 100, 115
    )
    newstate = incrementscore(initialstate, examplegame.p₂, 8)
    @test newstate.dealer ≡ initialstate.dealer
    @test newstate.s₁ == 100
    @test newstate.s₂ == 121

end

@testset "Test isgameover" begin

    gs = TwoPlayerGameState(examplegame, examplegame.p₁, 112, 118)
    @test !isgameover(gs)
    gs = TwoPlayerGameState(examplegame, examplegame.p₁, 121, 118)
    @test isgameover(gs)
    gs = TwoPlayerGameState(examplegame, examplegame.p₁, 100, 121)
    @test isgameover(gs)

end

@testset "Test reportscore" begin

    gs = TwoPlayerGameState(examplegame, examplegame.p₁, 121, 118)
    @test isa(reportscore(gs), String)

end

@testset "Test rotatedealer" begin
    initialstate = TwoPlayerGameState(examplegame, examplegame.p₁, 1, 2)
    newstate = rotatedealer(initialstate)
    @test newstate.dealer ≡ examplegame.p₂
    @test newstate.s₁ == initialstate.s₁
    @test newstate.s₂ == initialstate.s₂
    newstate = rotatedealer(newstate)
    @test newstate == initialstate
end

end  # module TwoPlayerTest


text:
module TestSeason
  using DataFrames
  using Base.Test
  include("../src/season.jl")

  # Construction
  s10 = Season(2010)
  s10s = Season("2010")
  s10c = Season(s10.off, s10.def, s10.k)
  s11 = Season(2011)
  s12 = Season(2012)
  s13 = Season(2013)
  s14 = Season(2014)

  @test typeof(s10) == Season
  @test typeof(s10s) == Season
  @test typeof(s10.off.qb) == DataFrame
  @test typeof(s10s.off.qb) == DataFrame
  @test typeof(s10.off.rb) == DataFrame
  @test typeof(s10.off.te) == DataFrame
  @test typeof(s10.off.wr) == DataFrame
  @test typeof(s10.def.def) == DataFrame
  @test typeof(s10.def.wr) == DataFrame
  @test typeof(s10.k.k) == DataFrame
  @test s10c.off == s10.off
  @test s10c.def == s10.def
  @test s10c.k == s10.k

  # filter
  nwe14 = filter(s14, "NWE")
  @test nwe14.k.k[:Name] == ["Stephen Gostkowski"]
  @test nwe14.def.def[:team_abbr] == ["NWE"]
  @test nwe14.def.wr[:team_abbr] == ["NWE", "NWE", "NWE", "NWE", "NWE", "NWE"]
  @test nwe14.off.qb[:Name] == ["Jimmy Garoppolo", "Tom Brady"]
  @test nwe14.off.te[:team_abbr] == ["NWE", "NWE", "NWE", "NWE"]
  @test nwe14.off.wr[:team_abbr] == ["NWE", "NWE", "NWE", "NWE", "NWE", "NWE"]
  @test nwe14.off.rb[:team_abbr] == ["NWE", "NWE", "NWE", "NWE", "NWE", "NWE", "NWE"]

  # filter - specific dispatch
  nwe14k = filter(s14.k, "NWE")
  @test nwe14k.k[:Name] == ["Stephen Gostkowski"]
  nwe14d = filter(s14.def, "NWE")
  @test nwe14d.def[:team_abbr] == ["NWE"]
  @test nwe14d.wr[:Name] == [
    "Julian Edelman",
    "Brandon LaFell",
    "Danny Amendola",
    "Brian Tyms",
    "Aaron Dobson",
    "Matthew Slater"
  ]
  nwe14o = filter(s14.off, "NWE")
  @test nwe14o.qb[:Name] == [
    "Jimmy Garoppolo",
    "Tom Brady"
  ]
  @test nwe14o.wr[:Name] == [
    "Julian Edelman",
    "Brandon LaFell",
    "Danny Amendola",
    "Brian Tyms",
    "Aaron Dobson",
    "Matthew Slater"
  ]
  @test nwe14o.rb[:Name] == [
    "LeGarrette Blount",
    "Jonas Gray",
    "Shane Vereen",
    "Stevan Ridley",
    "Brandon Bolden",
    "James White",
    "James Develin"
  ]
  @test nwe14o.te[:Name] == [
    "Rob Gronkowski",
    "Tim Wright",
    "Michael Hoomanawanui",
    "Steve Maneri"
  ]

  # score
  nwe14s = score(nwe14)
  @test sort(nwe14s.off.qb[:score]) == [11, 277]
  @test sort(nwe14s.off.wr[:score]) == [0, 3, 14, 26, 121, 135]
  @test sort(nwe14s.off.te[:score]) == [0, 4, 49, 160]
  @test sort(nwe14s.off.rb[:score]) == [4, 5, 14, 48, 71, 87, 113]
  @test sort(nwe14s.k.k[:score]) == [49]
  @test sort(nwe14s.def.def[:score]) == [138]
  @test sort(nwe14s.def.def[:score_def]) == [132]
  @test sort(nwe14s.def.wr[:score]) == [0, 0, 0, 0, 0, 6]

  # Score - specific dispatch
  nwe14os = score(nwe14.off)
  @test sort(nwe14os.qb[:score]) == [11, 277]
  @test sort(nwe14os.wr[:score]) == [0, 3, 14, 26, 121, 135]
  @test sort(nwe14os.te[:score]) == [0, 4, 49, 160]
  @test sort(nwe14os.rb[:score]) == [4, 5, 14, 48, 71, 87, 113]
  nwe14ks = score(nwe14.k)
  @test sort(nwe14ks.k[:score]) == [49]
  nwe14ds = score(nwe14.def)
  @test sort(nwe14ds.def[:score]) == [138]
  @test sort(nwe14ds.def[:score_def]) == [132]
  @test sort(nwe14ds.wr[:score]) == [0, 0, 0, 0, 0, 6]

  # Source data resolver helpers
  # FIXME *nix specific paths
  @test rsearch(table_path("2014", "DEF"), "/src/../data/2014/DEF.csv").start > 0
  @test rsearch(table_logical_path("2014", "DEF"), "/src/../data/2014/DEF").start > 0
end

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


text:
load_unit_system("default")
#load_unit_system("cgs")

# to use cgs instead of default, just comment out line 1 and uncomment line 2
# other systems can be added easily, just make a new folder, fill it with
# .jl files that define your prefered system and constants, then add
# load_unit_system("FOLDER_NAME") to this file and comment out everything else


text:
@testset "LinearOperator" begin
    coeffs = rand(3*3, 2+3)
    A = LinearOperator(Taylor(1) × Chebyshev(2), Fourier(1, 1.0)^3, coeffs)
    @test domain(A) == Taylor(1) × Chebyshev(2)
    @test codomain(A) == Fourier(1, 1.0)^3
    @test coefficients(A) == coeffs
    @test firstindex(A, 1) == 1
    @test firstindex(A, 2) == 1
    @test lastindex(A, 1) == 9
    @test lastindex(A, 2) == 5
    @test size(A) == (9, 5)
    @test length(A) == 45
end


text:
@coeff function c(reactor::Reactor)
  this_term = log((1+epsilon_B)/(1-epsilon_B))/(1-epsilon_B)

  that_term = 4*epsilon_B/(1-epsilon_B^2)

  cur_value = (
    ( this_term + that_term ) *
    R_0 * B_0 ^ 2 / ( 4 * mu_0 * sigma_max )
  )

  cur_value
end


text:
module Poirot

using Reexport, Printf, IRTools.All, QuadGK
@reexport using Distributions
@reexport using Statistics
using MacroTools: @forward
using Mjolnir: @trace

export Rejection
export infer, observe

include("compiler/trace.jl")
include("compiler/simplify.jl")
include("compiler/logprob.jl")

include("inference/distributions.jl")
include("inference/inference.jl")
include("inference/trivial.jl")
include("inference/quad.jl")
include("inference/rejection.jl")

end # module


text:
####################
## construct vine ##
####################

## from second layer on
function vineFrom2(vn::Vine)
    ## build vine from given first layer
    ##
    ## Inputs:
    ## 	vn			Vine with only first layer specified
    ## Outputs:
    ## 	finalVn	completely specified Vine

    nVars = size(vn.trees, 2)
    nLayers = nVars-1

    for layerNb=2:nLayers
        ## get current vine in tree notation
        trs = convert(CTreePaths, vn)

        ## get new conditional edges
        newEdges = selectLinks(trs, layerNb)

        ## add new edges to vine
        for e in newEdges
            addEdge!(vn, e)
        end
    end
    return vn
end

function addEdge!(vn::Vine, edge::CondEdge)
    v1, v2 = edge.vertices
    p1, p2 = edge.parents
    vn.trees[v1, v2] = p2
    vn.trees[v2, v1] = p1
    return vn
end

##################
## select edges ##
##################

## select edges out of list of possible new edges
function selectLinks(trs::Array{CTreePaths, 1}, nrLayer::Int)
    ## some edges are automatically selected, because they need to be
    ## picked
    ##
    ## some edges may only be selected together
    ##
    ## some edges are selected manually

    nVars = length(trs)
    nToMatch = nVars - nrLayer          # number of required edges

    ## find incomplete trees
    verticesToBeLinked = incomplTrees(trs)

    ## get possible matches
    ##---------------------
    possibleEdges = findPairMatch(trs, nrLayer, verticesToBeLinked)
    ## displayAllPossibleEdges(possibleEdges)

    nMatches = length(possibleEdges)

    ## validity checks
    if nMatches < nToMatch
        error("Too few matching conditional sets!")
    end

    reachedVs = reachedVertices(possibleEdges)
    if !(reachedVs == verticesToBeLinked)
        error("Not each required vertex can be connected!")
    end

    ## automatic selection: single possible edge
    ##------------------------------------------
    chosenEdges, stillToBeLinked = mandatoryEdges(possibleEdges,
                                                  verticesToBeLinked,
                                                  nVars)

    ## match remaining links
    ##----------------------

    nRemain = nToMatch - sum(chosenEdges)

    if nRemain > 0
        displayCurrentLayer(nrLayer)
        println("Current tree:")
        display(trs)
        displayAutomChosenEdges(possibleEdges, chosenEdges)
    end

    while nRemain > 0
        displayLayerState(possibleEdges, chosenEdges,
                          stillToBeLinked, nRemain)
        selectedInd = getUserInput("Which edge(s) should be chosen?")

        indConversion = find(!chosenEdges)
        indInList = indConversion[selectedInd]

        chosenEdges[indInList] = true

        ## rm vertices of chosen edge from stillToBeLinked list
        for eInd in indInList
            strikeOffVs = possibleEdges[eInd].vertices
            stillToBeLinked = setdiff(stillToBeLinked, strikeOffVs)
        end

        ## new remaining
        nRemain = nToMatch - sum(chosenEdges)
    end

    ## backtest: are all required edges connected through set of
    ## chosen edges?
    chosenEdgesList = possibleEdges[chosenEdges]
    return chosenEdgesList
end

function setFirstLayer(vertices::Array{Int, 2})
    nEdges = size(vertices, 1)
    nVars = nEdges + 1
    parNot = -1*ones(Int, nVars, nVars)
    for ii=1:nEdges
        v1, v2 = vertices[ii, :]
        parNot[v2, v1] = v1
        parNot[v1, v2] = v2
    end
    for ii=1:nVars
        parNot[ii, ii] = 0
    end
    return Vine(parNot)
end

## function getUserInput(prompt::String="")
##     print(prompt)
##     return chomp(readline())
## end
function getUserInput(prompt::String="")
    print(prompt)
    inputStr = chomp(readline())
    individualNumbers = split(inputStr, " ")
    inds = Int[int(nbStr) for nbStr in individualNumbers]
    return inds
end

function displayLayerState(possibleEdges::Array{CondEdge, 1},
                           chosenEdges::Array{Bool, 1},
                           stillToBeLinked::Array{Int, 1},
                           nRemain::Int)

    println("Variables that still need to be connected:")
    display(stillToBeLinked)
    println("-------------")

    println("$nRemain out of the following edges must be chosen:")
    remainingEdges = possibleEdges[!chosenEdges]
    for ii=1:length(remainingEdges)
        e = remainingEdges[ii]
        condSetString = "$(e.condSet[1])"
        for jj=2:length(e.condSet)
            condSetString =
                join([condSetString, ",", e.condSet[jj]])
        end
        edgeString = "$(e.vertices[1]),$(e.vertices[2]) | $condSetString"
        println("$ii |   $edgeString")
    end
    ## display(remainingEdges)
    println("-------------")
end

function displayCurrentLayer(nrLayer::Int)
    println("\n-------------")
    println("Manually select links at layer $nrLayer:")
    println("-------------")
end

function displayAllPossibleEdges(possibleEdges::Array{CondEdge, 1})
    println("Set of all possible new edges:")
    display(possibleEdges)
    println("-------------")
end

function displayAutomChosenEdges(possibleEdges::Array{CondEdge, 1},
                                 chosenEdges::Array{Bool, 1})
    println("Automatically chosen edges are:")
    display(possibleEdges[chosenEdges])
    println("-------------")
end

function mandatoryEdges(possibleEdges::Array{CondEdge, 1},
                        verticesToBeLinked::Array{Int, 1},
                        nVars::Int)

    nMatches = length(possibleEdges)
    chosenEdges = zeros(Bool, nMatches)
    stillToBeLinked = verticesToBeLinked

    ## calculate vertex table for given edge list
    vTable = vertexTable(possibleEdges, nVars)

    ## find vertices that need to be matched and only have one
    ## possible edge
    lookForUniqueEdge = true
    while lookForUniqueEdge
        eInd = findUniqueEdgeForVertex(vTable, stillToBeLinked)
        if isempty(eInd)
            lookForUniqueEdge = false
        else
            ## mark edge as already selected
            chosenEdges[eInd] = true

            ## remove edge vertices from toBeLinked list
            strikeOffVs = possibleEdges[eInd].vertices
            stillToBeLinked = setdiff(stillToBeLinked, strikeOffVs)
        end
    end
    return (chosenEdges, stillToBeLinked)
end

## tests exist
function incomplTrees(trs::Array{CTreePaths, 1})
    ## only incomplete trees need to be considered
    ## Input:
    ## 	trs			nx1 Array of Trees
    ## Out:
    ## 	incompl		Array{Int, 1} of indices of incomplete trees,
    ## 					chronologically sorted

    nTrees = length(trs)
    incompl = Array(Int, 0)
    for ii=1:nTrees
        if !(isempty(unlinkedVars(trs[ii], nTrees)))
            push!(incompl, ii)
        end
    end
    return incompl
end

## tests exist
function equSets(arr1::Array{Int, 1}, arr2::Array{Int, 1})
    ## return unique(sort(arr1)) == unique(sort(arr2))
    return Set(arr1) == Set(arr2)
end

function findUniqueEdgeForVertex(vTable::Array{Array{Int, 1}, 1},
                                toBeMatched::Array{Int, 1})
    edgeNumber = []
    for var in toBeMatched
        if length(vTable[var]) == 1
            edgeNumber = vTable[var][1]
            break
        end
    end
    return edgeNumber
end

#########################
## find possible edges ##
#########################

## for two given trees, find matches
function pairMatch(tr1::CTreeParRef, tr2::CTreeParRef, nrLayer::Int)
    ## check for correct layer?
    ## check whether already matched?
    ## nrLayer > 1
    ## Note:
    ## allows for multiple matching condition sets, but in reality
    ## this probably may not occur!

    if !(nrLayer > 1)
        error("First layer needs to be built elsewhere!")
    end

    condSetLen = nrLayer - 1
    validLen1 = find(getDepths(tr1) .== condSetLen)
    validLen2 = find(getDepths(tr2) .== condSetLen)

    matchingInd1 = Array(Int, 0)
    matchingInd2 = Array(Int, 0)
    for ind1 in validLen1
        for ind2 in validLen2
            if equSets(tr1[ind1], tr2[ind2])
                push!(matchingInd1, ind1)
                push!(matchingInd2, ind2)
            end
        end
    end
    return [matchingInd1 matchingInd2]
end

## for given layer, find all matches in tree array
## tests exist
function findPairMatch(trs::Array{CTreePaths, 1}, nrLayer::Int,
                       verticesToBeLinked::Array{Int, 1})
    ## constraints:
    ## - individual trees are sorted: tree at index ii has root ii
    ## - if ii is unlinked with jj, so is jj with ii
    ## - nVars == nTrees

    nTrees = length(trs)

    ## preallocate list of edges
    matchingEdges = Array(CondEdge, 0)

    for ii=1:(nTrees-1)
        for jj=(ii+1):nTrees
            ## exclude already complete trees
            isRequired1 = ii in verticesToBeLinked
            isRequired2 = jj in verticesToBeLinked

            if isRequired1 && isRequired2
                matches = pairMatch(trs[ii], trs[jj], nrLayer)
                if !isempty(matches)
                    e = CondEdge(ii, jj,
                                 trs[ii][matches[1]],
                                 trs[jj][matches[2]])
                    push!(matchingEdges, e)
                end
            end
        end
    end
    return matchingEdges
end

function unlinkedVars(tr::AbstractCTree, nVars::Int)
    ## only variables not yet linked need to be tested
    return setdiff([1:nVars], allNodes(tr))
end

function reachedVertices(edges::Array{CondEdge, 1})
    ## get chronological list of all connected vertices
    reachedVs = Array(Int, 0)
    nEdges = length(edges)
    for ii=1:nEdges
        v1, v2 = edges[ii].vertices
        push!(reachedVs, v1)
        push!(reachedVs, v2)
    end
    return sort(unique(reachedVs))
end


text:
function [Grid] = PolygonGrid(np,OrientFace)
dphi=2*pi/np;

phi=0;
NumNodes=np;
Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
for i=1:np
  x=sin(phi);
  y=-cos(phi);
  z=0;
  Nodes(NodeNumber)=Node([x;y;z],NodeNumber);
  NodeNumber=NodeNumber+1;
  phi=phi+dphi;
end
Grid.Nodes=Nodes;
Grid.NumNodes=NumNodes;

NumEdges=np;
Edges(1:NumEdges)=Edge([1 2],Grid,0,0);
EdgeNumber=1;
N1=1;
N2=2;
for i=1:np-1
  Edges(EdgeNumber)=Edge([N1 N2],Grid,EdgeNumber,EdgeNumber);
  EdgeNumber=EdgeNumber+1;
  N1=N2;
  N2=N2+1;
end
Edges(EdgeNumber)=Edge([N1 1],Grid,EdgeNumber,EdgeNumber);
Grid.Edges=Edges;
Grid.NumEdges=NumEdges;
Grid.NumEdgesI=0;
Grid.NumEdgesB=np;

NumFaces=1;
Faces(1:NumFaces)=Face([0 0 0 0],Grid,0,'x');
FaceNumber=1;
Type='o';
Faces(1)=Face(1:np,Grid,FaceNumber,Type,OrientFace);

Grid.Faces=Faces;
Grid.NumFaces=NumFaces;

end


text:
function calculate_transcriptional_control_array(t::Float64,x::Array{Float64,1}, problem::Dict{String,Any})::Float64

    # initialize -
    u_variable = 0.0

    # alias elements of the species vector -
    mRNA = x[1]
    G = x[2]
    σ70 = x[3]

    # get stuff from the problem dictionary -
    E1 = problem["E1"] #kJ/mol
    E2 = problem["E2"] #kJ/mol
    R = problem["ideal_gas_constant_R"] #kJ/mol-K
    T_K = problem["temperature_K"] #K
    KD = problem["inducer_dissociation_constant"] #nM
    n = problem["inducer_cooperativity_parameter"] #dimensionless
    #k = problem["boltzmann_factor"] #kJ/K

    # TODO: write u-varible function here
    #u_variable = (W1 + W2*fi)/(1 + W1 + W2*fi)
    #W1 = exp(-E1/(k*T_K))
    #W2 = exp(-E2/(k*T_K))
    #fi = (I^n) / ((KD^n) + (I^n))
    u_variable = (exp(-E1/(R*T_K)) + exp(-E2/(R*T_K)) * ((σ70^n) / ((KD^n) + (σ70^n)))) / (1 + exp(-E1/(R*T_K)) + exp(-E2/(R*T_K)) * ((σ70^n) / ((KD^n) + (σ70^n))))

    # return -
    return u_variable
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/LICENSE

"""
    calculate_bandwidth(G::Dict{Int, Vector{Int}})

Calculate bandwidth of graph ``G``. The bandwidth of the matrix is number ``k``
such that ``A_{ij} = 0`` if ``|i-j| > k``.

# Examples

```jldoctest
julia> G = Dict(
           1 => [3, 8, 9],
           2 => [3, 8, 7],
           3 => [1, 2],
           4 => [8, 9],
           5 => [7, 8],
           6 => [2, 7],
           7 => [5, 2, 6],
           8 => [1, 2, 4, 5],
           9 => [1, 4]);

julia> calculate_bandwidth(G)
17

```

# References

* Wikipedia contributors. "Band matrix." Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 28 Apr. 2017. Web. 15 Jul. 2017. https://en.wikipedia.org/wiki/Band_matrix.
"""
function calculate_bandwidth(G::Dict{Int, Vector{Int}})
    return 2*maximum([maximum(j-g) for (j,g) in G]) + 1
end