text:
struct Wilson_Dirac_1storder_operator{Dim,T,fermion} <: Dirac_operator{Dim} where T <: AbstractGaugefields
    U::Array{T,1}
    μ::Int64
    boundarycondition::Vector{Int8}
    _temporary_fermi::Vector{fermion}
end

function Wilson_Dirac_1storder_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},μ,x,boundarycondition) where {NC,Dim}
    xtype = typeof(x)
    num = 3
    _temporary_fermi = Array{xtype,1}(undef,num)
    for i=1:num
        _temporary_fermi[i] = similar(x)
    end
    T = eltype(U)

    return Wilson_Dirac_1storder_operator{Dim,T,xtype}(U,μ,boundarycondition,_temporary_fermi )
end

struct Wilson_Dirac_operator_faster{Dim,T,fermion} <: Wilson_Dirac_operators{Dim}  where T <: AbstractGaugefields
    U::Array{T,1}
    D::Vector{Wilson_Dirac_1storder_operator{Dim,T,fermion}}
    κ::Float64 #Hopping parameter
    _temporary_fermi::Vector{fermion}
    factor::Float64
    boundarycondition::Vector{Int8}
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    verbose_print::Verbose_print
    _temporary_fermion_forCG::Vector{fermion}

    γ::Array{ComplexF64,3}
    rplusγ::Array{ComplexF64,3}
    rminusγ::Array{ComplexF64,3}
    r::Float64 #Wilson term
    hopp::Array{ComplexF64,1}
    hopm::Array{ComplexF64,1}

end

#=
function Wilson_Dirac_operator_faster(U::Array{<: AbstractGaugefields{NC,Dim},1},x,κ,boundarycondition;factor=1) where {NC,Dim}
    T = eltype(U)
    xtype = typeof(x)
    D = Vector{Wilson_Dirac_1storder_operator{Dim,T,xtype}}(undef,Dim)
    for μ=1:Dim
        D[μ] = Wilson_Dirac_1storder_operator(U,μ,x,boundarycondition)
    end
    num = 5
    _temporary_fermi = Array{xtype,1}(undef,num)
    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    numcg = 7
    #numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    eps_CG = default_eps_CG
    #println("eps_CG = ",eps_CG)
    MaxCGstep = default_MaxCGstep

    verbose_level = 2
    verbose_print = Verbose_print(verbose_level)

    method_CG = "bicg"

    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    return Wilson_Dirac_operator_faster{Dim,T,xtype}(U,D,κ,_temporary_fermi,factor,
        boundarycondition,
        eps_CG,MaxCGstep,verbose_level,
        method_CG,
        verbose_print,
        _temporary_fermion_forCG)
end
=#

function Wilson_Dirac_operator_faster(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where {NC,Dim}
    T = eltype(U)
    xtype = typeof(x)
    @assert haskey(parameters,"κ") "parameters should have the keyword κ"

    κ = parameters["κ"]

    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end

    D = Vector{Wilson_Dirac_1storder_operator{Dim,T,xtype}}(undef,Dim)
    for μ=1:Dim
        D[μ] = Wilson_Dirac_1storder_operator(U,μ,x,boundarycondition)
    end

    num = 5
    _temporary_fermi = Array{xtype,1}(undef,num)
    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    factor = check_parameters(parameters,"factor",1)#1/2κ)

    numcg = check_parameters(parameters,"numtempvec_CG",7)
    #numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    #println("eps_CG = ",eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)
    verbose_print = Verbose_print(verbose_level)

    method_CG = check_parameters(parameters,"method_CG","bicg")

    r = check_parameters(parameters,"r",1.0)

    if Dim==4
        γ,rplusγ,rminusγ = mk_gamma(r)
        hopp = zeros(ComplexF64,4)
        hopm = zeros(ComplexF64,4)
        hopp .= κ
        hopm .= κ
    elseif Dim == 2
        γ,rplusγ,rminusγ = mk_sigma(r)
        hopp = zeros(ComplexF64,2)
        hopm = zeros(ComplexF64,2)
        hopp .= κ
        hopm .= κ
    end

    return Wilson_Dirac_operator_faster{Dim,T,xtype}(U,D,κ,_temporary_fermi,factor,
        boundarycondition,
        eps_CG,MaxCGstep,verbose_level,
        method_CG,
        verbose_print,
        _temporary_fermion_forCG,
        γ,#::Array{ComplexF64,3}
        rplusγ,#::Array{ComplexF64,3}
        rminusγ,#::Array{ComplexF64,3}
        r,#::Float64 #Wilson term
        hopp,#::Array{ComplexF64,1}
        hopm#::Array{ComplexF64,1})
    )
end

#=
struct DdagD_Wilson_operator_faster{Dim,T,fermion} <: DdagD_operator
    dirac::Wilson_Dirac_operator_faster{Dim,T,fermion}
    function DdagD_Wilson_operator_faster(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
        return new{Dim,eltype(U),typeof(x)}(Wilson_Dirac_operator_faster(U,x,parameters))
    end

    function DdagD_Wilson_operator_faster(D::Wilson_Dirac_operator{Dim,T,fermion}) where {Dim,T,fermion}
        return new{Dim,T,fermion}(D)
    end
end
=#

function (D::Wilson_Dirac_1storder_operator{Dim,T,fermion})(U) where {Dim,T,fermion}
    return Wilson_Dirac_1storder_operator{Dim,T,fermion}(
        U,
        D.μ,
        D.boundarycondition,
        D._temporary_fermi)
end

function (D::Vector{Wilson_Dirac_1storder_operator{Dim,T,fermion}})(U) where {Dim,T,fermion}
    for μ=1:Dim
        D[μ] = D[μ](U)
    end
    return D
end

function (D::Wilson_Dirac_operator_faster{Dim,T,fermion})(U) where {Dim,T,fermion}
    return Wilson_Dirac_operator_faster{Dim,T,fermion}(
        U,
        D.D(U),
        D.κ,
        D._temporary_fermi,
        D.factor,
        D.boundarycondition,
        D.eps_CG,D.MaxCGstep,D.verbose_level,
        D.method_CG,
        D.verbose_print,
        D._temporary_fermion_forCG,
        D.γ,#::Array{ComplexF64,3}
        D.rplusγ,#::Array{ComplexF64,3}
        D.rminusγ,#::Array{ComplexF64,3}
        D.r,#::Float64 #Wilson term
        D.hopp,#::Array{ComplexF64,1}
        D.hopm,#::Array{ComplexF64,1})
        )
end

struct Adjoint_Wilson_operator_faster{T} <: Adjoint_Dirac_operator
    parent::T
end

function Base.adjoint(A::T) where T <: Wilson_Dirac_operator_faster
    Adjoint_Wilson_operator_faster{typeof(A)}(A)
end

struct Adjoint_Wilson_Dirac_1storder_operator{T} <: Adjoint_Dirac_operator
    parent::T
end

function Base.adjoint(A::T) where T <: Wilson_Dirac_1storder_operator
    Adjoint_Wilson_Dirac_1storder_operator{typeof(A)}(A)
end

"""
ψ_n - κ sum_ν U_n[ν](1 - γν)*ψ_{n+ν} + U_{n-ν}[-ν]^+ (1 + γν)*ψ_{n-ν}
"""
function LinearAlgebra.mul!(y::T1,A::Wilson_Dirac_operator_faster{Dim,T,fermion},x::T3) where {T1 <:AbstractFermionfields,T, Dim,fermion, T3 <:AbstractFermionfields}
    clear_fermion!(y)
    add_fermion!(y,A.factor,x)
    for μ=1:Dim
        mul!(A._temporary_fermi[1],A.D[μ],x)
        add_fermion!(y,-A.factor*A.κ,A._temporary_fermi[1])
    end
end

function LinearAlgebra.mul!(y::T1,A::Adjoint_Wilson_operator_faster{Wilson_Dirac_operator_faster{Dim,T,fermion}},x::T3) where {T1 <:AbstractFermionfields,T, Dim,fermion, T3 <:AbstractFermionfields}
    clear_fermion!(y)
    add_fermion!(y,A.parent.factor,x)
    for μ=1:Dim
        mul!(A.parent._temporary_fermi[1],A.parent.D[μ]',x)
        add_fermion!(y,-A.parent.factor*A.parent.κ,A.parent._temporary_fermi[1])
    end
end

function LinearAlgebra.mul!(y::T1,A::Wilson_Dirac_1storder_operator{Dim,T,fermion},x::T3) where {T1 <:AbstractFermionfields,T, Dim,fermion, T3 <:AbstractFermionfields}
    if A.μ == 1
        apply_Dirac_1storder_1!(y,x,A.U,A.boundarycondition,A._temporary_fermi)
    elseif A.μ == 2
        apply_Dirac_1storder_2!(y,x,A.U,A.boundarycondition,A._temporary_fermi)
    elseif A.μ == 3
        apply_Dirac_1storder_3!(y,x,A.U,A.boundarycondition,A._temporary_fermi)
    elseif A.μ == 4
        apply_Dirac_1storder_4!(y,x,A.U,A.boundarycondition,A._temporary_fermi)
    else
        error("μ = $(A.μ) is not supported!!")
    end
end

function LinearAlgebra.mul!(y::T1,A::Adjoint_Wilson_Dirac_1storder_operator{Wilson_Dirac_1storder_operator{Dim,T,fermion}},x::T3) where {T1 <:AbstractFermionfields,T, Dim,fermion, T3 <:AbstractFermionfields}
    if A.parent.μ == 1
        apply_Dirac_1storder_1_dagger!(y,x,A.parent.U,A.parent.boundarycondition,A.parent._temporary_fermi)

    elseif A.parent.μ == 2
        apply_Dirac_1storder_2_dagger!(y,x,A.parent.U,A.parent.boundarycondition,A.parent._temporary_fermi)
    elseif A.parent.μ == 3
        apply_Dirac_1storder_3_dagger!(y,x,A.parent.U,A.parent.boundarycondition,A.parent._temporary_fermi)
    elseif A.parent.μ == 4
        apply_Dirac_1storder_4_dagger!(y,x,A.parent.U,A.parent.boundarycondition,A.parent._temporary_fermi)
    else
        error("μ = $(A.parent.μ) is not supported!!")
    end
end

const γ,rplusγ,rminusγ = mk_gamma(1)

"""
U_n[ν](1 - γν)*ψ_{n+ν} + U_{n-ν}[-ν]^+ (1 + γν)*ψ_{n-ν}
"""
function apply_Dirac_1storder_1!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 1
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1minusγ1x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1plusγ1x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_2!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 2
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1minusγ2x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1plusγ2x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_3!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 3
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1minusγ3x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1plusγ3x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_4!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 4
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1minusγ4x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1plusγ4x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_1_dagger!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 1
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1plusγ1x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1minusγ1x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_2_dagger!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 2
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1plusγ2x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1minusγ2x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_3_dagger!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 3
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1plusγ3x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1minusγ3x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function apply_Dirac_1storder_4_dagger!(y,x,U,boundarycondition,_temporary_fermi)

    ν = 4
    temp1 = _temporary_fermi[1]
    Udagx= _temporary_fermi[2]
    #Ux = _temporary_fermi[3]
    clear_fermion!(y)
    mul!(Udagx,U[ν]',x)
    set_wing_fermion!(Udagx)

    xplus = shift_fermion(x,ν)
    mul!(temp1,U[ν],xplus)
    mul_1plusγ4x!(y,temp1)
    #mul!(y,view(rminusγ,:,:,ν),temp1)

    xminus = shift_fermion(Udagx,-ν)
    #xminus = shift_fermion(x,-ν)
    #Uminus = shift_U(U[ν],-ν)
    #mul!(temp1,Uminus',xminus)

    mul_1minusγ4x!(temp1,xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),xminus)
    #mul!(temp1,view(rplusγ,:,:,ν),temp1)
    add_fermion!(y,1,temp1)
    set_wing_fermion!(y,boundarycondition)

    #display(xout)
    #    exit()
    return
end

function D4x!(xout::T1,U::Array{G,1},x::T2,A::T,Dim) where  {T1,T2,
                G <: AbstractGaugefields,T <: Wilson_Dirac_operator_faster}
    clear_fermion!(xout)
    for μ=1:Dim
        mul!(A._temporary_fermi[1],A.D[μ],x)
        add_fermion!(xout,0.5,A._temporary_fermi[1])
    end
    set_wing_fermion!(xout)
end

function D4dagx!(xout::T1,U::Array{G,1},x::T2,A::T,Dim) where  {T1,T2,
    G <: AbstractGaugefields,T <: Wilson_Dirac_operator_faster}
    clear_fermion!(xout)
    for μ=1:Dim
        mul!(A._temporary_fermi[1],A.D[μ]',x)
        add_fermion!(xout,0.5,A._temporary_fermi[1])
    end
    set_wing_fermion!(xout)
end


text:
using Comonicon
using Test
using Pkg
using SimpleMock
using PackageCompiler
using Comonicon.PATH
using Comonicon.BuildTools
using Comonicon.BuildTools: write_path, contain_comonicon_path, contain_comonicon_fpath, read_toml

Pkg.activate(PATH.project("test", "Foo"))
Pkg.develop(PackageSpec(path = PATH.project()))
Pkg.activate(PATH.project("test"))
Pkg.develop(PackageSpec(path = PATH.project("test", "Foo")))

using Foo

@test Foo.CASTED_COMMANDS["main"].version == v"0.1.0"

rcfile_content = """
some words
# generated by Comonicon
# Julia bin PATH
export PATH="$(homedir())/.julia/bin:\$PATH"\n
# generated by Comonicon
# Julia autocompletion PATH
export FPATH="$(homedir())/.julia/completions:\$FPATH"
autoload -Uz compinit && compinit
"""

rcfile = tempname()
rm(rcfile, force = true)
write(rcfile, "some words")
write_path(rcfile, true, Dict())
@test strip(read(rcfile, String)) == strip(rcfile_content)

# using Comonicon.BuildTools: read_toml

# using Foo

# Foo.command_main(String[])

# read_toml(PATH.project("test", "Comonicon.toml"))

d = Dict(
    "name" => "foo",
    "download" => Dict("repo" => "Foo.jl", "host" => "github.com", "user" => "Roger-luo"),
    "install" =>
        Dict("optimize" => 2, "quiet" => false, "completion" => true, "compile" => "min"),
    "sysimg" => Dict(
        "filter_stdlibs" => false,
        "cpu_target" => "native",
        "incremental" => true,
        "path" => "deps/lib",
    ),
)

@test d == read_toml(Foo)

mock(create_sysimage) do plus
    @assert plus isa Mock
    Comonicon.build(Foo, true; bin = PATH.project("test", "bin"), quiet = true)

    Comonicon.build(Foo, false; bin = PATH.project("test", "bin"), quiet = true)
    Comonicon.install(Foo)
    @test isfile(PATH.project("test", "bin", "foo"))
    @test isfile(PATH.project("test", "bin", "foo.jl"))
end

@test ispath(PATH.project("test", "Foo", "deps"))

empty!(ARGS)
push!(ARGS, "sysimg")
mock(create_sysimage) do plus
    @assert plus isa Mock
    Comonicon.install(Foo; bin = PATH.project("test", "bin"), quiet = true)
end

@test isfile(PATH.project("test", "Foo", "deps", BuildTools.tarball_name("foo")))

Pkg.rm(PackageSpec(name = "Foo"))

# PR #48: check if the name of entry stays the same as config
module PR48
using ..Comonicon
@cast foo(x) = x
@cast goo(y) = y

@main name = "PR48"

end

configs =
    Dict("name" => "foo", "completion" => true, "quiet" => false, "compile" => "min", "optimize" => 2)

@testset "#48" begin
    @test_throws ErrorException Comonicon.BuildTools.validate_toml(PR48, configs)
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Entwine"
version = v"2.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/connormanning/entwine/archive/refs/tags/$version.tar.gz", "c12989e417182aa43a574d58eac284d4ea422abd5ac041504f55f635905e9b32"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/entwine-*

#this is needed for v2.1.0, but has been fixed on master, so can probably get rid of next version release
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/add-pdal-lasheader.patch

if [[ ${target} == *mingw* ]]; then
    #this is needed for v2.1.0, but has been fixed on master, so can probably get rid of next version release
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-lowercase-shlwapi.patch

    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-replace-msvc-find-shlwapi.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-remove-msvc-warning-options.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-remove-dupenv-branch.patch
fi

if [[ ${target} == *-linux-musl* ]]; then

    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/musl-remove-backtrace-detection.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/musl-hardcode-curl-link.patch

    if [[ ${target} == x86_64-* ]]; then
        # otherwise get "undefined reference to getrandom" error message on x86_64
        rm /usr/lib/libexpat*
    fi
fi

mkdir build
cd build

cmake .. -G Ninja \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_LIBRARY_PATH:FILEPATH="${libdir}" \
    -DCMAKE_INCLUDE_PATH:FILEPATH="${includedir}" \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DWITH_TESTS=OFF

ninja -j${nproc}
ninja install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libentwine", :libentwine),
    ExecutableProduct("entwine", :entwine)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="PDAL_jll", uuid="a8197b14-d70b-5660-b10f-8b1ebb62825c"))
    Dependency(PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7.1.0", julia_compat="1.6")

max_stars_repo_name: remi-garcia/Yggdrasil

max_issues_repo_name: remi-garcia/Yggdrasil

max_forks_repo_name: remi-garcia/Yggdrasil


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/21/2014

#ACASX implementation: based on CCAS interface to libCAS

module ACASX_CCAS_Impl

export
    addObserver,

    initialize,
    update,

    ACASX_CCAS,
    ACASXInput,
    ACASXOutput

import Compat.ASCIIString

using AbstractCollisionAvoidanceSystemImpl
using AbstractCollisionAvoidanceSystemInterfaces
using CommonInterfaces

using AbstractCASCoordImpl

using ACASXCommonImpl
using CASCoordination
using CASInterface
using CCAS
using RLESUtils, Observers

import CommonInterfaces.initialize
import CommonInterfaces.update

type ACASX_CCAS <: AbstractCollisionAvoidanceSystem
  my_id::Int64 #aircraft number
  max_intruders::Int64
  quant::Int64
  constants::Constants
  casShared::CASShared
  version::ASCIIString
  equipage::EQUIPAGE
  input::Input
  output::Output
  coord::AbstractCASCoord

  function ACASX_CCAS(aircraft_id::Int64, libcas::AbstractString, config_file::AbstractString,
                 num_aircraft::Int, coord::AbstractCASCoord, equipage::EQUIPAGE=EQUIPAGE_TCAS)
    cas = new()
    cas.my_id = aircraft_id
    cas.max_intruders = num_aircraft - 1
    cas.quant = equipage == EQUIPAGE_TCAS ? 25 : 100
    cas.constants = Constants(cas.quant, config_file, cas.max_intruders)
    cas.casShared = CASShared(libcas, cas.constants)
    cas.version = version(cas.casShared)
    cas.equipage = equipage

    @assert cas.max_intruders == max_intruders(cas.casShared) #Will fail if library did not open properly

    cas.coord = coord
    cas.input = Input(cas.max_intruders)
    cas.output = Output(cas.max_intruders)
    setRecord(cas.coord, cas.my_id,
              ACASXCoordRecord(cas.my_id, equipage, cas.quant, cas.max_intruders))
    reset(cas.casShared)
    return cas
  end
end

function update(cas::ACASX_CCAS, input::ACASXInput)
  ACASXCommonImpl.update_from_coord!(input, cas.coord, cas.my_id)
  if cas.equipage == EQUIPAGE_TCAS
    update!(cas.casShared, input, cas.output) #cas.output is modified in place
    ACASXCommonImpl.update_to_coord!(cas.coord, cas.my_id, cas.output)
  end

  cas.input = input #keep a record
  return cas.output
end

function initialize(cas::ACASX_CCAS)
  ACASXCommonImpl.initialize(cas)
  reset(cas.casShared)
end

end #module

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
function LinearGaussian(μ::Float64, σ::Float64)
    α = 5.0
    x = rand(:x, Normal(μ, σ))
    y = rand(:y, Normal(α * x, 1.0))
    return y
end

function LinearGaussianProposal()
    α = 10.0
    x = rand(:x, Normal(α * 3.0, 3.0))
end

@testset "Importance sampling" begin
    ret, cl = simulate(LinearGaussian, 0.0, 1.0)
    sel = Jaynes.selection([(:x, )])

    @testset "Linear Gaussian model" begin
        new, discard = Jaynes.metropolis_hastings(sel, cl)
    end

    @testset "Linear Gaussian proposal" begin
    end
end


text:
struct AffiliatedKeywordsWrapper{E <: Element} <: Element
    element::E
    keywords::Vector{AffiliatedKeyword}
end


text:
module StreamStats

using OffsetArrays
using Statistics
using StatsBase
using JSON

function monte_carlo_pi(data)
    """
    Calculate Monte Carlo Pi approximation for a given byte array.
    Keyword arguments:
    data -- data bytes
    """
    set_length::Int32 = round((length(data) / 2) - 1)
    r_square::Int16 = 128^2
    circle_surface::Float16 = 0

    for i in 1:set_length
        (data[i*2] - 128)
        if ((data[i*2] - 128) ^ 2 + (data[i*2+1] - 128) ^ 2) <= r_square
            circle_surface += 1
        end
    end
    return 4 * circle_surface / set_length
end

function pi_deviation(pi_value)
    """
    Returns an absolute percentage of difference between the provided Pi
    value and canonic.
    Keyword arguments:
    pi_value -- Pi value to be tested for difference
    """
    return abs(100 - (pi_value * 100 / pi))
end

function chi_squared(data)
    """Calculate Chi-Squared value for a given byte array.
    Keyword arguments:
    data -- data bytes
    """
    expected = length(data) / 256
    observed = OffsetVector([0 for i in 0:255], 0:255)

    for b in data
        observed[b] += 1
    end

    chi_squared = 0
    for o in observed
        chi_squared += (o - expected) ^ 2 / expected
    end

    return chi_squared
end

function get_all(data)
    """
    Returns a set of Statistical Information about the data
    Keywird arguments:
    data -- data bytes
    """
    carlo = monte_carlo_pi(data)
    dist = StatsBase.countmap(data)
    ddist = Dict{Int32,Int32}()
    for i in 0:255
        if Base.haskey(dist,i)
          ddist[i] = dist[i]
        else
          ddist[i] = 0
        end
    end

    orderd_dist =  sort(collect(dist), by=x->x[1])

    return JSON.json(Dict{String, Any}([
                "length"         => length(data),
                "chi-squared"    => chi_squared(data),
                "min"            => Int(Statistics.minimum(data)),
                "max"            => Int(Statistics.maximum(data)),
                "mean"           => Statistics.mean(data),
                "median"         => Statistics.median(data),
                "middle"         => Statistics.middle(data),
                "standard-dev"   => StatsBase.std(data),
                "kurtosis"       => StatsBase.kurtosis(data),
                "cor"            => Int(StatsBase.cor(data)),
                "skewness"       => StatsBase.skewness(data),
                "variation"      => StatsBase.variation(data),
                "sem"            => StatsBase.sem(data),
                "entropy"        => StatsBase.entropy(data),
                "cov"            => StatsBase.cov(data),
                "mad"            => StatsBase.mad(data, normalize=false),
                "monte-carlo-pi" => carlo,
                "pi-deviation"   => pi_deviation(carlo),
                "distribution"   => orderd_dist
                ])
            )
end

Base.@ccallable function julia_main()::Cint
    try
        main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function main()
    println(ARGS)
    for file in ARGS
        if !isfile(file)
            error("could not find file $file")
        end
        println(get_all(read(file)))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

end


text:
"""
A generic interface for two-players, symmetric, zero-sum games.
"""
module GameInterface

export AbstractPlayer

import ..Util
using ..Util: @unimplemented

#####
##### Types
#####

"""
    AbstractGame

Abstract base type for a game state.

# Constructors

Any subtype `Game` must implement `Base.copy` along with
the following constructors:

    Game()

Return the initial state of the game.

    Game(board, white_playing=true)

Return the unique state specified by a board and a current player.
By convention, the first player to play is called _white_ and
the other is called _black_.
"""
abstract type AbstractGame end

"""
    Board(Game::Type{<:AbstractGame})

Return the board type corresponding to `Game`.

Board objects must be persistent or appear as such.
"""
function Board(::Type{<:AbstractGame})
  @unimplemented
end

"""
    Action(Game::Type{<:AbstractGame})

Return the action type corresponding to `Game`.

Actions must be "symmetric" in the following sense:

```
available_actions(s) ==
  available_actions(Game(board_symmetric(s), !white_playing(s)))
```
"""
function Action(::Type{<:AbstractGame})
  @unimplemented
end

"""
    actions(::Type{<:AbstractGame})

Return the vector of all game actions.
"""
function actions(::Type{<:AbstractGame})
  @unimplemented
end

function Base.copy(::AbstractGame)
  @unimplemented
end

#####
##### Game functions
#####

"""
    white_playing(state::AbstractGame) :: Bool

Return `true` if white is to play and `false` otherwise.
"""
function white_playing(::AbstractGame)
  @unimplemented
end

"""
    white_reward(state::AbstractGame)

Return `nothing` if the game hasn't ended. Otherwise, return a
reward for the white player as a number between -1 and 1.
"""
function white_reward(::AbstractGame)
  @unimplemented
end

"""
    board(state::AbstractGame)

Return the game board.
"""
function board(::AbstractGame)
  @unimplemented
end

"""
    board_symmetric(state::AbstractGame)

Return the symmetric of the game board
(where the roles of black and white are swapped).
"""
function board_symmetric(::AbstractGame)
  @unimplemented
end

"""
    actions_mask(state::AbstractGame)

Return a boolean mask indicating what actions are available from `state`.

The following identities must hold:

    game_terminated(state) || any(actions_mask(state))
    length(actions_mask(state)) == length(actions(typeof(state)))
"""
function actions_mask(state::AbstractGame)
  @unimplemented
end

"""
    play!(state::AbstractGame, action)

Update the game state by making the current player perform `action`.
"""
function play!(state::AbstractGame, action)
  @unimplemented
end

"""
    heuristic_value(state::AbstractGame)

Return a heuristic estimate of the state value for the current player.

The given state must be nonfinal and returned values must belong to the
``(-∞, ∞)`` interval. Also, implementations of this function must be
antisymmetric in the sense that:
```
heuristic_value(s) ==
  - heuristic_value(Game(board_symmetric(s), white_playing(s)))
```

This function is not needed by AlphaZero but it is useful for building
baselines such as minmax players.
"""
function heuristic_value(state::AbstractGame)
  @unimplemented
end

#####
##### Symmetries
#####

"""
    symmetries(::Type{G}, board) where {G <: AbstractGame}

Return the vector of all pairs `(b, σ)` where:
  - `b` is the image of `board` by a nonidentical symmetry
  - `σ` is the associated actions permutation, as an integer vector of
     size `num_actions(Game)`.

A default implementation is provided that returns an empty vector.
"""
function symmetries(::Type{G}, board) where {G <: AbstractGame}
  return Tuple{Board(G), Vector{Int}}[]
end

function test_symmetry(Game, board, (symboard, aperm))
  syms = symmetries
  mask = actions_mask(Game(board))
  symmask = actions_mask(Game(symboard))
  v = falses(length(symmask))
  v[mask] .= true
  v = v[aperm]
  return all(v[symmask]) && !any(v[.~symmask])
end

#####
##### Machine learning interface
#####

"""
    vectorize_board(::Type{<:AbstractGame}, board) :: Vector{Float32}

Return a vectorized representation of a board.
"""
function vectorize_board(::Type{<:AbstractGame}, board)
  @unimplemented
end

#####
##### Interface for interactive exploratory tools
#####

"""
    action_string(::Type{<:AbstractGame}, action) :: String

Return a human-readable string representing the provided action.
"""
function action_string(::Type{<:AbstractGame}, action)
  @unimplemented
end

"""
    parse_action(::Type{<:AbstractGame}, str::String)

Return the action described by string `str` or `nothing`
if `str` does not denote a valid action.
"""
function parse_action(::Type{<:AbstractGame}, ::String)
  @unimplemented
end

"""
    read_state(::Type{G}) where G <: AbstractGame :: Union{G, Nothing}

Read a state description from the standard input.
Return the corresponding state or `nothing` in case of an invalid input.
"""
function read_state(::Type{<:AbstractGame})
  @unimplemented
end

"""
    print_state(state::AbstractGame)

Print a state on the standard output.
"""
function print_state(::AbstractGame)
  @unimplemented
end

#####
##### Derived functions
#####

"""
    game_terminated(::AbstractGame)

Return a boolean indicating whether or not a game is in a terminal state.
"""
game_terminated(state) = !isnothing(white_reward(state))

"""
    num_actions(::Type{G})

Return the total number of actions associated with a game.
"""
num_actions(::Type{G}) where G = length(actions(G))

"""
    available_actions(state::AbstractGame)

Return the vector of all available actions in a given state.
"""
function available_actions(state::AbstractGame)
  Game = typeof(state)
  mask = actions_mask(state)
  return actions(Game)[mask]
end

"""
    board_dim(::Type{G})

Return a tuple that indicates the shape of a vectorized board representation.
"""
board_dim(::Type{G}) where G = size(vectorize_board(G, board(G())))

function canonical_board(state)
  return white_playing(state) ? board(state) : board_symmetric(state)
end

"""
    random_symmetric_state(::AbstractGame)

Return a fresh new state that is the image of the given state by a random
symmetry (see [`symmetries`](@ref)).
"""
function random_symmetric_state(state::Game) where {Game <: AbstractGame}
  bsym, _ = rand(symmetries(Game, board(state)))
  return Game(bsym, white_playing(state))
end

function board_memsize(::Type{G}) where G
  return Base.summarysize(board(G()))
end

symmetric_reward(r::Real) = -r

end

"""
    GameType(T)

Return the `Game` type associated with an object.
"""
function GameType(T)
  @unimplemented
end


text:
using Econometrics, Test
using Econometrics: ConvergenceException, Hermitian
using CSV, RDatasets
using MLJBase: MLJBase
using MLJModelInterface: MLJModelInterface

for file in ["exceptions", "linear_models", "mlogit", "ologit", "formula_display", "mlj", "docs"]
    include("$file.jl")
end


text:
using StressTest
using Test, Random

@testset "busywait(seconds) duration within 1% of target duration" begin
    rng = Random.seed!(1234)
    for duration in (rand(5) .* 2)
        e = @elapsed busywait(duration)
        @test isapprox(e,duration,rtol=0.01)
    end
end


text:
module FluidTensors

import LinearAlgebra
using LinearAlgebra

export Vec, SymTen, symouter, antisymouter, Lie, traceless, AbstractVecArray, VecArray, AbstractSymTenArray, AbstractSymTrTenArray, SymTenArray, SymTrTenArray, AbstractAntiSymTenArray, AntiSymTenArray, HomogeneousArray, eig, eigvec, anglecos, stress_state
export prop_decomp, inph_decomp, Ten, outer, ⊗, AntiSymTen, square, symmetric, antisymmetric

abstract type AbstractVec{T} <: AbstractArray{T,1} end
abstract type AbstractTen{T} <: AbstractArray{T,2} end

include("vectors.jl")
include("tensors.jl")
include("utils.jl")
include("eigen.jl")
include("decompositions.jl")

end


text:
function centroid(
    verts::AbstractVector{<:AbstractPoint{3,T1}},
) where {T1 <: Real}
    centroid_point = reduce(+, verts) / length(verts)
    return centroid_point
end
function centroid(mesh)
    verts = coordinates(mesh)
    centroid_point = centroid(verts)
end
function decompose_vertex_and_face(mesh::Mesh)
    verts = coordinates(mesh)
    triangles = faces(mesh)
    return verts, triangles
end
function angle_between_vector(v1, v2)
    v1 = v1 / norm(v1)
    v2 = v2 / norm(v2)
    cos_theta = dot(v1, v2)
    if cos_theta >= 1.0
        return 0.0
    end
    theta = acos(cos_theta)
    return theta
end

function calculate_area_bins(verts::AbstractVector{<:AbstractPoint{3,T1}},
            triangles::AbstractVector{<:TriangleFace},
            ) where T1 <: AbstractFloat
    areas_array = area.(Ref(verts), triangles)
    cumsum_area = similar(areas_array)
    cumsum_area = cumsum!(cumsum_area, areas_array)
    total_area = cumsum_area[end]
    bins = cumsum_area / total_area
    return bins
end

function triangle_area(p1::T, p2::T, p3::T) where {T <: AbstractVector}
    p1p2 = p2 - p1
    p1p3 = p3 - p1
    area = 0.5 * norm(cross(p1p2, p1p3))
    return area
end

function tetrahedron_volume(v1::T, v2::T, v3::T, v4::T) where {T <: AbstractVector}
    a = v1 - v4
    b = v2 - v4
    c = v3 - v4
    vol = norm(a ⋅ cross(b, c)) / 6
    return vol
end


text:
__precompile__()

module TestModule
	export special_add

	include("add.jl")

	function special_add(a, b)
		return add(a, b) + 1
	end
end


text:
@testset "extensions" begin
    c = SVector(99.0, 99.0, 99.0)
    e = SVector(1.0, 2.0, 3.0)

    aabb = OBB(c, e, one(SMatrix{3,3,Float64,9}))
    @test 8 * 6.0 == volume(aabb)
    @test 4 * 22.0 == area(aabb)
end


text:
"""
    DNNDetector(d = 0,
                metric = Euclidean(),
                algorithm = :kdtree,
                leafsize = 10,
                reorder = true,
                parallel = false)

Anomaly score based on the number of neighbors in a hypersphere of radius `d`. Knorr et al. [1] directly converted the
resulting outlier scores to labels, thus this implementation does not fully reflect the approach from the paper.

Parameters
----------
    d::Real
The hypersphere radius used to calculate the global density of an instance.

$KNN_PARAMS

Examples
--------
$(SCORE_UNSUPERVISED("DNNDetector"))

References
----------
[1] Knorr, Edwin M.; Ng, Raymond T. (1998): Algorithms for Mining Distance-Based Outliers in Large Datasets.
"""
OD.@detector mutable struct DNNDetector <: UnsupervisedDetector
    metric::DI.Metric = DI.Euclidean()
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :balltree))
    static::Union{Bool,Symbol} = :auto::(_ in (true, false, :auto))
    leafsize::Integer = 10::(_ ≥ 0)
    reorder::Bool = true
    parallel::Bool = false
    d::Real = 0::(_ > 0) # warns if `d` is not set
end

struct DNNModel <: DetectorModel
    tree::NN.NNTree
end

function OD.fit(detector::DNNDetector, X::Data; verbosity)::Fit
    X = prepare_data(X, detector.static)

    # create the specified tree
    tree = @tree detector X

    # use tree to calculate distances
    scores = dnn_others(NN.inrange(tree, X, detector.d))
    return DNNModel(tree), scores
end

function OD.transform(detector::DNNDetector, model::DNNModel, X::Data)::Scores
    X = prepare_data(X, detector.static)

    if detector.parallel
        return dnn_parallel(model.tree, X, detector.d)
    else
        return dnn(NN.inrange(model.tree, X, detector.d))
    end
end

@inline function dnn(idxs::AbstractVector{<:AbstractVector})::Scores
    # Helper function to reduce the instances to a global density score.
    1 ./ (length.(idxs) .+ 0.1) # min score = 0, max_score = 10
end

@inline function dnn_others(idxs::AbstractVector{<:AbstractVector})::Scores
    # Remove the (self) point previously added when fitting the tree, otherwise during `fit`, that point would always
    # be included in the density estimation
    1 ./ (length.(idxs) .- 0.9) # 1 - 0.1
end


text:
#### NullNode ####

type NullRenderer <: AudioRenderer end
typealias NullNode AudioNode{NullRenderer}
export NullNode

function render(node::NullRenderer, device_input::AudioBuf, info::DeviceInfo)
    # TODO: preallocate buffer
    return zeros(info.buf_size)
end

#### SinOsc ####

# Generates a sin tone at the given frequency

type SinOscRenderer{T<:Union(Float32, AudioNode)} <: AudioRenderer
    freq::T
    phase::Float32
    buf::AudioBuf

    function SinOscRenderer(freq)
        new(freq, 0.0, AudioSample[])
    end
end

typealias SinOsc AudioNode{SinOscRenderer}
SinOsc(freq::Real) = SinOsc(SinOscRenderer{Float32}(freq))
SinOsc(freq::AudioNode) = SinOsc(SinOscRenderer{AudioNode}(freq))
SinOsc() = SinOsc(440)
export SinOsc

function render(node::SinOscRenderer{Float32}, device_input::AudioBuf,
        info::DeviceInfo)
    if length(node.buf) != info.buf_size
        resize!(node.buf, info.buf_size)
    end
    outbuf = node.buf
    phase = node.phase
    freq = node.freq
    # make sure these are Float32s so that we don't allocate doing conversions
    # in the tight loop
    pi2::Float32 = 2pi
    phase_inc::Float32 = 2pi * freq / info.sample_rate
    i::Int = 1
    while i <= info.buf_size
        outbuf[i] = sin(phase)
        phase = (phase + phase_inc) % pi2
        i += 1
    end
    node.phase = phase
    return outbuf
end

function render(node::SinOscRenderer{AudioNode}, device_input::AudioBuf,
        info::DeviceInfo)
    freq = render(node.freq, device_input, info)::AudioBuf
    block_size = min(length(freq), info.buf_size)
    if(length(node.buf) != block_size)
        resize!(node.buf, block_size)
    end
    outbuf = node.buf

    phase::Float32 = node.phase
    pi2::Float32 = 2pi
    phase_step::Float32 = 2pi/(info.sample_rate)
    i::Int = 1
    while i <= block_size
        outbuf[i] = sin(phase)
        phase = (phase + phase_step*freq[i]) % pi2
        i += 1
    end
    node.phase = phase
    return outbuf
end

#### AudioMixer ####

# Mixes a set of inputs equally

type MixRenderer <: AudioRenderer
    inputs::Vector{AudioNode}
    buf::AudioBuf

    MixRenderer(inputs) = new(inputs, AudioSample[])
    MixRenderer() = MixRenderer(AudioNode[])
end

typealias AudioMixer AudioNode{MixRenderer}
export AudioMixer

function render(node::MixRenderer, device_input::AudioBuf, info::DeviceInfo)
    if length(node.buf) != info.buf_size
        resize!(node.buf, info.buf_size)
    end
    mix_buffer = node.buf
    n_inputs = length(node.inputs)
    i = 1
    max_samples = 0
    fill!(mix_buffer, 0)
    while i <= n_inputs
        rendered = render(node.inputs[i], device_input, info)::AudioBuf
        nsamples = length(rendered)
        max_samples = max(max_samples, nsamples)
        j::Int = 1
        while j <= nsamples
            mix_buffer[j] += rendered[j]
            j += 1
        end
        if nsamples < info.buf_size
            deleteat!(node.inputs, i)
            n_inputs -= 1
        else
            i += 1
        end
    end
    if max_samples < length(mix_buffer)
        return mix_buffer[1:max_samples]
    else
        # save the allocate and copy if we don't need to
        return mix_buffer
    end
end

Base.push!(mixer::AudioMixer, node::AudioNode) = push!(mixer.renderer.inputs, node)

#### Gain ####
type GainRenderer{T<:Union(Float32, AudioNode)} <: AudioRenderer
    in1::AudioNode
    in2::T
    buf::AudioBuf

    GainRenderer(in1, in2) = new(in1, in2, AudioSample[])
end

function render(node::GainRenderer{Float32},
                device_input::AudioBuf,
                info::DeviceInfo)
    input = render(node.in1, device_input, info)::AudioBuf
    if length(node.buf) != length(input)
        resize!(node.buf, length(input))
    end
    i = 1
    while i <= length(input)
        node.buf[i] = input[i] * node.in2
        i += 1
    end
    return node.buf
end

function render(node::GainRenderer{AudioNode},
                device_input::AudioBuf,
                info::DeviceInfo)
    in1_data = render(node.in1, device_input, info)::AudioBuf
    in2_data = render(node.in2, device_input, info)::AudioBuf
    block_size = min(length(in1_data), length(in2_data))
    if length(node.buf) != block_size
        resize!(node.buf, block_size)
    end
    i = 1
    while i <= block_size
        node.buf[i] = in1_data[i] * in2_data[i]
        i += 1
    end
    return node.buf
end

typealias Gain AudioNode{GainRenderer}
Gain(in1::AudioNode, in2::Real) = Gain(GainRenderer{Float32}(in1, in2))
Gain(in1::AudioNode, in2::AudioNode) = Gain(GainRenderer{AudioNode}(in1, in2))
export Gain

#### Offset ####
type OffsetRenderer <: AudioRenderer
    in_node::AudioNode
    offset::Float32
    buf::AudioBuf

    OffsetRenderer(in_node, offset) = new(in_node, offset, AudioSample[])
end

function render(node::OffsetRenderer, device_input::AudioBuf, info::DeviceInfo)
    input = render(node.in_node, device_input, info)::AudioBuf
    if length(node.buf) != length(input)
        resize!(node.buf, length(input))
    end
    i = 1
    while i <= length(input)
        node.buf[i] = input[i] + node.offset
        i += 1
    end
    return node.buf
end

typealias Offset AudioNode{OffsetRenderer}
export Offset

#### Array Player ####

# Plays a AudioBuf by rendering it out piece-by-piece

type ArrayRenderer <: AudioRenderer
    arr::AudioBuf
    arr_index::Int
    buf::AudioBuf

    ArrayRenderer(arr::AudioBuf) = new(arr, 1, AudioSample[])
end

typealias ArrayPlayer AudioNode{ArrayRenderer}
export ArrayPlayer

function render(node::ArrayRenderer, device_input::AudioBuf, info::DeviceInfo)
    range_end = min(node.arr_index + info.buf_size-1, length(node.arr))
    block_size = range_end - node.arr_index + 1
    if length(node.buf) != block_size
        resize!(node.buf, block_size)
    end
    copy!(node.buf, 1, node.arr, node.arr_index, block_size)
    node.arr_index = range_end + 1
    return node.buf
end

# Allow users to play a raw array by wrapping it in an ArrayPlayer
function play(arr::AudioBuf, args...)
    player = ArrayPlayer(arr)
    play(player, args...)
end

# If the array is the wrong floating type, convert it
function play{T <: FloatingPoint}(arr::Array{T}, args...)
    arr = convert(AudioBuf, arr)
    play(arr, args...)
end

# If the array is an integer type, scale to [-1, 1] floating point

# integer audio can be slightly (by 1) more negative than positive,
# so we just scale so that +/- typemax(T) becomes +/- 1
function play{T <: Signed}(arr::Array{T}, args...)
    arr = arr / typemax(T)
    play(arr, args...)
end

function play{T <: Unsigned}(arr::Array{T}, args...)
    zero = (typemax(T) + 1) / 2
    range = floor(typemax(T) / 2)
    arr = (arr .- zero) / range
    play(arr, args...)
end

#### Noise ####

type WhiteNoiseRenderer <: AudioRenderer end
typealias WhiteNoise AudioNode{WhiteNoiseRenderer}
export WhiteNoise

function render(node::WhiteNoiseRenderer, device_input::AudioBuf, info::DeviceInfo)
    return rand(AudioSample, info.buf_size) .* 2 .- 1
end

#### AudioInput ####

# Renders incoming audio input from the hardware

type InputRenderer <: AudioRenderer
    channel::Int
    InputRenderer(channel::Integer) = new(channel)
    InputRenderer() = new(1)
end

function render(node::InputRenderer, device_input::AudioBuf, info::DeviceInfo)
    @assert size(device_input, 1) == info.buf_size
    return device_input[:, node.channel]
end

typealias AudioInput AudioNode{InputRenderer}
export AudioInput

#### LinRamp ####

type LinRampRenderer <: AudioRenderer
    key_samples::Array{AudioSample}
    key_durations::Array{Float32}

    duration::Float32
    buf::AudioBuf

    LinRampRenderer(start, finish, dur) = LinRampRenderer([start,finish], [dur])

    LinRampRenderer(key_samples, key_durations) =
        LinRampRenderer(
            [convert(AudioSample,s) for s in key_samples],
            [convert(Float32,d) for d in key_durations]
        )

    function LinRampRenderer(key_samples::Array{AudioSample}, key_durations::Array{Float32})
        @assert length(key_samples) == length(key_durations) + 1
        new(key_samples, key_durations, sum(key_durations), AudioSample[])
    end
end

typealias LinRamp AudioNode{LinRampRenderer}
export LinRamp

function render(node::LinRampRenderer, device_input::AudioBuf, info::DeviceInfo)
    # Resize buffer if (1) it's too small or (2) we've hit the end of the ramp
    ramp_samples::Int = int(node.duration * info.sample_rate)
    block_samples = min(ramp_samples, info.buf_size)
    if length(node.buf) != block_samples
        resize!(node.buf, block_samples)
    end

    # Fill the buffer as long as there are more segments
    dt::Float32 = 1/info.sample_rate
    i::Int = 1
    while i <= length(node.buf) && length(node.key_samples) > 1

        # Fill as much of the buffer as we can with the current segment
        ds::Float32 = (node.key_samples[2] - node.key_samples[1]) / node.key_durations[1] / info.sample_rate
        while i <= length(node.buf)
            node.buf[i] = node.key_samples[1]
            node.key_samples[1] += ds
            node.key_durations[1] -= dt
            node.duration -= dt
            i += 1

            # Discard segment if we're finished
            if node.key_durations[1] <= 0
                if length(node.key_durations) > 1
                    node.key_durations[2] -= node.key_durations[1]
                end
                shift!(node.key_samples)
                shift!(node.key_durations)
                break
            end
        end
    end

    return node.buf
end


text:
module Jazz
using MIDI, MusicManipulations

"""
```julia
average_swing_ratio(notes::Notes, asr_method::String)
```
Calculate the average swing ratio of given `notes` array. Return the average
swing ratio and the associated standard deviation.
## Methods
### `"swung8s"`
Classify each note with the grid `[0, 2//3, 1]`. Use the notes with
index 2 to calculate the swing ratio, and then average for all.
### `"triplets"`
Classify each note with the grid `[0, 1//3, 2//3, 1]`. Use the notes with
index 3 to calculate the swing ratio, and then average for all.
"""
function average_swing_ratio(notes::Notes, asr_method::String)

    tpq = notes.tpq
    sr = Float64[]

    # Create array of Swing Notes
    if asr_method == "swung8s"
        s = classify(notes, [0, 2//3, 1])
        swingnotes = notes.notes[s .== 2]
    elseif asr_method == "triplets"
        s = classify(notes, [0, 1//3, 2//3, 1])
        swingnotes = notes.notes[s .== 3]
    end

    # Create Swing ratios array
    for i in 1:length(swingnotes)
        pos = Int64(swingnotes[i].position)
        posmod = mod(pos, tpq)
        push!(sr, posmod/(tpq-posmod))
    end

    return mean(sr), std(sr)
end

function inbetween_portion(notes::Notes)
  clas = classify(notes, [0, 1//3, 2//3, 1])
  count(i -> i == 2, clas)/length(notes)
end

end


text:
@testset "$N" for (N, FR, TI) in (("FastQ32", FastQ32, Int32), ("FastQ64", FastQ64, Int64), ("FastQBig", FastQBig, BigInt))
  @eval begin
    @test $FR(1//1) == 1
    @test $FR(2//2) == 1
    @test $FR(1//1) == 1//1
    @test $FR(2//2) == 1//1
    @test $FR(2//4) == $FR(8//16)
    @test $FR(1//2) + $FR(1//2) == 1
    @test $FR((-1)//3) == $FR(-(1//3)) == -(1//3)
    @test $FR(1//2) + $FR(3//4) == $FR(5//4) == 5//4
    @test $FR(1//3) * $FR(3//4) == $FR(1//4) == 1//4
    @test $FR(1//2) / $FR(3//4) == $FR(2//3)
    @test $FR(-1//1) / $FR(-1//1) == $FR(1//1)

    #=
    @test_throws OverflowError -FR(typemin(TI)//1)
    @test_throws OverflowError FR(typemax(TI)//3) + 1
    @test_throws OverflowError FR(typemax(TI)//3) * 2
    @test_throws OverflowError FR(1//2)^63
    =#

    for a = -5:5, b = -5:5
        if b == 0
            # @test_throws ArgumentError round(Integer, $FR(a, b))
            continue
        end
        if ispow2(b)
            @test $FR(a//b) == $FR(a/b)
            @test convert(FastRational{$TI}, a/b) == $FR(a//b)
        end
        # @test rationalize(a/b) == a//b
        @test $FR(a//b) == a//b
        @test round(Integer,$FR(a//b)) == round(Integer,a/b)
        for c = -5:5
            @test ($FR(a//b) == c) == (a/b == c)
            @test ($FR(a//b) != c) == (a/b != c)
            @test ($FR(a//b) <= c) == (a/b <= c)
            @test ($FR(a//b) <  c) == (a/b <  c)
            @test ($FR(a//b) >= c) == (a/b >= c)
            @test ($FR(a//b) >  c) == (a/b >  c)
            for d = -5:5
                if d == 0; continue; end
                @test ($FR(a//b) == $FR(c//d)) == (a/b == c/d)
                @test ($FR(a//b) != $FR(c//d)) == (a/b != c/d)
                @test ($FR(a//b) <= $FR(c//d)) == (a/b <= c/d)
                @test ($FR(a//b) <  $FR(c//d)) == (a/b <  c/d)
                @test ($FR(a//b) >= $FR(c//d)) == (a/b >= c/d)
                @test ($FR(a//b) >  $FR(c//d)) == (a/b >  c/d)
            end
        end
    end

    @test 0.5 == $FR(1//2)
    @test 0.1 ≈ float($FR(1//10))
    @test 1/3 ≈ float($FR(1//3))
    @test !($FR(1//3) < 1/3)
    @test -1/3 < $FR(1//3)
    @test -1/3 ≈ float($FR(-1//3))
    @test 1/3 > $FR(-1//3)
    @test 1/5 ≈ float($FR(1//5))

    # PR 29561
    @test abs(one(FastRational{$TI})) == one(FastRational{$TI})
    @test abs(-one(FastRational{$TI})) == one(FastRational{$TI})
  end
end

@testset "$N" for (N, FR, TI) in (("FastQ32 typemax/min", FastQ32, Int32), ("FastQ64 typemax/min", FastQ64, Int64))
  @eval begin

    @test $FR(typemax($TI)//one($TI)) * $FR(one($TI)//typemax($TI)) == 1
    @test $FR(typemax($TI)//one($TI)) / $FR(typemax($TI)//one($TI)) == 1
    @test $FR(one($TI)//typemax($TI)) / $FR(one($TI)//typemax($TI)) == 1

    #>>>FIXME @test_throws InexactError -FR(typemin(TI)//1)
    @test_throws InexactError $FR(typemax($TI)//3) + $TI(1)
    @test_throws InexactError $FR(typemax($TI)//3) * $TI(2)
    @test_throws InexactError $FR(1//2)^63

  end
end


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD

function ColoringTheNode(i, j, Network, Num)
    check = false
    if Network[i, j] == 0
        Network[i, j] = Num
        check = true
        Num += 1
    end
    return Network, Num, check
end

function FindNeighbors(i, j, Network, dim)
    Numbereds = 0
    colorList = []
    for neighbor in ([0,1], [1,0], [0,-1], [-1,0])
        if j + neighbor[2] >= 1 && j + neighbor[2] <= dim && Network[([i, j] + neighbor)...] > 0
            push!(colorList,Network[([i, j] + neighbor)...])
            Numbereds += 1
        end
    end
    push!(colorList, Network[i, j])
    return colorList, Numbereds #colorList contain the node and its neighbors
end

function JoinColors(Network, dim, colorList)
    mincolor = min(colorList...)
    for i in 1:dim
        for j in 1:dim
            if Network[i, j] in colorList
                Network[i, j] = mincolor
            end
        end
    end
    return Network
end

function InitialNetwork(dim, P)
    Network = sample([-1,0], Weights([1-P, P]),(dim,dim))
    return Network
end

@inline function PercolationCheck(dim, P)
    Num = 2
    maxnum = 100000
    Network = InitialNetwork(dim, P)
    Network[1,:] = ones(dim)
    Network[end,:] = fill(maxnum,dim)
    for i in 2:dim-1
        for j in 1:dim
            Network, Num, check = ColoringTheNode(i, j, Network, Num)
            if check == true
                colorList, Numbereds = FindNeighbors(i, j, Network, dim)
                if Numbereds == 0
                    continue
                elseif Numbereds == 1
                    Network[i, j] = colorList[1]
                else
                    Network = JoinColors(Network, dim, colorList)
                end
            end
        end
    end
    if 1 in Network[end,:]
        return 1
    else
        return 0
    end
end

a = time()
PercolationCheck(100,0.6)
time()-a
Plist = hcat(0:0.02:1)
Avglist = []
STDlist = []
TotData = []
# PercolationCheck(40,0.55,1)
for P in Plist
    print("\r$P")
    AllRuns = []
    for i in 1:5000
        push!(AllRuns, PercolationCheck(10,P))
    end
    push!(TotData, AllRuns)
    push!(STDlist, std(AllRuns))
    push!(Avglist, mean(AllRuns))
end

save("10x10.jld", "data", TotData)

scatter(Plist, Avglist, yerr = STDlist, legend = nothing, xlabel = L"P", ylabel = L"Avg\ Q", title = L"Percolation \ for \ 10\times10 \ Network\ (5000 \ runs)")
savefig("../../Figs/Q5/10x10.pdf")


text:
using Base.Test
using RNN

#
# test initialisation
#
rnn = RNN_t(1,2,3)
@test size(rnn.w) == (6,6)
@test size(rnn.a) == (6,1)
@test size(rnn.b) == (6,1)
@test size(rnn.o) == (6,1)
@test rnn.ni      == 1
@test rnn.no      == 2
@test rnn.nh      == 3
@test rnn.n       == 6

#
# test modify synapse on a net without synapses
#
rnn = RNN_t(1,2,3)
rnn_modify_synapse!(rnn, 1.1, 100.0, 10.0)
for i=1:rnn.n
  for j=1:rnn.n
    @test rnn.w[i,j] == 0.0
  end
end

#
# test add synapses
#
rnn = RNN_t(1,2,3)
rnn_add_synapse!(rnn, 1.1, 1.0)
for i=1:rnn.n
  for j=1:rnn.n
    if i <= rnn.ni
      @test rnn.w[i,j] ==  0.0
    else
      @test rnn.w[i,j] !=  0.0
    end
    @test rnn.w[i,j] <=  1.0
    @test rnn.w[i,j] >= -1.0
  end
end

#
# test remove synapse
#
rnn = RNN_t(1,2,3)
rnn_remove_synapse!(rnn, 1.1)
for i=1:rnn.n
  for j=1:rnn.n
    @test rnn.w[i,j] == 0.0
  end
end

#
# test modify synapse
#
rnn1 = RNN_t(1,2,3)
rnn_add_synapse!(rnn1, 1.1, 1.0)
rnn2 = rnn_copy(rnn1)
rnn_modify_synapse!(rnn1, 1.1, 5.0, 10.0)
for i=rnn.ni+1:rnn.n
  for j=1:rnn.n
    @test rnn1.w[i,j] != rnn2.w[i,j]
  end
end
for i=1:rnn.ni
  for j=1:rnn.n
    @test rnn1.w[i,j] == 0.0
    @test rnn2.w[i,j] == 0.0
  end
end

#
# test add neuron
#
rnn = RNN_t(1,2,3)
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 4)
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 4)
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 4)
@test rnn.n       == 7
@test rnn.ni      == 1
@test rnn.no      == 2
@test rnn.nh      == 4
@test size(rnn.w) == (7,7)
@test size(rnn.a) == (7,1)
@test size(rnn.b) == (7,1)
@test size(rnn.o) == (7,1)
for i=rnn.ni+1:rnn.n
  @test rnn.w[i,rnn.n] != 0.0
end
for i=1:rnn.ni
  @test rnn.w[i,rnn.n] == 0.0
end
for i=1:rnn.n
  @test rnn.w[rnn.n,i] != 0.0
end
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 8)
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 8)
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 8)
rnn_add_neuron!(rnn, 1.0, 1.0, 2.0, 1.0, 8)
@test rnn.n       == 11
@test rnn.ni      == 1
@test rnn.no      == 2
@test rnn.nh      == 8
@test size(rnn.w) == (11,11)
@test size(rnn.a) == (11,1)
@test size(rnn.b) == (11,1)
@test size(rnn.o) == (11,1)

#
# test modify bias
#
rnn = RNN_t(1,2,3)
rnn_modify_bias!(rnn, 1.0, 1.0, 2.0)
for i=1:rnn.ni
  @test rnn.b[i] == 0.0
end
for i=rnn.ni+1:rnn.n
  @test rnn.b[i] != 0.0
end

#
# test remove neuron
#
rnn1 = RNN_t(1,2,3)
rnn_modify_synapse!(rnn1, 1.1, 5.0, 10.0)
rnn_modify_bias!(rnn1, 1.0, 1.0, 2.0)
rnn2 = rnn_copy(rnn1)
rnn_remove_neuron!(rnn1, 1.1)
rnn_add_neuron!(rnn1, 1.0, 1.0, 2.0, 1.0, 20)
rnn_add_neuron!(rnn1, 1.0, 1.0, 2.0, 1.0, 20)
rnn_add_neuron!(rnn1, 1.0, 1.0, 2.0, 1.0, 20)
rnn_add_neuron!(rnn1, 1.0, 1.0, 2.0, 1.0, 20)
rnn_add_neuron!(rnn1, 1.0, 1.0, 2.0, 1.0, 20)
rnn_remove_neuron!(rnn1, 1.1)
rnn_remove_neuron!(rnn1, 1.1)
rnn_remove_neuron!(rnn1, 1.1)
@test rnn1.n          == 3
@test rnn1.ni         == 1
@test rnn1.no         == 2
@test rnn1.nh         == 0
@test size(rnn1.w)    == (3,3)
@test size(rnn1.a)    == (3,1)
@test size(rnn1.b)    == (3,1)
@test size(rnn1.o)    == (3,1)
@test rnn1.w[1:3,1:3] == rnn2.w[1:3,1:3]
@test rnn1.b[1:3]     == rnn2.b[1:3]

#
# test update function
#
rnn = RNN_t(1,2,0)
rnn.w[1,1] = 0.0
rnn.w[2,1] = 2.0
rnn.w[2,2] = 0.0
rnn.w[1,2] = 0.0
rnn.b[1]   = 1.0
rnn.b[2]   = 0.0
rnn.b[3]   = 0.0
rnn.a[1]   = 0.0
rnn.a[2]   = 0.0
rnn.a[3]   = 0.0
rnn.o[1]   = 0.0
rnn.o[2]   = 0.0
rnn.o[3]   = 0.0
rnn_update!(rnn)
@test rnn.b[1] == 1.0
@test rnn.b[2] == 0.0
@test rnn.b[3] == 0.0
@test rnn.a[1] == 1.0
@test rnn.a[2] == 0.0
@test rnn.a[3] == 0.0
@test rnn.o[1] == 1.0
@test rnn.o[2] == 0.0
@test rnn.o[3] == 0.0
rnn_update!(rnn)
@test rnn.b[1] == 1.0
@test rnn.b[2] == 0.0
@test rnn.b[3] == 0.0
@test rnn.a[1] == 1.0
@test rnn.a[2] == 2.0
@test rnn.a[3] == 0.0
@test rnn.o[1] == 1.0
@test rnn.o[2] == tanh(2.0)
@test rnn.o[3] == 0.0

ni1 = 20 + int(rand() + 50)
no1 = 20 + int(rand() + 50)
nh1 = 20 + int(rand() + 50)
ni2 = 20 + int(rand() + 50)
no2 = 20 + int(rand() + 50)
nh2 = 20 + int(rand() + 50)

rnn1   = RNN_t(ni1, no1, nh1)
rnn1.w = rand(size(rnn1.w))
rnn2   = RNN_t(ni2, no2, nh2)
rnn2.w = rand(size(rnn2.w))

rnn = rnn_merge_two(rnn1, rnn2)

@test rnn.ni == ni1 + ni2
@test rnn.no == no1 + no2
@test rnn.nh == nh1 + nh2

for src=1:rnn.n
  for dest=1:rnn.n

    # inputs to input neurons must all be zero
    if dest <= rnn.ni && src <= rnn1.ni + rnn2.ni
      @test rnn.w[dest, src] == 0.0
    end

    #
    # RNN 1 Tests
    #

    # rnn1 input to output
    if dest > rnn.ni && dest <= rnn.ni + rnn1.no && src <= rnn1.ni
      d = dest - rnn.ni + rnn1.ni
      @test rnn.w[dest, src] == rnn1.w[d, src]
    end

    # rnn1 input to hidden
    if dest > rnn.ni + rnn.no && dest <= rnn.ni + rnn.no + rnn1.nh && src <= rnn1.ni
      d = dest - rnn.ni + rnn1.ni - rnn.no + rnn1.no
      @test rnn.w[dest, src] == rnn1.w[d, src]
    end

    # rnn1 output to output
    if dest > rnn.ni && dest <= rnn.ni + rnn1.no &&
       src  > rnn.ni && src  <= rnn.ni + rnn1.no
      d = dest - rnn.ni + rnn1.ni
      s = src  - rnn.ni + rnn1.ni
      @test rnn.w[dest, src] == rnn1.w[d, s]
    end

    # rnn1 output to hidden
    if dest > rnn.ni + rnn.no && dest <= rnn.ni + rnn.no + rnn1.nh &&
       src  > rnn.ni          && src  <= rnn.ni + rnn1.no
      d = dest - rnn.ni + rnn1.ni - rnn.no + rnn1.no
      s = src  - rnn.ni + rnn1.ni
      @test rnn.w[dest, src] == rnn1.w[d, s]
    end

    # rnn1 hidden to output
    if dest > rnn.ni          && dest <= rnn.ni + rnn1.no &&
       src  > rnn.ni + rnn.no && src  <= rnn.ni + rnn.no + rnn1.nh
      d = dest - rnn.ni + rnn1.ni
      s = src  - rnn.ni + rnn1.ni - rnn.no + rnn1.no
      @test rnn.w[dest, src] == rnn1.w[d, s]
    end

    #= # rnn1 hidden to hidden =#
    if dest > rnn.ni + rnn.no && dest <= rnn.ni + rnn.no + rnn1.nh &&
       src  > rnn.ni + rnn.no && src  <= rnn.ni + rnn.no + rnn1.nh
      d = dest - rnn.ni - rnn.no + rnn1.ni + rnn1.no
      s = src  - rnn.ni - rnn.no + rnn1.ni + rnn1.no
      @test rnn.w[dest, src] == rnn1.w[d, s]
    end

    #
    # RNN 2 Tests
    #

    # rnn2 input to output
    if dest > rnn.ni + rnn1.no && dest <= rnn.ni + rnn.no &&
      src   > rnn1.ni          && src  <= rnn.ni
      d = dest - rnn.ni - rnn1.no + rnn2.ni
      s = src  - rnn1.ni
      @test rnn.w[dest, src] == rnn2.w[d, s]
    end

    #= # rnn2 input to hidden =#
    if dest >  rnn.ni + rnn.no + rnn1.nh &&
       dest <= rnn.ni + rnn.no + rnn.nh  &&
       src  >  rnn1.ni &&
       src  <= rnn.ni
      d = dest - rnn.ni - rnn.no - rnn1.nh + rnn2.ni + rnn2.no
      s = src  - rnn1.ni
      @test rnn.w[dest, src] == rnn2.w[d, s]
    end

    # rnn2 output to output
    if dest > rnn.ni + rnn1.no && dest <= rnn.ni + rnn.no &&
       src  > rnn.ni + rnn1.no && src  <= rnn.ni + rnn.no
      d = dest - rnn.ni - rnn1.no + rnn2.ni
      s = src  - rnn.ni - rnn1.no + rnn2.ni
      @test rnn.w[dest, src] == rnn2.w[d, s]
    end

    # rnn2 output to hidden
    if dest > rnn.ni + rnn.no + rnn1.nh && dest <= rnn.ni + rnn.no + rnn.nh &&
       src  > rnn.ni + rnn1.no          && src  <= rnn.ni + rnn.no
      d = dest - rnn.ni - rnn.no - rnn1.nh + rnn2.ni + rnn2.no
      s = src  - rnn.ni - rnn1.no + rnn2.ni
      @test rnn.w[dest, src] == rnn2.w[d, s]
    end

    # rnn2 hidden to output
    if dest > rnn.ni + rnn1.no && dest <= rnn.ni + rnn.no &&
       src  > rnn.ni + rnn.no + rnn1.nh && src  <= rnn.ni + rnn.no + rnn.nh
      d = dest - rnn.ni - rnn1.no + rnn2.ni
      s = src  - rnn.ni - rnn.no - rnn1.nh + rnn2.ni + rnn2.no
      @test rnn.w[dest, src] == rnn2.w[d, s]
    end

    # rnn2 hidden to hidden
    if dest > rnn.ni + rnn.no + rnn1.nh && dest <= rnn.ni + rnn.no + rnn.nh &&
       src  > rnn.ni + rnn.no + rnn1.nh && src  <= rnn.ni + rnn.no + rnn.nh
      d = dest - rnn.ni - rnn.no - rnn1.nh + rnn2.ni + rnn2.no
      s = src  - rnn.ni - rnn.no - rnn1.nh + rnn2.ni + rnn2.no
      @test rnn.w[dest, src] == rnn2.w[d, s]
    end

  end
end

rnn   = RNN_t(5, 5, 10)
rnn.w = ones(size(rnn.w))
rnn.b = ones(size(rnn.b)) .* 2.0

rnn_add_input_neurons!(rnn,3)
@test rnn.ni == 8
@test rnn.no == 5
@test rnn.nh == 10
@test rnn.n  == 23
for i=1:23
  if i in [6,7,8]
    @test rnn.b[i] == 0.0
  else
    @test rnn.b[i] == 2.0
  end
end

for i=1:23
  for j=1:23
    if i in [6:8]
      @test rnn.w[i,j] == 0.0
    elseif j in [6:8]
      @test rnn.w[i,j] == 0.0
    else
      @test rnn.w[i,j] == 1.0
    end
  end
end

rnn_add_output_neurons!(rnn,3)
@test rnn.ni == 8
@test rnn.no == 8
@test rnn.nh == 10
@test rnn.n  == 26
for i=1:26
  if i in [6:8]
    @test rnn.b[i] == 0.0
  elseif i in [14:16]
    @test rnn.b[i] == 0.0
  else
    @test rnn.b[i] == 2.0
  end
end

for i=1:26
  for j=1:26
    if i in [6:8]
      @test rnn.w[i,j] == 0.0
    elseif i in [14:16]
      @test rnn.w[i,j] == 0.0
    elseif j in [6:8]
      @test rnn.w[i,j] == 0.0
    elseif j in [14:16]
      @test rnn.w[i,j] == 0.0
    else
      @test rnn.w[i,j] == 1.0
    end
  end
end

rnn_add_hidden_neurons!(rnn,5)
@test rnn.ni == 8
@test rnn.no == 8
@test rnn.nh == 15
@test rnn.n  == 31
@test size(rnn.w) == (31,31)
@test size(rnn.b) == (31,1)
@test size(rnn.o) == (31,1)
@test size(rnn.a) == (31,1)

for i=1:31
  if i in [6:8]
    @test rnn.b[i] == 0.0
  elseif i in [14:16]
    @test rnn.b[i] == 0.0
  elseif i in [27:31]
    @test rnn.b[i] == 0.0
  else
    @test rnn.b[i] == 2.0
  end
end

for i=1:26
  for j=1:26
    if i in [6:8]
      @test rnn.w[i,j] == 0.0
    elseif i in [14:16]
      @test rnn.w[i,j] == 0.0
    elseif j in [6:8]
      @test rnn.w[i,j] == 0.0
    elseif j in [14:16]
      @test rnn.w[i,j] == 0.0
    elseif i in [27:31]
      @test rnn.w[i,j] == 0.0
    elseif j in [27:31]
      @test rnn.w[i,j] == 0.0
    else
      @test rnn.w[i,j] == 1.0
    end
  end
end


text:
# ---
# title: 41. First Missing Positive
# id: problem41
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array
# link: <https://leetcode.com/problems/first-missing-positive/description/>
# hidden: true
# ---
#
# Given an unsorted integer array `nums`, find the smallest missing positive
# integer.
#
# **Follow up:**  Could you implement an algorithm that runs in `O(n)` time and
# uses constant extra space.?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,0]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,4,-1,1]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: nums = [7,8,9,11,12]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 300`
#   * `-231 <= nums[i] <= 231 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end