text:
using OpenSpecFun
using Base.Test

const SF = OpenSpecFun

# useful test functions for relative error, which differ from isapprox
# in that errc separately looks at the real and imaginay parts
err(z, x) = z == x ? 0.0 : abs(z - x) / abs(x)
errc(z, x) = max(err(real(z),real(x)), err(imag(z),imag(x)))
≅(a,b) = errc(a,b) ≤ 1e-13

# error functions
@testset "error functions" begin
    @test SF.erf(Float16(1)) ≈ 0.84270079294971486934
    @test SF.erf(1) ≈ 0.84270079294971486934
    @test SF.erfc(1) ≈ 0.15729920705028513066
    @test SF.erfc(Float16(1)) ≈ 0.15729920705028513066
    @test SF.erfcx(1) ≈ 0.42758357615580700442
    @test SF.erfcx(Float32(1)) ≈ 0.42758357615580700442
    @test SF.erfcx(Complex64(1)) ≈ 0.42758357615580700442
    @test SF.erfi(1) ≈ 1.6504257587975428760
    @test SF.erfinv(0.84270079294971486934) ≈ 1
    @test SF.erfcinv(0.15729920705028513066) ≈ 1
    @test SF.dawson(1) ≈ 0.53807950691276841914

    @test SF.erf(1+2im) ≈ -0.53664356577856503399-5.0491437034470346695im
    @test SF.erfc(1+2im) ≈ 1.5366435657785650340+5.0491437034470346695im
    @test SF.erfcx(1+2im) ≈ 0.14023958136627794370-0.22221344017989910261im
    @test SF.erfi(1+2im) ≈ -0.011259006028815025076+1.0036063427256517509im
    @test SF.dawson(1+2im) ≈ -13.388927316482919244-11.828715103889593303im

    for elty in [Float32,Float64]
        for x in logspace(-200, -0.01)
            @test_approx_eq_eps SF.erf(SF.erfinv(x)) x 1e-12*x
            @test_approx_eq_eps SF.erf(SF.erfinv(-x)) -x 1e-12*x
            @test_approx_eq_eps SF.erfc(SF.erfcinv(2*x)) 2*x 1e-12*x
            if x > 1e-20
                xf = Float32(x)
                @test_approx_eq_eps SF.erf(SF.erfinv(xf)) xf 1e-5*xf
                @test_approx_eq_eps SF.erf(SF.erfinv(-xf)) -xf 1e-5*xf
                @test_approx_eq_eps SF.erfc(SF.erfcinv(2xf)) 2xf 1e-5*xf
            end
        end
        @test SF.erfinv(one(elty)) == Inf
        @test SF.erfinv(-one(elty)) == -Inf
        @test_throws DomainError SF.erfinv(convert(elty,2.0))

        @test SF.erfcinv(zero(elty)) == Inf
        @test_throws DomainError SF.erfcinv(-one(elty))
    end

    @test SF.erfinv(one(Int)) == SF.erfinv(1.0)
    @test SF.erfcinv(one(Int)) == SF.erfcinv(1.0)
end

# airy
@testset "airy" begin
    @test SF.airy(1.8) ≈ SF.airyai(1.8)
    @test SF.airyprime(1.8) ≈ -0.0685247801186109345638
    @test SF.airyaiprime(1.8) ≈ SF.airyprime(1.8)
    @test SF.airybi(1.8) ≈ 2.595869356743906290060
    @test SF.airybiprime(1.8) ≈ 2.98554005084659907283
    @test_throws SF.AmosException SF.airy(200im)
    @test_throws SF.AmosException SF.airybi(200)
    @test_throws ArgumentError SF.airy(5,one(Complex128))
    z = 1.8 + 1.0im
    for elty in [Complex64,Complex128]
        @test SF.airy(convert(elty,1.8)) ≈ 0.0470362168668458052247
        z = convert(elty,z)
        @test SF.airyx(z) ≈ SF.airyx(0,z)
        @test SF.airyx(0, z) ≈ SF.airy(0, z) * exp(2/3 * z * sqrt(z))
        @test SF.airyx(1, z) ≈ SF.airy(1, z) * exp(2/3 * z * sqrt(z))
        @test SF.airyx(2, z) ≈ SF.airy(2, z) * exp(-abs(real(2/3 * z * sqrt(z))))
        @test SF.airyx(3, z) ≈ SF.airy(3, z) * exp(-abs(real(2/3 * z * sqrt(z))))
        @test_throws ArgumentError SF.airyx(5,z)
    end
    @test_throws MethodError SF.airy(complex(big(1.0)))
end

# bessely0, bessely1, besselj0, besselj1
@testset "bessel" begin
    @test SF.besselj0(Float32(2.0)) ≈ SF.besselj0(Float64(2.0))
    @test SF.besselj1(Float32(2.0)) ≈ SF.besselj1(Float64(2.0))
    @test SF.bessely0(Float32(2.0)) ≈ SF.bessely0(Float64(2.0))
    @test SF.bessely1(Float32(2.0)) ≈ SF.bessely1(Float64(2.0))
    @test SF.besselj0(2) ≈ SF.besselj0(2.0)
    @test SF.besselj1(2) ≈ SF.besselj1(2.0)
    @test SF.bessely0(2) ≈ SF.bessely0(2.0)
    @test SF.bessely1(2) ≈ SF.bessely1(2.0)
    @test SF.besselj0(2.0 + im) ≈ SF.besselj(0, 2.0 + im)
    @test SF.besselj1(2.0 + im) ≈ SF.besselj(1, 2.0 + im)
    @test SF.bessely0(2.0 + im) ≈ SF.bessely(0, 2.0 + im)
    @test SF.bessely1(2.0 + im) ≈ SF.bessely(1, 2.0 + im)

    @test_throws MethodError SF.besselj(1.2,big(1.0))
    @test_throws MethodError SF.besselj(1,complex(big(1.0)))
    @test_throws MethodError SF.besseljx(1,big(1.0))
    @test_throws MethodError SF.besseljx(1,complex(big(1.0)))

    # besselh
    true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
    @test SF.besselh(3,1,3) ≈ true_h133
    @test SF.besselh(-3,1,3) ≈ -true_h133
    @test SF.besselh(3,2,3) ≈ conj(true_h133)
    @test SF.besselh(-3,2,3) ≈ -conj(true_h133)
    @test_throws SF.AmosException SF.besselh(1,0)

    @test_throws MethodError SF.besselh(1,big(1.0))
    @test_throws MethodError SF.besselh(1,complex(big(1.0)))
    @test_throws MethodError SF.besselhx(1,big(1.0))
    @test_throws MethodError SF.besselhx(1,complex(big(1.0)))

    # besseli
    true_i33 = 0.95975362949600785698
    @test SF.besseli(3,3) ≈ true_i33
    @test SF.besseli(-3,3) ≈ true_i33
    @test SF.besseli(3,-3) ≈ -true_i33
    @test SF.besseli(-3,-3) ≈ -true_i33
    @test SF.besseli(Float32(-3),Complex64(-3,0)) ≈ -true_i33
    @test_throws SF.AmosException SF.besseli(1,1000)
    @test_throws DomainError SF.besseli(0.4,-1.0)

    @test_throws MethodError SF.besseli(1,big(1.0))
    @test_throws MethodError SF.besseli(1,complex(big(1.0)))
    @test_throws MethodError SF.besselix(1,big(1.0))
    @test_throws MethodError SF.besselix(1,complex(big(1.0)))

    # besselj
    @test SF.besselj(0,0) == 1
    for i = 1:5
        @test SF.besselj(i,0) == 0
        @test SF.besselj(-i,0) == 0
        @test SF.besselj(-i,Float32(0)) == 0
        @test SF.besselj(-i,Float32(0)) == 0
    end

    j33 = SF.besselj(3,3.)
    @test SF.besselj(3,3) == j33
    @test SF.besselj(-3,-3) == j33
    @test SF.besselj(-3,3) == -j33
    @test SF.besselj(3,-3) == -j33
    @test SF.besselj(3,3f0) ≈ j33
    @test SF.besselj(3,complex(3.)) ≈ j33
    @test SF.besselj(3,complex(3f0)) ≈ j33
    @test SF.besselj(3,complex(3)) ≈ j33

    j43 = SF.besselj(4,3.)
    @test SF.besselj(4,3) == j43
    @test SF.besselj(-4,-3) == j43
    @test SF.besselj(-4,3) == j43
    @test SF.besselj(4,-3) == j43
    @test SF.besselj(4,3f0) ≈ j43
    @test SF.besselj(4,complex(3.)) ≈ j43
    @test SF.besselj(4,complex(3f0)) ≈ j43
    @test SF.besselj(4,complex(3)) ≈ j43

    @test j33 ≈ 0.30906272225525164362
    @test j43 ≈ 0.13203418392461221033
    @test_throws DomainError    SF.besselj(0.1, -0.4)
    @test SF.besselj(0.1, complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im
    @test SF.besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im
    @test SF.besselj(1, 3im) ≈ 3.953370217402609396im
    @test SF.besselj(1.0,3im) ≈ SF.besselj(1,3im)
    @test_throws SF.AmosException SF.besselj(20,1000im)
    @test_throws MethodError SF.besselj(big(1.0),3im)

    # besselk
    true_k33 = 0.12217037575718356792
    @test SF.besselk(3,3) ≈ true_k33
    @test SF.besselk(-3,3) ≈ true_k33
    true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
    @test_throws DomainError SF.besselk(3,-3)
    @test SF.besselk(3,complex(-3)) ≈ true_k3m3
    @test SF.besselk(-3,complex(-3)) ≈ true_k3m3
    @test_throws SF.AmosException SF.besselk(200,0.01)
    # base julia issue #6564
    @test SF.besselk(1.0,0.0) == Inf

    @test_throws MethodError SF.besselk(1,big(1.0))
    @test_throws MethodError SF.besselk(1,complex(big(1.0)))
    @test_throws MethodError SF.besselkx(1,big(1.0))
    @test_throws MethodError SF.besselkx(1,complex(big(1.0)))

    # bessely
    y33 = SF.bessely(3,3.)
    @test SF.bessely(3,3) == y33
    @test SF.bessely(3.,3.) == y33
    @test SF.bessely(3,Float32(3.)) ≈ y33
    @test SF.bessely(-3,3) ≈ -y33
    @test y33 ≈ -0.53854161610503161800
    @test_throws DomainError SF.bessely(3,-3)
    @test SF.bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im
    @test_throws SF.AmosException SF.bessely(200.5,0.1)
    @test_throws DomainError SF.bessely(0.4,-1.0)
    @test_throws DomainError SF.bessely(0.4,Float32(-1.0))
    @test_throws DomainError SF.bessely(1,Float32(-1.0))
    @test_throws DomainError SF.bessely(Cint(3),Float32(-3.))
    @test_throws DomainError SF.bessely(Cint(3),Float64(-3.))

    @test_throws MethodError SF.bessely(1.2,big(1.0))
    @test_throws MethodError SF.bessely(1,complex(big(1.0)))
    @test_throws MethodError SF.besselyx(1,big(1.0))
    @test_throws MethodError SF.besselyx(1,complex(big(1.0)))

    # besselhx
    for elty in [Complex64,Complex128]
        z = convert(elty, 1.0 + 1.9im)
        @test SF.besselhx(1.0, 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
        @test SF.besselhx(Float32(1.0), 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
    end

    @test_throws MethodError SF.besselh(1,1,big(1.0))
    @test_throws MethodError SF.besselh(1,1,complex(big(1.0)))
    @test_throws MethodError SF.besselhx(1,1,big(1.0))
    @test_throws MethodError SF.besselhx(1,1,complex(big(1.0)))

    # base julia issue #6653
    for f in (SF.besselj,SF.bessely,SF.besseli,SF.besselk,SF.hankelh1,SF.hankelh2)
        @test f(0,1) ≈ f(0,Complex128(1))
        @test f(0,1) ≈ f(0,Complex64(1))
    end

    # scaled bessel[ijky] and hankelh[12]
    for x in (1.0, 0.0, -1.0), y in (1.0, 0.0, -1.0), nu in (1.0, 0.0, -1.0)
        z = Complex128(x + y * im)
        z == zero(z) || @test SF.hankelh1x(nu, z) ≈ SF.hankelh1(nu, z) * exp(-z * im)
        z == zero(z) || @test SF.hankelh2x(nu, z) ≈ SF.hankelh2(nu, z) * exp(z * im)
        (nu < 0 && z == zero(z)) || @test SF.besselix(nu, z) ≈ SF.besseli(nu, z) * exp(-abs(real(z)))
        (nu < 0 && z == zero(z)) || @test SF.besseljx(nu, z) ≈ SF.besselj(nu, z) * exp(-abs(imag(z)))
        z == zero(z) || @test SF.besselkx(nu, z) ≈ SF.besselk(nu, z) * exp(z)
        z == zero(z) || @test SF.besselyx(nu, z) ≈ SF.bessely(nu, z) * exp(-abs(imag(z)))
    end
    @test_throws SF.AmosException SF.hankelh1x(1, 0)
    @test_throws SF.AmosException SF.hankelh2x(1, 0)
    @test_throws SF.AmosException SF.besselix(-1, 0)
    @test_throws SF.AmosException SF.besseljx(-1, 0)
    @test SF.besselkx(1, 0) == Inf
    @test_throws SF.AmosException SF.besselyx(1, 0)
    @test_throws DomainError SF.besselix(0.4,-1.0)
    @test_throws DomainError SF.besseljx(0.4, -1.0)
    @test_throws DomainError SF.besselkx(0.4,-1.0)
    @test_throws DomainError SF.besselyx(0.4,-1.0)
end

# SF.beta, SF.lbeta
@testset "beta, lbeta" begin
    @test SF.beta(3/2,7/2) ≈ 5π/128
    @test SF.beta(3,5) ≈ 1/105
    @test SF.lbeta(5,4) ≈ log(SF.beta(5,4))
    @test SF.beta(5,4) ≈ SF.beta(4,5)
    @test SF.beta(-1/2, 3) ≈ SF.beta(-1/2 + 0im, 3 + 0im) ≈ -16/3
    @test SF.lbeta(-1/2, 3) ≈ log(16/3)
    @test SF.beta(Float32(5),Float32(4)) == SF.beta(Float32(4),Float32(5))
    @test SF.beta(3,5) ≈ SF.beta(3+0im,5+0im)
    @test(SF.beta(3.2+0.1im,5.3+0.3im) ≈ exp(SF.lbeta(3.2+0.1im,5.3+0.3im)) ≈
          0.00634645247782269506319336871208405439180447035257028310080 -
          0.00169495384841964531409376316336552555952269360134349446910im)
end

# SF.gamma, SF.lgamma (complex argument)
@testset "gamma, lgamma" begin
    if Base.Math.libm == "libopenlibm"
        @test SF.gamma(Float64[1:25;]) == SF.gamma(1:25)
    else
        @test SF.gamma(Float64[1:25;]) ≈ SF.gamma(1:25)
    end
    for elty in (Float32, Float64)
        @test SF.gamma(convert(elty,1/2)) ≈ convert(elty,sqrt(π))
        @test SF.gamma(convert(elty,-1/2)) ≈ convert(elty,-2sqrt(π))
        @test SF.lgamma(convert(elty,-1/2)) ≈ convert(elty,log(abs(SF.gamma(-1/2))))
    end
    @test SF.lgamma(1.4+3.7im) ≈ -3.7094025330996841898 + 2.4568090502768651184im
    @test SF.lgamma(1.4+3.7im) ≈ log(SF.gamma(1.4+3.7im))
    @test SF.lgamma(-4.2+0im) ≈ SF.lgamma(-4.2)-pi*im
    @test factorial(3.0) == SF.gamma(4.0) == factorial(3)
    for x in (3.2, 2+1im, 3//2, 3.2+0.1im)
        @test factorial(x) == SF.gamma(1+x)
    end
    @test lfact(1) == 0
    @test lfact(2) == SF.lgamma(3)
end

# SF.digamma
@testset "digamma" begin
    for elty in (Float32, Float64)
        @test SF.digamma(convert(elty, 9)) ≈ convert(elty, 2.140641477955609996536345)
        @test SF.digamma(convert(elty, 2.5)) ≈ convert(elty, 0.7031566406452431872257)
        @test SF.digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041107679516822)
        @test SF.digamma(convert(elty, 7e-4)) ≈ convert(elty, -1429.147493371120205005198)
        @test SF.digamma(convert(elty, 7e-5)) ≈ convert(elty, -14286.29138623969227538398)
        @test SF.digamma(convert(elty, 7e-6)) ≈ convert(elty, -142857.7200612932791081972)
        @test SF.digamma(convert(elty, 2e-6)) ≈ convert(elty, -500000.5772123750382073831)
        @test SF.digamma(convert(elty, 1e-6)) ≈ convert(elty, -1000000.577214019968668068)
        @test SF.digamma(convert(elty, 7e-7)) ≈ convert(elty, -1428572.005785942019703646)
        @test SF.digamma(convert(elty, -0.5)) ≈ convert(elty, .03648997397857652055902367)
        @test SF.digamma(convert(elty, -1.1)) ≈ convert(elty,  10.15416395914385769902271)

        @test SF.digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041108)
        @test SF.digamma(convert(elty, 1/2)) ≈ convert(elty, -γ - log(4))
        @test SF.digamma(convert(elty, 1)) ≈ convert(elty, -γ)
        @test SF.digamma(convert(elty, 2)) ≈ convert(elty, 1 - γ)
        @test SF.digamma(convert(elty, 3)) ≈ convert(elty, 3/2 - γ)
        @test SF.digamma(convert(elty, 4)) ≈ convert(elty, 11/6 - γ)
        @test SF.digamma(convert(elty, 5)) ≈ convert(elty, 25/12 - γ)
        @test SF.digamma(convert(elty, 10)) ≈ convert(elty, 7129/2520 - γ)
    end

    for x in -10.2:0.3456:50
        @test 1e-12 > err(SF.digamma(x+0im), SF.digamma(x))
    end

    # Test against Wolfram Alpha
    @test SF.digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232
    @test SF.digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im
    @test SF.digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im

    @test isa([SF.digamma(x) for x in [1.0]], Vector{Float64})
end

# SF.trigamma
@testset "trigamma" begin
    for elty in (Float32, Float64)
        @test SF.trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
        @test SF.trigamma(convert(elty, 1/2)) ≈ convert(elty, π^2/2)
        @test SF.trigamma(convert(elty, 1)) ≈ convert(elty, π^2/6)
        @test SF.trigamma(convert(elty, 2)) ≈ convert(elty, π^2/6 - 1)
        @test SF.trigamma(convert(elty, 3)) ≈ convert(elty, π^2/6 - 5/4)
        @test SF.trigamma(convert(elty, 4)) ≈ convert(elty, π^2/6 - 49/36)
        @test SF.trigamma(convert(elty, 5)) ≈ convert(elty, π^2/6 - 205/144)
        @test SF.trigamma(convert(elty, 10)) ≈ convert(elty, π^2/6 - 9778141/6350400)
    end

    # Test against Wolfram Alpha
    @test SF.trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546
    @test SF.trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im
    @test SF.trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im

    @test isa([SF.trigamma(x) for x in [1.0]], Vector{Float64})
end

# SF.invdigamma
@testset "invdigamma" begin
    for elty in (Float32, Float64)
        for val in [0.001, 0.01, 0.1, 1.0, 10.0]
            @test abs(SF.invdigamma(SF.digamma(convert(elty, val))) - convert(elty, val)) < 1e-8
        end
    end
    @test abs(SF.invdigamma(2)) == abs(SF.invdigamma(2.))
end

@testset "polygamma" begin
    @test SF.polygamma(20, 7.) ≈ -4.644616027240543262561198814998587152547
    @test SF.polygamma(20, Float16(7.)) ≈ -4.644616027240543262561198814998587152547

    # Test against Wolfram Alpha
    @test SF.polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600
    @test SF.polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im
    @test SF.polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im
    @test 1e-11 > err(SF.polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im)
    @test SF.polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im
    @test 1e-11 > err(SF.polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291)

    @test SF.polygamma(4, -3.1+Inf*im) == SF.polygamma(4, 3.1+Inf*im) == 0
    @test SF.polygamma(4, -0.0) == Inf == -SF.polygamma(4, +0.0)
    @test isa([SF.polygamma(3,x) for x in [1.0]], Vector{Float64})
    @test SF.polygamma(3,5) ≈ SF.polygamma(3,5.)
end

# SF.eta, SF.zeta
@testset "eta, zeta" begin
    @test SF.eta(1) ≈ log(2)
    @test SF.eta(2) ≈ pi^2/12
    @test SF.eta(Float32(2)) ≈ SF.eta(2)
    @test SF.eta(Complex64(2)) ≈ SF.eta(2)
    @test SF.zeta(0) ≈ -0.5
    @test SF.zeta(2) ≈ pi^2/6
    @test SF.zeta(Complex64(2)) ≈ SF.zeta(2)
    @test SF.zeta(4) ≈ pi^4/90
    @test SF.zeta(1,Float16(2.)) ≈ SF.zeta(1,2.)
    @test SF.zeta(1.,Float16(2.)) ≈ SF.zeta(1,2.)
    @test SF.zeta(Float16(1.),Float16(2.)) ≈ SF.zeta(1,2.)
    @test isnan(SF.zeta(NaN))
    @test isnan(SF.zeta(1.0e0))
    @test isnan(SF.zeta(1.0f0))
    @test isnan(SF.zeta(complex(0,Inf)))
    @test isnan(SF.zeta(complex(-Inf,0)))

    @test SF.zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im
    @test SF.zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im
    @test SF.zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im
    @test 1e-12 > errc(SF.zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im)
    @test SF.zeta(3.1,-4.2) ≅ SF.zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295
    @test 1e-15 > errc(SF.zeta(3.1+0im,-4.2), SF.zeta(3.1,-4.2+0im))
    @test SF.zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234
    @test SF.zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14
    @test SF.zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9
    @test 1e-12 > err(SF.zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14)

    @test SF.zeta(4, +0.0) == SF.zeta(4, -0.0) ≅ pi^4 / 90
    @test SF.zeta(5, +0.0) == SF.zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974
    @test SF.zeta(Inf, 1.) == 1
    @test SF.zeta(Inf, 2.) == 0
    @test isnan(SF.zeta(NaN, 1.))

    @test SF.zeta(2 + 1im, -1.1) ≅ SF.zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im

    @test SF.zeta(-3.0, 7.0) ≅ -52919/120
    @test SF.zeta(-3.0, -7.0) ≅ 94081/120
    @test SF.zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719
    @test SF.zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576
    @test SF.zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748
    @test SF.zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161
    @test SF.zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im
    @test SF.zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im
    @test SF.zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im
    @test SF.zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im
    @test SF.zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im
    @test SF.zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im
    @test SF.zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im
    @test SF.zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6
    @test SF.zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078
    @test SF.zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977
    @test SF.zeta(-3.2+0.1im, 0.0) ≅ SF.zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im
    @test SF.zeta(-3.2, 0.0) ≅ SF.zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378

    @test 1e-14 > err(SF.eta(1+1e-9), 0.693147180719814213126976796937244130533478392539154928250926)
    @test 1e-14 > err(SF.eta(1+5e-3), 0.693945708117842473436705502427198307157819636785324430166786)
    @test 1e-13 > err(SF.eta(1+7.1e-3), 0.694280602623782381522315484518617968911346216413679911124758)
    @test 1e-13 > err(SF.eta(1+8.1e-3), 0.694439974969407464789106040237272613286958025383030083792151)
    @test 1e-13 > err(SF.eta(1 - 2.1e-3 + 2e-3 * im), 0.69281144248566007063525513903467244218447562492555491581+0.00032001240133205689782368277733081683574922990400416791019im)
    @test 1e-13 > err(SF.eta(1 + 5e-3 + 5e-3 * im), 0.69394652468453741050544512825906295778565788963009705146+0.00079771059614865948716292388790427833787298296229354721960im)
    @test 1e-12 > errc(SF.zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im)
    @test 1e-13 > errc(SF.zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im)

    # Base Julia issue #7169: (TODO: better accuracy should be possible?)
    @test 1e-9 > errc(SF.zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
    @test 1e-12 > errc(SF.zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im)
    @test 1e-9 > errc(SF.zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
    @test 1e-13 > errc(SF.zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im)
end


text:
import Dagger.Sch: SchedulerOptions, ThunkOptions

function inc(x)
    x+1
end
@everywhere begin
function checkwid(x...)
    @assert myid() == 1
    return 1
end
end

@testset "Scheduler" begin
    @testset "Scheduler options: single worker" begin
        a = delayed(checkwid)(1)
        b = delayed(checkwid)(2)
        c = delayed(checkwid)(a,b)

        @test collect(Context(), c; options=SchedulerOptions(1)) == 1
    end
    @testset "Thunk options: single worker" begin
        options = ThunkOptions(1)
        a = delayed(checkwid; options=options)(1)

        @test collect(Context(), a) == 1
    end
end

ext: jl


text:
using InvariantCausalPrediction
using Test

@testset "InvariantCausalPrediction.jl" begin
    # Write your tests here.
end


text:
using ONNX
using Test

include("readwrite.jl")
include("ort.jl")
include("saveload.jl")


text:
"""
    precedence_graph(jobs)
Construct a SimpleDiGraph with as vertices job-ids and an edge between vertice v_i and v_j if a precedence i ≺ j exists.
"""
function precedence_graph(jobs)
    location_dict = group(loc, jobs)
    foreach( k -> sort!(location_dict[k]) , keys(location_dict))
    pg = SimpleDiGraph(length(jobs))
    for pjobs in values(location_dict)
        iter = Iterators.zip(pjobs[1:end-1],pjobs[2:end])
        for (s,d) in iter
            add_edge!(pg, id(s), id(d))
        end
    end
    return pg
end

"""
    precedence_matrix(jobs)
Construct a Matrix{Bool} containing entries for each pair of jobs, with 'true' entries signifying the presence of a precedence relation.
"""
precedence_matrix(jobs) = map(J -> ≺(J...), Iterators.product(jobs,jobs))

"""
    assignment_matrix(jobs, cranes)
Construct a Matrix{Bool} containing entries for each job-crane-pair, with 'true' entries signfifying the possibility of an assignment of the given job to the given crane.
"""
assignment_matrix(jobs, cranes) = map( x -> loc(x[1]) ∈ zone(x[2]), Iterators.product(jobs, cranes))

"""
    compare(jobs, property)
Perform a pairwise comparison of `property` on the `jobs`.
"""
compare(jobs, property) = map(J -> property(J...), Iterators.product(jobs, jobs))

function instance_to_lb_solution(instance)
    Ω = jobs(instance)
    Q = cranes(instance)
    L = jobs_by_loc(Ω)
    jd = _jobdata(L)
    qd = _cranedata(Q)
    @show jd
    @show qd
    return jd, qd
end

"""
    jobs_by_loc(jobs)

Returns a vector of vectors of jobs, with on each index `i` in the top vector the _ordered_ list of jobs on location `i`
"""
function jobs_by_loc(jobs)
    pg = precedence_graph(jobs)
    L = [Job[] for _ in 1:length(jobs)]
    for v in topological_sort_by_dfs(pg)
        j = jobs[v]
        push!(L[loc(j)], j)
    end
    return L
end

"""
    loctojobdict
"""
function loctojobdict(jobs)
    L = jobs_by_loc(jobs)
    D = Dict{Int, Vector{Job}}()
    for (l,jobs) in enumerate(L)
        if !isempty(jobs)
            D[l] = jobs
        end
    end
    return D
end

function _jobdata(L)
    jd = Tuple{Int,Int}[]
    for l in L
        t = 0
        for j in l
            ts = max(t, t_arrival(j))
            push!(jd, (ts,id(j)))
            t = ts + t_processing(j)
        end
    end
    return jd
end

_cranedata(Q) = map(q -> (0,starting_position(q)), Q)

function lb_obj(jobs, jd)
    cmax = maximum(j -> first(j) + t_processing(jobs[last(j)]), jd)
    twt = 0
    for (ts, jid) in jd
        j = jobs[jid]
        if istruck(j)
            twt += t_processing(j) + ts - t_arrival(j)
        end
    end
    return cmax, twt
end

"""
    assignment_ranges(Ω, Q)

Return a vector of ranges designating which craens can be assigned to each job.
"""
function assignment_ranges(Ω, Q)
    R = UnitRange{Int}[]
    for j in Ω
        cmin = findfirst( c -> loc(j) ∈ zone(c), Q)
        cmax = findlast( c -> loc(j) ∈ zone(c), Q)
        push!(R, cmin:cmax)
    end
    return R
end


text:
using Distributions, Random
using Test, LinearAlgebra, PDMats

p = 4
n1 = p + abs(10randn())
n2 = p + abs(10randn())

B = MatrixBeta(p, n1, n2)
C = MatrixBeta(1, n1, n2)
c = Beta(n1/2, n2/2)

Z = rand(C)
z = Z[1]

@testset "MatrixBeta promotion during construction" begin
    R = MatrixBeta(p, Float32(n1), BigFloat(n2))
    @test partype(R) == BigFloat
end

@testset "MatrixBeta construction errors" begin
    @test_throws ArgumentError MatrixBeta(-2, n1, n2)
    @test_throws ErrorException MatrixBeta(p, p - 2, n2)
    @test_throws ErrorException MatrixBeta(p, n1, p - 2)
end

@testset "MatrixBeta params" begin
    df1, df2 = params(B)
    @test df1 == n1
    @test df2 == n2
end

@testset "MatrixBeta dim" begin
    @test dim(B) == p
    @test dim(C) == 1
end

@testset "MatrixBeta size" begin
    @test size(B) == (p, p)
    @test size(C) == (1, 1)
end

@testset "MatrixBeta rank" begin
    @test rank(B) == p
    @test rank(C) == 1
    @test rank(B) == rank(rand(B))
    @test rank(C) == rank(rand(C))
end

@testset "MatrixBeta insupport" begin
    @test insupport(B, rand(B))
    @test insupport(C, rand(C))
    @test all(0 .< eigvals(rand(B)) .< 1)

    @test !insupport(B, rand(B) + rand(B) * im)  #  not real
    @test !insupport(B, rand(C))                 #  wrong dims
    @test !insupport(C, rand(B))                 #  wrong dims
    @test !insupport(B, rand(p, p))              #  U isn't posdef
    @test !insupport(B, rand(B) + I)             #  I - U isn't posdef
end

@testset "MatrixBeta logpdf" begin
    @test logpdf(C, Z) ≈ logpdf(c, z)
end

@testset "MatrixBeta cov and var" begin
    @test vec(var(B)) ≈ diag(cov(B))
    @test reshape(cov(B), p, p, p, p) ≈ cov(B, Val(false))
    @test cov(C)[1] ≈ var(c)
end

@testset "MatrixBeta sample moments" begin
    @test isapprox(mean(rand(B, 10000)), mean(B) , atol = 0.1)
    @test isapprox(mean(rand(C, 10000))[1, 1], mean(c) , atol = 0.1)

    @test isapprox(cov(hcat(vec.(rand(B, 100000))...)'), cov(B) , atol = 0.1)
    @test isapprox(cov(hcat(vec.(rand(C, 100000))...)'), cov(C) , atol = 0.1)
end

@testset "MatrixBeta conversion" for elty in (Float32, Float64, BigFloat)
    Bel1 = convert(MatrixBeta{elty}, B)
    Bel2 = convert(MatrixBeta{elty}, B.W1, B.W2, B.logc0)

    @test partype(Bel1) == elty
    @test partype(Bel2) == elty
end


text:
using BinDeps, UrlDownload

export download_datasets

# tuples of (dataset_name, filename, url)
# not a const to avoid errors when updating urls
function get_dataset_info()
    [
    # https://globalwindatlas.info/api/gis/global/wind-speed/10
    # https://globalwindatlas.info/api/gis/global/wind-speed/50
    # https://globalwindatlas.info/api/gis/global/wind-speed/100
    # https://globalwindatlas.info/api/gis/global/wind-speed/150
    # https://globalwindatlas.info/api/gis/global/wind-speed/200
    "GWA100"       ("Global Wind Atlas", "Global Wind Atlas v3 - 100m wind speed.tif",
        "https://chalmersuniversity.box.com/shared/static/wfr6dm9bcmj0mcqtdn0uimhg0otd4ht1.tif")
    "GWA150"       ("Global Wind Atlas", "Global Wind Atlas v3 - 150m wind speed.tif",
        "https://chalmersuniversity.box.com/shared/static/ghexnwa7crukl58nkwric8v9kmlymvj2.tif")
    "GWA200"       ("Global Wind Atlas", "Global Wind Atlas v3 - 200m wind speed.tif",
        "https://chalmersuniversity.box.com/shared/static/7ib2jdni6hu9uwe1hp1mqoeeyqtg5601.tif")
    "WDPA"         ("WDPA (protected areas):", "WDPA.zip", url_WDPA("Aug2021"))
        # "https://chalmersuniversity.box.com/shared/static/wn1kznvy7qh1issqcxdlsq64kgtkaayi.zip")
    "GADM"         ("GADM (global administrative areas)", "gadm36.zip",
        "https://biogeo.ucdavis.edu/data/gadm3.6/gadm36_shp.zip")
    "NUTS"         ("NUTS (administrative areas in Europe)", "nuts-2016-01m.shp.zip",
        "https://ec.europa.eu/eurostat/cache/GISCO/distribution/v2/nuts/download/ref-nuts-2016-01m.shp.zip")
    "landcover"    ("Land Cover", "Landcover - USGS MODIS.tif",
        "https://chalmersuniversity.box.com/shared/static/zm8zdkv1mz0wns0u77afkna95wbo7ggo.tif")
    "topography"   ("ETOPO1 Topography", "ETOPO1_Ice_c_geotiff.zip",
        "https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/cell_registered/georeferenced_tiff/ETOPO1_Ice_c_geotiff.zip")
    "ssppop1"      ("1-km downscaled spatial population scenarios (Gao et al), file 1/3", "temp_ssp1.tar",
        "http://www.cgd.ucar.edu/iam/modeling/data/SSP1_1km_netcdf.tar")
    "ssppop2"      ("1-km downscaled spatial population scenarios (Gao et al), file 2/3", "temp_ssp2.tar",
        "http://www.cgd.ucar.edu/iam/modeling/data/SSP2_1km_netcdf.tar")
    "ssppop3"      ("1-km downscaled spatial population scenarios (Gao et al), file 3/3", "temp_ssp3.tar",
        "http://www.cgd.ucar.edu/iam/modeling/data/SSP3_1km_netcdf.tar")
    "gdppop"       ("Global population & GDP", "temp_popgdp.zip",
        "https://osf.io/7jv3n/download")
    "powerplants"  ("WRI Global Power Plant Database", "WRI - Global Power Plant Database v1.10.zip",
        "https://chalmersuniversity.box.com/shared/static/ss6gycw7hf10e1fiicbxl5rgk08q5xr9.zip")
        # switch to the official v1.2 link later (some plant cleanup is hardcoded for v1.1)
        # "http://datasets.wri.org/dataset/540dcf46-f287-47ac-985d-269b04bea4c6/resource/c240ed2e-1190-4d7e-b1da-c66b72e08858/download/globalpowerplantdatabasev120")
    "timezones"    ("Time zone shape file", "timezones-with-oceans.shapefile.zip",
        "https://github.com/evansiroky/timezone-boundary-builder/releases/download/2019b/timezones-with-oceans.shapefile.zip")
    "monthlywind"  ("Average monthly wind speeds 1979-2019", "era5monthlywind.h5",
        "https://chalmersuniversity.box.com/shared/static/otvb5nq0lz5ntqx0e65kocos2afo7gas.h5")
    "monthlysolar" ("Average monthly solar insolation 1979-2019", "era5monthlysolar.h5",
        "https://chalmersuniversity.box.com/shared/static/jlpspmp9ou96hk7xno46rf79wg3d5hqc.h5")
    "various"      ("Various smaller datasets", "Various_smaller_datasets.zip",
        "https://chalmersuniversity.box.com/shared/static/w3pmx4xhgorgd6jejv23gn4ycsnza8s6.zip")
]
end
download_datasets(startfrom::Int) = download_datasets(get_dataset_info()[startfrom:end, 1]...)

url_WDPA(monthyear=Dates.format(now(), "uyyyy")) = "https://d1gam3xoknrgr2.cloudfront.net/current/WDPA_$(monthyear)_Public_shp.zip"

function download_datasets(shortnames::String...)
    datafolder = getconfig("datafolder")
    if !isfile(datafolder)
        mkpath(datafolder)
    end

    dataset_info = get_dataset_info()
    datasets = Dict(r[1] => r[2] for r in eachrow(dataset_info))
    shortnames = isempty(shortnames) ? dataset_info[:, 1] : shortnames

    for (i, shortname) in enumerate(shortnames)
        fullname, filename, url = datasets[shortname]

        println("\nDownloading dataset $i: $fullname")
        try
            download_progressbar(url, joinpath(datafolder, filename))
        catch
            guess_url = url_WDPA()
            if shortname == "WDPA" && url != guess_url
                println("\nHardcoded WDPA url no longer working, probably due to an end-of-month update at www.protectedplanet.net.")
                println("Retrying using a guessed url for the current month...")
                download_progressbar(guess_url, joinpath(datafolder, filename))
            else
                rethrow()
            end
        end
        unpack_and_cleanup(shortname, filename, datafolder, dataset_info)
    end
    println("\nDownloads complete.")
end

function unpack_and_cleanup(shortname, filename, datafolder, dataset_info)
    WDPA_filename = splitext(basename(dataset_info[4, 2][3]))[1]

    foldername, extension = splitext(filename)
    fullpath = joinpath(datafolder, filename)

    if extension in [".zip", ".tar"] && isfile(fullpath)
        println("\nUnpacking archive: $filename")
        unpack(fullpath, joinpath(datafolder, foldername), extension)
    end

    function renameWDPAfiles(WDPAfolder)
        for filename in readdir(WDPAfolder)
            newname = replace(filename, WDPA_filename => "WDPA-shapefile")
            if newname != filename
                mv(joinpath(WDPAfolder, filename), joinpath(WDPAfolder, newname))
            end
        end
    end

    if shortname == "WDPA"
        renameWDPAfiles(joinpath(datafolder, "WDPA"))
        for i = 0:2
            foldername = "WDPA_Oct2021_Public_shp_$i"
            println("\nUnpacking archive: $foldername.zip")
            unpack(joinpath(datafolder, "WDPA", "$foldername.zip"),
                    joinpath(datafolder, "WDPA", foldername), ".zip")
            renameWDPAfiles(joinpath(datafolder, "WDPA", foldername))
        end
    elseif shortname == "NUTS"
        filename = "NUTS_RG_01M_2016_4326_LEVL_3.shp.zip"
        println("\nUnpacking archive: $filename")
        unpack(joinpath(datafolder, "nuts-2016-01m.shp", filename), joinpath(datafolder, "nuts2016-level3"), ".zip")
        rm(joinpath(datafolder, "nuts-2016-01m.shp"), force=true, recursive=true)
    elseif shortname == "topography"
        mv(joinpath(datafolder, "ETOPO1_Ice_c_geotiff", "ETOPO1_Ice_c_geotiff.tif"), joinpath(datafolder, "ETOPO1_Ice_c_geotiff.tif"))
        rm(joinpath(datafolder, "ETOPO1_Ice_c_geotiff"))
    elseif shortname[1:end-1] == "ssppop"
        sspfolder = joinpath(datafolder, "SSP_1km")
        !isdir(sspfolder) && mkdir(sspfolder)
        for ssp = 1:3
            for y = 2010:10:2100
                sspfile = joinpath(datafolder, "temp_ssp$ssp", "SSP$(ssp)_1km", "ssp$(ssp)_total_$y.nc4")
                isfile(sspfile) && mv(sspfile, joinpath(sspfolder, "ssp$(ssp)_total_$y.nc4"))
            end
            rm(joinpath(datafolder, "temp_ssp$ssp"), force=true, recursive=true)
        end
    elseif shortname == "gdppop"
        mv(joinpath(datafolder, "temp_popgdp", "data"), joinpath(datafolder, "global_population_and_gdp"))
        rm(joinpath(datafolder, "temp_popgdp"))
    elseif shortname == "powerplants"
        mv(joinpath(datafolder, "WRI - Global Power Plant Database v1.10"), joinpath(datafolder, "tempWRI"))
        mv(joinpath(datafolder, "tempWRI", "WRI - Global Power Plant Database v1.10"),
            joinpath(datafolder, "WRI - Global Power Plant Database v1.10"))
        rm(joinpath(datafolder, "tempWRI"))
    elseif shortname == "various"
        mixeddir = joinpath(datafolder, "Various_smaller_datasets")
        for file in readdir(mixeddir)
            mv(joinpath(mixeddir, file), joinpath(datafolder, file), force=true)
        end
        rm(mixeddir)
    end

    if extension in [".zip", ".tar"] && isfile(fullpath)
        rm(fullpath)
    end
end

function download_progressbar(url::AbstractString, filename::AbstractString)
    println("Downloading to $filename...")
    # UrlDownload can take care of unpacking too, look into this.
    urldownload(url, true, compress=:none, parser=identity, save_raw=filename)
end

function unpack(inputfilename, outputpath, extension)
    !isdir(outputpath) && mkdir(outputpath)
    run(BinDeps.unpack_cmd(inputfilename, outputpath, extension, ""))
end

max_stars_repo_name: janohlen/GlobalEnergyGIS

max_issues_repo_name: janohlen/GlobalEnergyGIS

max_forks_repo_name: janohlen/GlobalEnergyGIS


text:
using StochasticDelayDiffEq
using Test

@testset "discrete callback (#34)" begin
    f(u, h, p, t) = 0.3 * h(p, t - p)
    g(u, h, p, t) = 0.1
    h(p, t) = 0.0
    prob = SDDEProblem(f, g, h, (0.0, 4.0), 0.5; constant_lags=(0.5,))

    # event at `t = 3`
    cb = DiscreteCallback(
        (u, t, integrator) -> t == 3,
        integrator -> (integrator.u = -integrator.u)
    )

    sol = solve(prob, RKMil(), tstops=(3,), callback=cb)
    ts = findall(x -> x == 3, sol.t)
    @test length(ts) == 2
    @test sol.u[ts[1]] == -sol.u[ts[2]]
    @test sol(3.0; continuity=:right) == -sol(3.0; continuity=:left)
end


text:
module BitwiseKoans

function not(x)
    ~x
end

function and(x, y)
    x & y
end

function or(x, y)
    x | y
end

function xor_koan(x, y)
    xor(x, y)
end

function right_shift(x, y)
    x >> y
end

function left_shift(x, y)
    x << y
end

end


text:
using DataFrames, CSV, Plots, StatsPlots

ProjDir = @__DIR__

nts = [1, 2, 4, 8]
N = 6

arm_no_tbb = CSV.read(joinpath(ProjDir, "results", "arm_no_tbb.csv"), DataFrame)
intel_tbb = CSV.read(joinpath(ProjDir, "results", "intel_tbb.csv"), DataFrame)

#arm_no_tbb |> display
#intel_tbb |> display

# Sort on (:num_threads,:num_chains0 instead on (:num_threads, :num_cpp_chains)
#sort(intel_tbb, [:num_threads, :num_chains])

function timings_plot(arm::DataFrame, intel::DataFrame; ts = [1, 8])
    ptype = ["arm", "intel"]
    colors = [:darkred, :darkblue, :darkgreen, :black]
    fig = plot(; xlim=(0, 9), ylim=(0, 100),
        xlab="chains (Julia or c++)", ylab="elapsed time [s]")
    for (indx, df) in enumerate([arm, intel])
        for t in ts[1]
            dft = df[df.num_threads .== t .&& df.num_cpp_chains .== 1, :]
            #println(dft)
            marksym = :cross
            plot!(dft.num_chains, dft.median;
                marker=marksym, lab="$(ptype[indx])_julia", color=colors[indx])
        end
        for t in ts[2]
            dft = df[df.num_threads .== t .&& df.num_chains .== 1, :]
            #println(dft)
            marksym = :xcross
            plot!(dft.num_cpp_chains, dft.median;
                marker=marksym, lab="$(ptype[indx])_c++", color=:grey)
        end
    end
    title!("Julia & c++\n(8 threads)")
    fig
end

f1 = timings_plot(arm_no_tbb, intel_tbb)

function threads_plot(df::DataFrame; ts = [1, 2, 4, 8])
    colors = [:darkred, :darkblue, :darkgreen, :black]
    fig = plot(; xlim=(0, 9), ylim=(0, 100),
        xlab="num_c++_chains", ylab="elapsed time [s]", leg=:topleft)
    for (indx, t) in enumerate(ts)
        dft = df[df.num_threads .== t .&& df.num_chains .== 1, :]
        #println(dft)
        marksym = :cross
        plot!(dft.num_cpp_chains, dft.median;
            marker=marksym, lab="threads = $(t)", color=colors[indx])
    end
    title!("TBB & c++\n(8 core Intel)")
    fig
end

f2 = threads_plot(intel_tbb)

fig = plot(f1, f2; layout=(1, 2))
savefig(joinpath(ProjDir, "graphs", "performance.png"))


text:
# Autogenerated wrapper script for Hanabi_jll for x86_64-linux-musl-cxx03
export compile_commands, libhanabi, libhanabi_h

JLLWrappers.@generate_wrapper_header("Hanabi")
JLLWrappers.@declare_file_product(compile_commands)
JLLWrappers.@declare_library_product(libhanabi, "libhanabi.so")
JLLWrappers.@declare_file_product(libhanabi_h)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        compile_commands,
        "share/compile_commands.json",
    )

    JLLWrappers.@init_library_product(
        libhanabi,
        "lib/libhanabi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_file_product(
        libhanabi_h,
        "include/pyhanabi.h",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
  HttpMethods

Http method types.
"""
module HttpMethods

export HttpMethod,
  Get,
  Post,
  Put,
  Patch,
  Delete

const _get = "GET"
const _post = "POST"
const _put = "PUT"
const _patch = "PATCH"
const _delete = "DELETE"

"""
  HttpMethod

Abstract type representing HTTP methods.
"""
abstract type HttpMethod end

Base.show(io::IO, x::HttpMethod) = write(io, string(x))

# http_method

"""
  http_method(::String)

Return type corresponding to string.
"""

function http_method(str::String)::HttpMethod
  cleaned = str |> strip |> uppercase
  if cleaned == _get
    return Get()
  elseif cleaned == _post
    return Post()
  elseif cleaned == _put
    return Put()
  elseif cleaned == _patch
    return Patch()
  elseif cleaned == _delete
    return Delete()
  end
end

# Get

struct Get <: HttpMethod end

Base.string(::Get) = "GET"

# Post

struct Post <: HttpMethod end

Base.string(::Post) = "POST"

# Put

struct Put <: HttpMethod end

Base.string(::Put) = "PUT"

# Patch

struct Patch <: HttpMethod end

Base.string(::Patch) = "PATCH"

# Delete

struct Delete <: HttpMethod end

Base.string(::Delete) = "DELETE"

end # module HttpMethods


text:
#!/usr/bin/env julia

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

using IntervalTrees
import IntervalTrees: Slice, InternalNode, LeafNode, Interval, IntervalBTree

# Getters
@testset "Getters" begin
    i = Interval(5, 6)
    iv = IntervalValue(10, 11, "FOO")
    @test first(i)  === 5
    @test last(i)   === 6
    @test first(iv) === 10
    @test last(iv)  === 11
    @test value(iv) == "FOO"
end

# Generating random intervals
srand(12345)
const maxend = round(Int, 1e9)

function randinterval(minstart=1, maxend=maxend)
    a = rand(minstart:maxend)
    b = rand(a:maxend)
    return (a, b)
end

# Verify that internal node keys and maxend values are correct
function validkeys{K, V, B}(t::IntervalTrees.IntervalBTree{K, V, B})
    minint = IntervalTrees.Interval{K}(0, 0)
    maxint = IntervalTrees.Interval{K}(maxend, maxend)
    return validkeys(t.root, minint, maxint)
end

function validkeys(node::IntervalTrees.InternalNode, minint, maxint)
    for child in node.children
        if child.maxend > node.maxend
            return false
        end
    end

    for i in 1:length(node.keys)
        k = node.keys[i]
        if IntervalTrees.minkey(node.children[i+1]) != k
            return false
        end

        if !validkeys(node.children[i], minint, k) ||
           !validkeys(node.children[i+1], k, maxint)
           return false
       end
    end
    return true
end

function validkeys(node::IntervalTrees.LeafNode, minint, maxint)
    for entry in node.entries
        if !(minint <= entry <= maxint) || last(entry) > node.maxend
            return false
        end
    end
    return true
end

# Verify that parent pointers are correct
function validparents(t::IntervalTrees.IntervalBTree)
    if !isnull(t.root.parent)
        return false
    end

    return validparents(t.root)
end

function validparents(node::IntervalTrees.InternalNode)
    for child in node.children
        if !validparents(child) || get(child.parent) != node
            return false
        end
    end
    return true
end

function validparents(node::IntervalTrees.LeafNode)
    return true
end

# Verify that sibling/cousin pointers are correct
function validsiblings(t::IntervalTrees.IntervalBTree)
    if !isnull(t.root.left) || !isnull(t.root.right)
        return false
    end

    # Do an in-order traversal, pushing nodes onto a stack indexed by
    # depth.
    tdepth = depth(t)
    nodestacks = [Array(IntervalTrees.Node, 0) for _ in 1:tdepth]

    function _visit(node::IntervalTrees.InternalNode, k)
        push!(nodestacks[k], node)
        for child in node.children
            _visit(child, k+1)
        end
    end

    function _visit(node::IntervalTrees.LeafNode, k)
        push!(nodestacks[k], node)
    end

    _visit(t.root, 1)

    for k in 1:tdepth
        stack = nodestacks[k]
        if isempty(stack)
            continue
        end

        if !isnull(stack[1].left) || !isnull(stack[end].right)
            return false
        end

        if length(stack) > 1
            if isnull(stack[1].right) || get(stack[1].right) != stack[2] ||
               isnull(stack[end].left) || get(stack[end].left) != stack[end-1]
                return false
            end
        end

        for i in 2:length(stack)-1
            if isnull(stack[i].left) || get(stack[i].left) != stack[i-1] ||
               isnull(stack[i].right) || get(stack[i].right) != stack[i+1]
                return false
            end
        end
    end

    return true
end

@testset "Search" begin
    t = IntervalMap{Int, Int}()
    @test !haskey(t, (1,2))

    n = 10000
    maxend = 1000000
    intervals = [randinterval(1, maxend) for i in 1:n]

    @testset "true positives" begin
        for (i, interval) in enumerate(intervals)
            t[interval] = i
        end

        shuffle!(intervals)
        @test all(Bool[haskey(t, interval) for interval in intervals])
        @test all(Bool[get!(t, interval, -1) != IntervalValue{Int, Int}(interval[1], interval[2], -1)
                       for interval in intervals])
        @test all(Bool[get(t, interval, -1) != -1 for interval in intervals])
    end

    @testset "true negatives" begin
        @test all(Bool[!haskey(t, interval)
                       for interval in [randinterval(maxend+1, 2 * maxend)
                                        for i in 1:n]])
        @test all(Bool[get!(t, interval, -1) == IntervalValue{Int, Int}(interval[1], interval[2], -1)
                       for interval in [randinterval(maxend+1, 2 * maxend)
                                        for i in 1:n]])
        @test all(Bool[get(t, interval, -1) == -1
                       for interval in [randinterval(maxend+1, 2 * maxend)
                                        for i in 1:n]])
    end
end

@testset "Iteration" begin
    t = IntervalMap{Int, Int}()
    @test isempty([x for x in t])

    @testset "from" begin
        n = 100
        maxend = 1000000
        intervals = [randinterval(1, maxend) for _ in 1:n]
        startpos = 50000
        expected_count = 0
        t = IntervalMap{Int, Int}()
        for (i, interval) in enumerate(intervals)
            t[interval] = i
            if interval[2] >= startpos
                expected_count += 1
            end
        end

        @test length(collect(from(t, startpos))) === expected_count
        @test length(collect(from(t, 0))) === n
        @test length(collect(from(t, maxend + 1))) === 0
        @test length(collect(from(IntervalTree{Int, Int}(), 0))) === 0
    end
end

@testset "Interval Intersection" begin
    # generate n end-to-end intervals
    t = IntervalMap{Int, Int}()
    intervals = Any[]
    a = 1
    for i in 1:10000
        b = a + rand(0:100)
        t[(a, b)] = i
        push!(intervals, (a, b))
        a = b + 1
    end

    # one
    x = rand(0:a-1)
    @test length(collect(intersect(t, (x, x)))) === 1
    @test hasintersection(t, x)

    # nothing
    @test length(collect(intersect(t, (a, a)))) === 0
    @test !hasintersection(t, a)

    # everything
    @test length(collect(intersect(t, (intervals[1][1], intervals[end][2])))) == length(t)

    @test all(Bool[hasintersection(t, interval[1]) for interval in intervals])
    @test all(Bool[hasintersection(t, interval[end]) for interval in intervals])

    # some random intersection queries
    function random_intersection_query()
        i = rand(1:length(intervals))
        j = rand(i:length(intervals))
        return length(collect(intersect(t, (intervals[i][1], intervals[j][2])))) == j - i + 1
    end

    @test all(Bool[random_intersection_query() for _ in 1:1000])

    # intervals separated by 1
    t = IntervalMap{Int, Int}()

    @test !hasintersection(t, 1)

    intervals = Any[]
    gaps = Any[]
    a = 1
    for i in 1:10000
        b = a + rand(0:100)
        if iseven(b)
            t[a, b] = i
            push!(intervals, (a, b))
        else
            push!(gaps, (a, b))
        end
        a = b + 1
    end

    @test all(Bool[hasintersection(t, interval[1]) for interval in intervals])
    @test all(Bool[hasintersection(t, interval[end]) for interval in intervals])
    @test !all(Bool[hasintersection(t, interval[1]) for interval in gaps])
    @test !all(Bool[hasintersection(t, interval[end]) for interval in gaps])
end

@testset "Nonunique" begin
    t = IntervalTree{Int, IntervalValue{Int, Int}}()
    n = 1000
    for i in 1:n
        push!(t, IntervalValue{Int, Int}(50, 100, i))
    end
    push!(t, IntervalValue{Int, Int}(1, 200, n + 1))
    push!(t, IntervalValue{Int, Int}(1, 1, n + 2))
    push!(t, IntervalValue{Int, Int}(200, 200, n + 3))

    @test (length(collect(t)) == length(t) == n + 3)
    @test issorted(collect(keys(t)))
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
end

@testset "Tree Intersection" begin
    n = 10000
    t1 = IntervalTrees.IntervalMap{Int, Int}()
    t2 = IntervalTrees.IntervalMap{Int, Int}()
    maxend = 1000000
    for k in 1:n
        # generate small-ish intervals so we avoid the worst case
        # of O(n^2) intersecting pairs and this test runs quickly
        u = rand(1:maxend)
        v = rand(u:u+1000)
        t1[(u,v)] = k

        u = rand(1:maxend)
        v = rand(u:u+1000)
        t2[(u,v)] = k
    end

    # Since IntervalTrees has two tree intersection algorithms, I'm
    # testing by checking that they are in agreement.
    a = intersect(t1, t2, method=:successive)
    b = intersect(t1, t2, method=:iterative)
    @test collect(a) == collect(b)

    ## handle a particular intersection case that may not otherwise get hit
    t1 = IntervalMap{Int, Int}()
    t1[(1, 2)] = 1
    t2 = IntervalMap{Int, Int}()
    t2[(1001, 1002)] = 2
    @test isempty(collect(intersect(t1, t2, method=:iterative)))
    @test isempty(collect(intersect(t1, t2, method=:successive)))
end

@testset "Insertion" begin
    t = IntervalMap{Int, Int}()
    n = 100000

    @testset "random insertions" begin
        for v in 1:n
            k = randinterval()
            t[k] = v
        end
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)
    end

    @testset "ordered insertions" begin
        for v in 1:n
            k = (v,rand(v:maxend))
            t[k] = v
        end
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)
    end

    @testset "reverse ordered insertions" begin
        for v in 1:n
            k = (maxend-v,rand((maxend-v):maxend))
            t[k] = v
        end
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)
    end

    @testset "bulk insertion" begin
        intervals = IntervalValue{Int, Int}[]
        for v in 1:n
            a, b = randinterval()
            push!(intervals, IntervalValue{Int, Int}(a, b, v))
        end
        sort!(intervals)
        t = IntervalMap{Int, Int}(intervals)

        @test collect(keys(t)) == [Interval{Int}(interval.first, interval.last)
                                   for interval in intervals]
        @test collect(values(t)) == [interval.value for interval in intervals]
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)

        # don't break on empty arrays
        t = IntervalMap{Int, Int}(IntervalValue{Int, Int}[])
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)

        @test_throws MethodError IntervalTree{Int, Int}(intervals, Int[1])
        shuffle!(intervals)
        @test_throws UndefVarError IntervalTree{Int, Int}(intervals, vals)
    end
end

@testset "Updates" begin
    t = IntervalMap{Int, Int}()
    n = 100000
    ks = [randinterval() for _ in 1:n]
    for (v, k) in enumerate(ks)
        t[k] = v
    end

    shuffle!(ks)
    for (v, k) in enumerate(ks)
        t[k] = v
    end

    @test issorted(collect(keys(t)))
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
end

@testset "Deletion" begin
    # Set B to smaller number to try to induce more weird merges and borrows and
    # such
    t = IntervalBTree{Int, IntervalValue{Int, Int}, 8}()
    n = 100000

    # insert n random intervals
    intervals = [randinterval() for _ in 1:n]
    for (i, interval) in enumerate(intervals)
        t[interval] = i
    end
    shuffle!(intervals)

    # delete non-existant is a nop
    @test delete!(t, (-1, -1)).n === n

    # delete one
    interval = pop!(intervals)
    delete!(t, interval)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
    @test issorted(collect(keys(t)))
    @test !haskey(t, interval)
    @test all(interval -> haskey(t, interval), intervals)

    # delete a random 50%
    for interval in 1:(n-1)/2
        interval = pop!(intervals)
        delete!(t, interval)
    end
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
    @test issorted(collect(keys(t)))
    @test all(interval -> haskey(t, interval), intervals)

    # delete the rest
    for interval in intervals
        delete!(t, interval)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # delete left-to-right
    shuffle!(intervals)
    for (i, interval) in enumerate(intervals)
        t[interval] = i
    end
    sort!(intervals)
    for interval in intervals
        delete!(t, interval)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # delete right-to-left
    shuffle!(intervals)
    for (i, interval) in enumerate(intervals)
        t[interval] = i
    end
    reverse!(intervals)
    for interval in intervals
        delete!(t, interval)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # delete from the middle
    shuffle!(intervals)
    for (i, interval) in enumerate(intervals)
        t[interval] = i
    end
    mid = div(length(intervals), 2)
    for i in mid+1:length(intervals)
        delete!(t, intervals[i])
    end
    for i in mid:-1:1
        delete!(t, intervals[i])
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # contrive a special case: merge with the left when we have a non-null
    # right
    t = IntervalBTree{Int, IntervalValue{Int, Int}, 4}()
    for i in 1:24
        t[(i, i)] = i
    end
    @test isa(t.root, InternalNode)
    @test length(t.root.children) === 3
    @test length(t.root.children[2].children) === 2
    delete!(t, (15, 15))
    keys_to_delete = collect(keys(t))
    for key in keys_to_delete
        delete!(t, key)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
end

@testset "Slices" begin
    xs = Slice{Int, 10}()

    @test_throws BoundsError xs[0]
    @test_throws BoundsError xs[11]
    @test_throws BoundsError xs[0] = 1
    @test_throws BoundsError xs[11] = 1
    @test_throws BoundsError pop!(xs)
    @test_throws BoundsError insert!(xs, 0, 1)

    for x in 1:10
        push!(xs, x)
    end

    @test_throws BoundsError push!(xs, 11)
    @test_throws BoundsError insert!(xs, 11, 1)
    @test_throws BoundsError resize!(xs, 11)
    @test_throws BoundsError splice!(xs, 0)
    @test_throws BoundsError splice!(xs, 11)
end

# Abuse low-level functions to test cases that wouldn't otherwise occur
@testset "Low Level" begin
    # findidx should return 0 when called on an empty node
    K = Int
    V = IntervalValue{Int}
    B = 32
    x = Interval{Int}(1, 1)
    node = InternalNode{K, V, B}()
    @test IntervalTrees.findidx(node, x) === 0
    @test IntervalTrees.firstfrom(node, 1) == (node, 0)
    node = LeafNode{K, V, B}()
    @test IntervalTrees.findidx(node, x) == 0

    result = IntervalTrees.Intersection{K, V, B}()
    IntervalTrees.firstintersection!(node, x, result)
    @test result.index == 0

    @test IntervalTrees.firstfrom(node, 1) == (node, 0)

    push!(node.entries, IntervalValue{Int, Int}(1, 1, 1))
    @test IntervalTrees.firstfrom(node, 2) == (node, 0)

    # test that delete! still works on a contrived tree with one internal and
    # one leaf node
    t = IntervalBTree{K, V, B}()
    t.root = InternalNode{K, V, B}()
    push!(t.root.children, LeafNode{K, V, B}())
    push!(t.root.keys, x)
    push!(t.root.children[1].entries, IntervalValue{Int, Int}(1, 1, 1))
    delete!(t, (1,1))
    @test isa(t.root, LeafNode)

    ## test that the right thing happens if you delete the last key in a non-root
    ## leaf node (which can't actually happen)
    node = LeafNode{Int, IntervalValue{Int}, 32}()
    push!(node.entries, IntervalValue{Int, Int}(1, 1, 1))
    @test IntervalTrees.findidx(node, x) == 1
end


text:
using RunTests
using Base.Test
import MetaTools.@commutative

@testmodule CommutativeTests begin
  using MetaTools

  module MetaTestModule1
    function baz()
    end
    function boo()
    end
  end

  function commutative_method()
    @commutative function foo(s::String, i::Integer)
      return "$s -> $i"
    end

    @test length(methods(foo)) == 2
    @test foo("hello", 42) == "hello -> 42"
    @test foo(42, "hello") == "hello -> 42"
  end

  function commutative_with_wrong_type_of_function()
    @test_throws eval(:(@commutative function foo(s::String) end))
    @test_throws eval(:(@commutative function foo(s1::String, s2::String) end))
    @test_throws eval(:(@commutative function foo(s::String, i::Integer...) end))
    @test_throws eval(:(@commutative function foo(s::String, i::Integer=10) end))
    @test_throws eval(:(@commutative function foo(s1::String, s2::Integer, s3::FloatingPoint) end))
    @test_throws eval(:(@commutative function foo(s::String, i::Integer; special=false) end))
    @test_throws eval(:(@commutative function(s::String, i::Integer) end))
    @test_throws eval(:(@commutative (s::String, i::Integer) -> begin end))
  end

  function commutative_shorthand_method()
    @commutative bar(s::String, i::Integer) = "$s -> $i"

    @test length(methods(bar)) == 2
    @test bar("hello", 42) == "hello -> 42"
    @test bar(42, "hello") == "hello -> 42"
  end

  function commutative_method_in_namespace()
    @commutative function MetaTestModule1.baz(s::String, i::Integer)
      return "$s -> $i"
    end

    @test MetaTestModule1.baz("hello", 42) == "hello -> 42"
    @test MetaTestModule1.baz(42, "hello") == "hello -> 42"
  end

  function commutative_shorthand_method_in_namespace()
    @commutative MetaTestModule1.boo(s::String, i::Integer) = "$s -> $i"

    @test MetaTestModule1.boo("hello", 42) == "hello -> 42"
    @test MetaTestModule1.boo(42, "hello") == "hello -> 42"
  end

  function commutative_method_with_an_any_arg_makes_3_methods()
    @commutative function hmmm(a::Any, s::String)
      return "$s -> $a"
    end

    @test length(methods(hmmm)) == 3
    @test hmmm("hello", "world") == "hello -> world"
    @test hmmm("hello", 42) == "hello -> 42"
    @test hmmm(42, "hello") == "hello -> 42"
  end

end


text:
@deprecate BullGraph bull_graph
@deprecate ChvatalGraph chvatal_graph
@deprecate CubicalGraph cubical_graph
@deprecate DesarguesGraph desargues_graph
@deprecate DiamondGraph diamond_graph
@deprecate DodecahedralGraph dodecahedral_graph
@deprecate FruchtGraph frucht_graph
@deprecate HeawoodGraph heawood_graph
@deprecate HouseGraph house_graph
@deprecate HouseXGraph house_x_graph
@deprecate IcosahedralGraph icosahedral_graph
@deprecate KarateGraph karate_graph
@deprecate KrackhardtKiteGraph krackhardt_kite_graph
@deprecate MoebiusKantorGraph moebius_kantor_graph
@deprecate OctahedralGraph octahedral_graph
@deprecate PappusGraph pappus_graph
@deprecate PetersenGraph petersen_graph
@deprecate SedgewickMazeGraph sedgewick_maze_graph
@deprecate TetrahedralGraph tetrahedral_graph
@deprecate TruncatedCubeGraph truncated_cube_graph
@deprecate TruncatedTetrahedronGraph truncated_tetrahedron_graph
@deprecate TruncatedTetrahedronDiGraph truncated_tetrahedron_digraph
@deprecate TutteGraph tutte_graph
@deprecate CompleteGraph complete_graph
@deprecate CompleteBipartiteGraph complete_bipartite_graph
@deprecate CompleteMultipartiteGraph complete_multipartite_graph
@deprecate TuranGraph turan_graph
@deprecate CompleteDiGraph complete_digraph
@deprecate StarGraph star_graph
@deprecate StarDigraph star_digraph
@deprecate PathGraph path_graph
@deprecate PathDiGraph path_digraph
@deprecate CycleGraph cycle_graph
@deprecate CycleDiGraph cycle_digraph
@deprecate WheelGraph wheel_graph
@deprecate WheelDiGraph wheel_digraph
@deprecate Grid grid
@deprecate BinaryTree binary_tree
@deprecate Doublebinary_tree double_binary_tree
@deprecate RoachGraph roach_graph
@deprecate CliqueGraph clique_graph
@deprecate LadderGraph ladder_graph
@deprecate Circularladder_graph circular_ladder_graph
@deprecate BarbellGraph barbell_graph
@deprecate LollipopGraph lollipop_graph


text:
function AmericanSpreadCallBinomial(K, T, S1, S2, sig1, sig2, div1, div2, rho, r, N)

  dt = T/N ## Time to Maturity/ Number of Steps
  nu1 = r - div1 - 0.5 * sig1^2
  nu2 = r - div2 - 0.5 * sig2^2

  dx1 = sig1 * sqrt(dt)
  dx2 = sig2 * sqrt(dt)

  disc = exp(-r * dt)

  puu = (dx1 * dx2 + (dx2 * nu1 + dx1 * nu2 + rho * sig1 * sig2) * dt) / (4 * dx1 * dx2) * disc

  pud = (dx1 * dx2 + (dx2 * nu1 - dx1 * nu2 - rho * sig1 * sig2) * dt) / (4 * dx1 * dx2) * disc

  pdu = (dx1 * dx2 + (-dx2 * nu1 + dx1 * nu2 - rho * sig1 * sig2) * dt) / (4 * dx1 * dx2) * disc

  pdd = (dx1 * dx2 + (-dx2 * nu1 - dx1 * nu2 + rho * sig1 * sig2) * dt) / (4 * dx1 * dx2) * disc

  edx1 = exp(dx1)
  edx2 = exp(dx2)

  S1t = linspace(-N, N, 2N+1)
  S2t = linspace(-N, N, 2N+1)

  S1t[1] = S1 * exp(-N * dx1)
  S2t[1] = S2 * exp(-N * dx2)

  for j in 2:length(-N:N)
    S1t[j] = S1t[j-1] * edx1
    S2t[j] = S2t[j-1] * edx2
  end

  C = Array(Float64, (2N+1,2N+1))

  for j = 1:2:length(-N:N)
    for k = 1:2:length(-N:N)
      C[k, j] = max(0.0, S1t[j] - S2t[k] - K)
    end
  end

  for i in range(N-1, -1, N)
    for j in range(-i, 2, i+1) # potentially this is i^2
      for k in range(-i, 2, i+1) # potentiall this is i^2 as well
        C[k+(N+1), j+(N+1)] = pdd * C[k+(N+1)-1, j+(N+1)-1] + pud * C[k+(N+1)-1, j+(N+1)+1] + pdu * C[k+(N+1)+1, j+(N+1)-1] + puu * C[k+(N+1)+1, j+(N+1)+1]
        C[k+(N+1), j+(N+1)] = max(C[k+(N+1), j+(N+1)], S1t[j+(N+1)] - S2t[k+(N+1)] - K)
      end
    end
  end

  C[N+1, N+1]

end


text:
# Functions for efficiently evaluating all non-zero B-splines (and/or their
# derivatives) at a given point.

using Base.Cartesian: @ntuple
using Base: @propagate_inbounds
using StaticArrays: MVector

"""
    evaluate_all(
        B::BSplineBasis, x::Real,
        [op = Derivative(0)], [T = float(typeof(x))];
        [ileft = nothing],
    ) -> i, bs

Evaluate all B-splines which are non-zero at coordinate `x`.

Returns a tuple `(i, bs)`, where `i` is an index identifying the basis functions
that were computed, and `bs` is a tuple with the actual values.

More precisely:

- `i` is the index of the first B-spline knot ``t_{i}`` when going from ``x``
  towards the left.
  In other words, it is such that ``t_{i} ≤ x < t_{i + 1}``.

  It is effectively computed as `i = searchsortedlast(knots(B), x)`.
  If the correct value of `i` is already known, one can avoid this computation by
  manually passing this index via the optional `ileft` keyword argument.

- `bs` is a tuple of B-splines evaluated at ``x``:

  ```math
  (b_i(x), b_{i - 1}(x), …, b_{i - k + 1}(x)).
  ```

  It contains ``k`` values, where ``k`` is the order of the B-spline basis.
  Note that values are returned in backwards order starting from the ``i``-th
  B-spline.

## Computing derivatives

One can pass the optional `op` argument to compute B-spline derivatives instead
of the actual B-spline values.

## Examples

See [`AbstractBSplineBasis`](@ref) for some examples using the alternative
evaluation syntax `B(x, [op], [T]; [ileft])`, which calls this function.
"""
@propagate_inbounds function evaluate_all(
        B::BSplineBasis, x::Real, op::Derivative, ::Type{T}; kws...,
    ) where {T <: Number}
    _evaluate_all(knots(B), x, BSplineOrder(order(B)), op, T; kws...)
end

@propagate_inbounds evaluate_all(
    B, x, op::AbstractDifferentialOp = Derivative(0); kws...,
) = evaluate_all(B, x, op, float(typeof(x)); kws...)

@propagate_inbounds evaluate_all(B, x, ::Type{T}; kws...) where {T <: Number} =
    evaluate_all(B, x, Derivative(0), T; kws...)

@propagate_inbounds function _knotdiff(x, ts, i, n)
    @boundscheck checkbounds(ts, i:(i + n))
    @inbounds ti = ts[i]
    @inbounds tj = ts[i + n]
    # @assert ti ≠ tj
    (x - ti) / (tj - ti)
end

# TODO
# - this is redundant with Splines.knot_interval...
"""
    find_knot_interval(ts::AbstractVector, x::Real) -> (i, zone)

Finds the index ``i`` corresponding to the knot interval ``[t_i, t_{i + 1}]``
that should be used to evaluate B-splines at location ``x``.

The knot vector is assumed to be sorted in non-decreasing order.

It also returns a `zone` integer, which is:

- `0`  if `x` is within the knot domain (`ts[begin] ≤ x ≤ ts[end]`),
- `-1` if `x < ts[begin]`,
- `1`  if `x > ts[end]`.

This function is functionally equivalent to de Boor's `INTERV` routine (de Boor
2001, p. 74).
"""
function find_knot_interval(ts::AbstractVector, x::Real)
    if x < first(ts)
        return firstindex(ts), -1
    end
    i = searchsortedlast(ts, x)
    Nt = lastindex(ts)
    if i == Nt
        tlast = ts[Nt]
        while true
            i -= 1
            ts[i] ≠ tlast && break
        end
        zone = (x > tlast) ? 1 : 0
        return i, zone
    else
        return i, 0  # usual case
    end
end

function _evaluate_all(
        ts::AbstractVector, x::Real, ::BSplineOrder{k},
        op::Derivative{0}, ::Type{T};
        ileft = nothing,
    ) where {k, T}
    if @generated
        @assert k ≥ 1
        ex = quote
            i, zone = _find_knot_interval(ts, x, ileft)
            if zone ≠ 0
                return (i, @ntuple($k, j -> zero($T)))
            end
            # We assume that the value of `i` corresponds to the good interval,
            # even if `ileft` was passed by the user.
            # This is the same as in de Boor's BSPLVB routine.
            # In particular, this allows plotting the extension of a polynomial
            # piece outside of its knot interval (as in de Boor's Fig. 6, p. 114).
            bs_1 = (one(T),)
        end
        for q ∈ 2:k
            bp = Symbol(:bs_, q - 1)
            bq = Symbol(:bs_, q)
            ex = quote
                $ex
                Δs = @ntuple(
                    $(q - 1),
                    j -> @inbounds($T(_knotdiff(x, ts, i - j + 1, $q - 1))),
                )
                $bq = _evaluate_step(Δs, $bp, BSplineOrder($q))
            end
        end
        bk = Symbol(:bs_, k)
        quote
            $ex
            return i, $bk
        end
    else
        _evaluate_all_alt(ts, x, BSplineOrder(k), op, T; ileft = ileft)
    end
end

# Derivatives
function _evaluate_all(
        ts::AbstractVector, x::Real, ::BSplineOrder{k},
        op::Derivative{n}, ::Type{T};
        ileft = nothing,
    ) where {k, n, T}
    if @generated
        n::Int
        @assert n ≥ 1
        # We first need to evaluate the B-splines of order p.
        p = k - n
        if p < 1
            # Derivatives are zero. The returned index is arbitrary...
            return :( firstindex(ts), ntuple(_ -> zero($T), Val($k)) )
        end
        bp = Symbol(:bs_, p)
        ex = quote
            i, $bp = _evaluate_all(ts, x, BSplineOrder($p), Derivative(0), $T; ileft)
        end
        for q ∈ (p + 1):k
            bp = Symbol(:bs_, q - 1)
            bq = Symbol(:bs_, q)
            ex = quote
                $ex
                $bq = _evaluate_step_deriv(ts, i, $bp, BSplineOrder($q), $T)
            end
        end
        bk = Symbol(:bs_, k)
        quote
            $ex
            return i, $bk
        end
    else
        _evaluate_all_alt(ts, x, BSplineOrder(k), op, T; ileft = ileft)
    end
end

# Non-@generated version
@inline function _evaluate_all_alt(
        ts::AbstractVector, x::Real, ::BSplineOrder{k},
        ::Derivative{0}, ::Type{T};
        ileft = nothing,
    ) where {k, T}
    @assert k ≥ 1
    i, zone = _find_knot_interval(ts, x, ileft)
    if zone ≠ 0
        return (i, ntuple(j -> zero(T), Val(k)))
    end
    bq = MVector{k, T}(undef)
    Δs = MVector{k - 1, T}(undef)
    bq[1] = one(T)
    @inbounds for q ∈ 2:k
        for j ∈ 1:(q - 1)
            Δs[j] = _knotdiff(x, ts, i - j + 1, q - 1)
        end
        bp = bq[1]
        Δp = Δs[1]
        bq[1] = Δp * bp
        for j = 2:(q - 1)
            bpp, bp = bp, bq[j]
            Δpp, Δp = Δp, Δs[j]
            bq[j] = Δp * bp + (1 - Δpp) * bpp
        end
        bq[q] = (1 - Δp) * bp
    end
    i, Tuple(bq)
end

# Derivatives / non-@generated variant
function _evaluate_all_alt(
        ts::AbstractVector, x::Real, ::BSplineOrder{k},
        ::Derivative{n}, ::Type{T};
        ileft = nothing,
    ) where {k, n, T}
    n::Int
    @assert n ≥ 1
    # We first need to evaluate the B-splines of order m.
    m = k - n
    if m < 1
        # Derivatives are zero. The returned index is arbitrary...
        return firstindex(ts), ntuple(_ -> zero(T), Val(k))
    end
    i, bp = _evaluate_all(ts, x, BSplineOrder(m), Derivative(0), T; ileft)
    bq = MVector{k, T}(undef)
    us = MVector{k - 1, T}(undef)
    bq[1:m] .= bp
    for q = (m + 1):k
        p = q - 1
        for δj = 1:p
            @inbounds us[δj] = bq[δj] / (ts[i + q - δj] - ts[i + 1 - δj])
        end
        @inbounds bq[1] = p * us[1]
        for j = 2:(q - 1)
            # Note: adding @inbounds here can slow down stuff (from 25ns to
            # 80ns), which is very strange!! (Julia 1.8-beta2)
            bq[j] = p * (us[j] - us[j - 1])
        end
        @inbounds bq[q] = -p * us[p]
    end
    i, Tuple(bq)
end

function _find_knot_interval(ts, x, ileft)
    if isnothing(ileft)
        i, zone = find_knot_interval(ts, x)
    else
        i = ileft
        zone = (x < first(ts)) ? -1 : (x > last(ts)) ? 1 : 0
    end
    i, zone
end

@inline @generated function _evaluate_step(Δs, bp, ::BSplineOrder{k}) where {k}
    ex = quote
        @inbounds b_1 = Δs[1] * bp[1]
    end
    for j = 2:(k - 1)
        bj = Symbol(:b_, j)
        ex = quote
            $ex
            @inbounds $bj = (1 - Δs[$j - 1]) * bp[$j - 1] + Δs[$j] * bp[$j]
        end
    end
    b_last = Symbol(:b_, k)
    quote
        $ex
        @inbounds $b_last = (1 - Δs[$k - 1]) * bp[$k - 1]
        @ntuple $k b
    end
end

@inline @generated function _evaluate_step_deriv(
        ts, i, bp, ::BSplineOrder{k}, ::Type{T},
    ) where {k, T}
    p = k - 1
    ex = quote
        us = @ntuple(
            $p,
            δj -> @inbounds($T(bp[δj] / (ts[i + $k - δj] - ts[i + 1 - δj]))),
        )
        @inbounds b_1 = $p * us[1]
    end
    for j = 2:(k - 1)
        bj = Symbol(:b_, j)
        ex = quote
            $ex
            @inbounds $bj = $p * (-us[$j - 1] + us[$j])
        end
    end
    b_last = Symbol(:b_, k)
    quote
        $ex
        @inbounds $b_last = -$p * us[$p]
        @ntuple $k b
    end
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
### documentation build script
debug = false
JDK_HOME = "E:\\Compilatori\\Java64.180u25"
PYTHON_HOME = "E:\\Compilatori\\Python2710_64"

### step 1: create .md file from documentation mixed into source files (with Lexicon)
# set path to JDK home into JAVA_HOME environment variable
if !("JAVA_HOME" in keys(ENV))
	ENV["JAVA_HOME"] = JDK_HOME
end
pdir = Pkg.dir()
include(pdir * "/WinReg/docs/build_md.jl")

### step 2: create static html site from .md files (with MkDocs)
# get path to PYTHON 2.7.x home
if "PYTHON_HOME" in keys(ENV)
	pyhome = ENV["PYTHON_HOME"]
else
	pyhome = PYTHON_HOME
end
pdir = Pkg.dir()
mkdocs = PYTHON_HOME * "/Scripts/mkdocs"
cd(pdir * "/WinReg")

if debug
	#serve documentation with local web server
	run(`$mkdocs serve`)
else
	#generate documentation
	run(`$mkdocs build --clean`)
end


text:
include("automata.jl")

const λ = 'λ'

mutable struct NFA <: AbstractFiniteAutomaton
    cur_states::Set{State}
    aut::Automaton
end

function NFA(aut::Automaton)
    NFA(Set([aut.q₀]), aut)
end

NFA(δ::AbstractDict,
    q₀::State,
    F::AbstractSet{State}) = NFA(Automaton(δ, q₀, F))

function try_λ_transitions!(M::NFA)
    states_to_try = copy(M.cur_states)
    while !isempty(states_to_try)
        state = pop!(states_to_try)
        key = (state, λ)
        if haskey(M.aut.δ, key)
            next_states = M.aut.δ[key]
            push!(M.cur_states, next_states...)
            push!(states_to_try, next_states...)
        end
    end
end

function transition!(M::NFA, a)::Bool
    old_states = copy(M.cur_states)
    empty!(M.cur_states)
    for state in old_states
        key = (state, a)
        if haskey(M.aut.δ, key)
            push!(M.cur_states, M.aut.δ[key]...)
        end
    end

    isempty(M.cur_states) ? false : true
end

function reset!(M::NFA)
    M.cur_states = Set([M.aut.q₀])
end

function accepting(M::NFA)
    !isdisjoint(M.cur_states, M.aut.F)
end

function evaluate!(M::NFA, w)
    reset!(M)

    try_λ_transitions!(M)

    for a in w
        transition!(M, a) || return false
        try_λ_transitions!(M)
    end

    accepting(M)
end

max_issues_repo_name: MadhavJivrajani/DFA-Sim

max_forks_repo_name: MadhavJivrajani/DFA-Sim


text:
using Arrhenius
using ForwardDiff
using LinearAlgebra
using Test

@testset "jl" begin
    # Write your tests here.

    gas = CreateSolution("../mechanism/gri30.yaml")
    ns = gas.n_species

    Y0 = ones(ns) ./ ns
    T0 = 1200.0
    P = one_atm
    wdot = set_states(gas, T0, P, Y0)
    @show size(wdot)
    @show wdot[1:3]

    u0 = vcat(Y0, T0)
    function f(u)
        T = u[end]
        Y = @view(u[1:ns])
        mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
        ρ_mass = P / R / T * mean_MW
        X = Y2X(gas, Y, mean_MW)
        C = Y2C(gas, Y, ρ_mass)
        cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
        h_mole = get_H(gas, T, Y, X)
        S0 = get_S(gas, T, P, X)
        wdot = wdot_func(gas.reaction, T, C, S0, h_mole)
        Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass
        return Tdot
    end
    grad = ForwardDiff.gradient(f, u0)
    @show size(grad)
    @show grad[1:3]

end


text:
@everywhere using SeisDownload
using Dates

#==================================================#
# Input Parameters
MAX_MEM_PER_CPU = 2.0 # [GB] maximum allocated memory for one cpu
DownloadType = "Noise" # Choise of "Noise" or "Earthquake"

network     = ["BP"]
station     = ["LCCB", "MMNB"]
#station = ["CCRB","EADB","FROB","GHIB","JCNB","JCSB","LCCB","MMNB","SCYB","SMNB","VARB","VCAB"]

location    = [""]
channel     = ["BP1"]
datacenter  = "FDSN" #Data center
src         = "NCEDC"

# Time info for Noise case
starttime   = DateTime(2004,9,25,0,0,0)
endtime     = DateTime(2004,9,25,12,0,0)

IsLocationBox = false
method  = "FDSN" # Method to download data.
datasource = "NCEDC" # currently, only one src can be specified.

DL_time_unit = 3600 * 2 #3600 * 24 # Download tiem unit [s] more than one day is better to avoid artifacts of response removal

IsResponseRemove = true #whether instrumental response is removed or not
pre_filt    = (0.001, 0.002, 10.0, 20.0) #prefilter tuple used obspy remove_response: taper between f1 and f2, f3 and f4 with obspy

fodir       = "./dataset"
foname      = "BPnetwork" # data is saved at ./dataset/$foname.jld2
#==================================================#

# store metadata in Dictionary
# This can be customized by users

stationlist       = String[]
stationmethod    = String[]
stationsrc        = String[]
for i=1:length(network)
    for j=1:length(station)
        for k=1:length(location)
            for l=1:length(channel)
                stationname = join([network[i], station[j], location[k], channel[l]], ".")
                push!(stationlist, stationname)

                #Here should be improved for multiple seismic network; we have to make
                #proper conbination of request station and data server.
                push!(stationmethod, method)
                push!(stationsrc, datasource)
            end
        end
    end
end

stationinfo = Dict(["stationlist" => stationlist, "stationmethod" => stationmethod, "stationsrc" => stationsrc])

mkpath(fodir)
fopath=joinpath(fodir, foname*".jld2")

#if lat-log box or not
IsLocationBox ? reg=locationbox : reg=[]
IsResponseRemove ? pre_filt = pre_filt : pre_filt = []

InputDictionary = Dict([
        "DownloadType"=> DownloadType,
        "stationinfo" => stationinfo,
        "starttime"   => starttime,
        "endtime"     => endtime,
        "DL_time_unit"=> DL_time_unit,
        "IsLocationBox"   => IsLocationBox,
        "reg"             => reg,
        "IsResponseRemove"=> IsResponseRemove,
        "pre_filt"        => pre_filt,
        "fopath"          => fopath
    ])

# mass request with input Dictionary
seisdownload(InputDictionary, MAX_MEM_PER_CPU=float(MAX_MEM_PER_CPU))


text:
@ghdef mutable struct Installation
    id::Union{Int, Nothing}
end

namefield(i::Installation) = i.id

@api_default function create_access_token(api::GitHubAPI, i::Installation, auth::JWTAuth; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    payload = gh_post_json(api, "/installations/$(i.id)/access_tokens", auth = auth,
        headers=headers, options...)
    OAuth2(payload["token"])
end

@api_default function installations(api::GitHubAPI, auth::JWTAuth; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    results, page_data = gh_get_paged_json(api, "/app/installations", auth = auth,
        headers=headers, options...)
    map(Installation, results), page_data
end

@api_default function repos(api::GitHubAPI, inst::Installation; headers = Dict(), options...)
    headers["Accept"] = "application/vnd.github.machine-man-preview+json"
    results, page_data = github_paged_get(api, "/installation/repositories";
        headers=headers, options...)
    mapreduce(x->map(Repo, JSON.parse(HTTP.payload(x, String))["repositories"]), vcat, results; init=Repo[]), page_data
end


text:
using Printf
using LinearAlgebra

include("../common/PWGrid_v01.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice.jl")

include("structure_factor.jl")
include("gen_dr.jl")
include("gen_rho.jl")
include("write_xsf.jl")

function test_main()
    Ns = [64,64,64]
    a = 16.0
    LatVecs = gen_lattice_sc(a)
    pwgrid = PWGrid(Ns,LatVecs)

    Npoints = prod(Ns)
    Ω = pwgrid.Ω
    G2 = pwgrid.G2

    # Atomic positions and nuclear charge
    #Xpos = reshape( [5.0, 5.0, 5.0,
    #                 6.5, 5.0, 5.0], (3,2) )
    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    # Calculate structure factor
    @printf("Npoints = %8d\n", Npoints)
    Sf = structure_factor( Xpos, pwgrid.G )

    center = 0.5*a*ones(3)
    println(center)
    dr = gen_dr( pwgrid.r, center )
    rho = gen_rho( Ns, dr, 0.25, Sf )

    sumSf = sum(Sf)
    @printf("sum(Sf) = (%15.10e,%15.10e)\n", real(sumSf), imag(sumSf))
    @printf("integRho = %18.10f\n", sum(rho)*Ω/Npoints)

    write_xsf( "rho.xsf", LatVecs, Xpos )
    write_xsf_3d_crystal( "rho.xsf", Ns, LatVecs, rho )
end

@time test_main()


text:
measures = [accuracy, confmat, misclassification_rate]
@test MLJTuning.measure_adjusted_weights([2, 3, 4], measures) == [-2, 0, 4]
@test MLJTuning.measure_adjusted_weights(nothing, measures) == [-1, 0, 0]
@test_throws(DimensionMismatch,
             MLJTuning.measure_adjusted_weights([2, 3], measures))


text:
function benchmark_cpu_fft(size, dims; FT=Float64, planner_flag=FFTW.PATIENT)
    A = zeros(complex(FT), size...)
    FFT! = FFTW.plan_fft!(A, dims, flags=planner_flag)
    trial = @benchmark ($FFT! * $A) samples=10
    return trial
end

# 3x1D: 61+111+98 = 270 ms
# 3D: 245 ms

function benchmark_cpu_dct(size, dims; FT=Float64, planner_flag=FFTW.PATIENT)
    A = zeros(FT, size...)
    DCT! = FFTW.plan_r2r!(A, FFTW.REDFT10, dims, flags=planner_flag)
    trial = @benchmark ($DCT! * $A) samples=10
    return trial
end

# 3x1D: 86+122+143=351 ms
# 3D: 328 ms

function benchmark_cpu_fft_y(N; FT=Float64, planner_flag=FFTW.PATIENT)
    A = zeros(complex(FT), (N, N, N))
    B = similar(A)
    FFT! = FFTW.plan_fft!(A, 1, flags=planner_flag)

    trial = @benchmark begin
        permutedims!($B, $A, (2, 1, 3))
        $FFT! * $B
        permutedims!($A, $B, (2, 1, 3))
    end

    return trial
end

# 250 ms vs. 111 ms

function benchmark_gpu_fft_y(N; FT=Float64)
    A = zeros(complex(FT), (N, N, N)) |> CuArray
    B = similar(A)
    FFT! = CUDA.CUFFT.plan_fft!(A, 1)

    trial = @benchmark begin
        CUDA.@sync begin
            permutedims!($B, $A, (2, 1, 3))
            $FFT! * $B
            permutedims!($A, $B, (2, 1, 3))
        end
    end

    return trial
end

# 5 ms vs. 3 ms


text:
module RecursiveLeastSquares

using LinearAlgebra
using Distributions

const REALVEC{T<:Real} = Vector{T}
const REALMAT{T<:Real} = Matrix{T}

include("types.jl")
include("kernels.jl")
include("RLS.jl")
include("KRLS.jl")
include("GP.jl")

# export types
export Linear, RBF, RationalQuadratic, Periodic, LocallyPeriodic

# Export functions
export RLS, KRLS, GP

end


text:
module ModflowInterface

using Printf
using Modflow6_jll
import BasicModelInterface as BMI

const BMI_LENCOMPONENTNAME = unsafe_load(cglobal((:BMI_LENCOMPONENTNAME, libmf6), Cint))
const BMI_LENVARADDRESS = unsafe_load(cglobal((:BMI_LENVARADDRESS, libmf6), Cint))
const BMI_LENVARTYPE = unsafe_load(cglobal((:BMI_LENVARTYPE, libmf6), Cint))
const BMI_LENGRIDTYPE = unsafe_load(cglobal((:BMI_LENGRIDTYPE, libmf6), Cint))
const BMI_LENERRMESSAGE = unsafe_load(cglobal((:BMI_LENERRMESSAGE, libmf6), Cint))

struct ModflowModel end

@enum Status begin
    success = 0
    failure = 1
end

@enum State begin
    uninitialized = 1
    initialized = 2
end

# Utilities

function trimmed_string(buffer)::String
    string_end = findfirst(iszero, buffer) - 1
    return String(buffer[1:string_end])
end

# TODO: Maybe check for non-ASCII chars
function prepare_string(s::String)
    return uppercase(s)
end

function parse_type(type::String)::Type
    type = lowercase(type)
    if startswith(type, "double")
        return Float64
    elseif startswith(type, "float")
        return Float32
    elseif startswith(type, "int")
        return Int32
    else
        error("unsupported type")
    end
    return
end

function execute_function(m::ModflowModel, f::Function, args...)
    result = f(args...)
    if result != success
        try
            error_message = get_last_bmi_error(m)
            component_name = get_component_name(m)
            @printf("--- Kernel message (%s) --- \n=> %s", component_name, error_message)
        catch
            println("--- Kernel Message --- \n=> no details ...")
        end
        error("BMI exception")
    end
    return
end

# BMI proper

function BMI.initialize(::Type{ModflowModel})
    @ccall libmf6.initialize()::Cint
    return ModflowModel()
end

function BMI.finalize(::ModflowModel)
    @ccall libmf6.finalize()::Cint
    return
end

function BMI.update(::ModflowModel)
    @ccall libmf6.update()::Cint
    return
end

function BMI.get_start_time(::ModflowModel)::Float64
    start_time = Ref{Float64}(0)
    @ccall libmf6.get_current_time(start_time::Ptr{Float64})::Cint
    return start_time[]
end

function BMI.get_current_time(::ModflowModel)::Float64
    current_time = Ref{Float64}(0)
    @ccall libmf6.get_current_time(current_time::Ptr{Float64})::Cint
    return current_time[]
end

function BMI.get_end_time(::ModflowModel)::Float64
    end_time = Ref{Float64}(0)
    @ccall libmf6.get_end_time(end_time::Ptr{Float64})::Cint
    return end_time[]
end

function BMI.get_time_step(::ModflowModel)::Float64
    time_step = Ref{Float64}(0)
    @ccall libmf6.get_end_time(time_step::Ptr{Float64})::Cint
    return time_step[]
end

"Returns the string `MODFLOW 6`"
function BMI.get_component_name(::ModflowModel)::String
    buffer = zeros(UInt8, BMI_LENCOMPONENTNAME)
    @ccall libmf6.get_component_name(buffer::Ptr{UInt8})::Cint
    return trimmed_string(buffer)
end

function BMI.get_input_item_count(::ModflowModel)::Int
    count = Ref{Cint}(0)
    @ccall libmf6.get_input_item_count(count::Ptr{Cint})::Cint
    return Int(count[])
end

function BMI.get_output_item_count(::ModflowModel)::Int
    count = Ref{Cint}(0)
    @ccall libmf6.get_output_item_count(count::Ptr{Cint})::Cint
    return Int(count[])
end

function BMI.get_input_var_names(m::ModflowModel)::Vector{String}
    shape = (BMI_LENVARADDRESS, BMI.get_input_item_count(m))
    buffer = zeros(UInt8, shape)
    @ccall libmf6.get_input_var_names(buffer::Ptr{UInt8})::Cint
    return [trimmed_string(part) for part in eachcol(buffer)]
end

function BMI.get_output_var_names(m::ModflowModel)::Vector{String}
    shape = (BMI_LENVARADDRESS, BMI.get_output_item_count(m))
    buffer = zeros(UInt8, shape)
    @ccall libmf6.get_output_var_names(buffer::Ptr{UInt8})::Cint
    return [trimmed_string(part) for part in eachcol(buffer)]
end

function BMI.get_var_grid(::ModflowModel, name::String)::Int
    grid_id = Ref{Cint}(0)
    @ccall libmf6.get_var_grid(name::Ptr{UInt8}, grid_id::Ptr{Cint})::Cint
    return Int(grid_id[])
end

function BMI.get_var_type(::ModflowModel, name::String)::String
    buffer = zeros(UInt8, BMI_LENVARTYPE)
    @ccall libmf6.get_var_type(name::Ptr{UInt8}, buffer::Ptr{UInt8})::Cint
    return trimmed_string(buffer)
end

function BMI.get_var_itemsize(::ModflowModel, name::String)::Int
    item_size = Ref{Cint}(0)
    @ccall libmf6.get_var_itemsize(name::Ptr{UInt8}, item_size::Ptr{Cint})::Cint
    return Int(item_size[])
end

function BMI.get_var_nbytes(::ModflowModel, name::String)::Int
    nbytes = Ref{Cint}(0)
    @ccall libmf6.get_var_nbytes(name::Ptr{UInt8}, nbytes::Ptr{Cint})::Cint
    return Int(nbytes[])
end

function BMI.get_value(m::ModflowModel, name::String)
    return copy(BMI.get_value_ptr(m, name))
end

#function BMI.get_value(m::ModflowModel, name::String, dest::Array{T} where T)
#    type = parse_type(BMI.get_var_type(m, name))
#    if type != T
#        error("var type does match destination type")
#    copyto!(BMI.get_value_ptr(m, name), dest)
#    return
#end

function BMI.get_value_ptr(m::ModflowModel, name::String)
    type = parse_type(BMI.get_var_type(m, name))
    shape = get_var_shape(m, name)

    null_pointer = Ref(C_NULL)
    if type == Int32
        @ccall libmf6.get_value_ptr_int(name::Ptr{UInt8}, null_pointer::Ptr{Cvoid})::Cint
    elseif type == Float32
        @ccall libmf6.get_value_ptr_float(name::Ptr{UInt8}, null_pointer::Ptr{Cvoid})::Cint
    elseif type == Float64
        @ccall libmf6.get_value_ptr_double(name::Ptr{UInt8}, null_pointer::Ptr{Cvoid})::Cint
    else
        error("unsupported type")
    end
    typed_pointer = Base.unsafe_convert(Ptr{type}, null_pointer[])

    values = unsafe_wrap(Array, typed_pointer, shape)
    return values
end

function BMI.get_grid_rank(::ModflowModel, grid::Int)::Int
    grid_rank = Ref{Cint}(0)
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_var_grid(grid_rank::Ptr{Cint}, c_grid::Ptr{Cint})::Cint
    return Int(grid_rank[])
end

function BMI.get_grid_size(::ModflowModel, grid::Int)::Int
    grid_size = Ref{Cint}(0)
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_var_grid(grid_size::Ptr{Cint}, c_grid::Ptr{Cint})::Cint
    return Int(grid_size[])
end

function BMI.get_grid_type(::ModflowModel, grid::Int)::String
    buffer = zeros(UInt8, BMI_LENGRIDTYPE)
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_var_type(c_grid::Ptr{Cint}, buffer::Ptr{UInt8})::Cint
    return trimmed_string(buffer)
end

function BMI.get_grid_shape(m::ModflowModel, grid::Int)::Tuple{Int}
    rank = BMI.get_grid_rank(m, grid)
    shape = zeros(Cint, rank)
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_shape(c_grid::Ptr{Cint}, shape::Ptr{Cint})::Cint
    return tuple(shape...)
end

function BMI.get_grid_x(::ModflowModel, grid::Int, x::Vector{Float64})
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_x(c_grid::Ptr{Cint}, x::Ptr{Float64})::Cint
    return
end

function BMI.get_grid_y(::ModflowModel, grid::Int, y::Vector{Float64})
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_y(c_grid::Ptr{Cint}, y::Ptr{Float64})::Cint
    return
end

function BMI.get_grid_z(::ModflowModel, grid::Int, z::Vector{Float64})
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_z(c_grid::Ptr{Cint}, z::Ptr{Float64})::Cint
    return
end

function BMI.get_grid_node_count(::ModflowModel, grid::Int)::Int
    grid_node_count = Ref{Cint}(0)
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_node_count(c_grid::Ptr{Cint}, grid_node_count::Ptr{Cint})::Cint
    return Int(grid_node_count[])
end

function BMI.get_grid_face_count(::ModflowModel, grid::Int)::Int
    grid_face_count = Ref{Cint}(0)
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_face_count(c_grid::Ptr{Cint}, grid_face_count::Ptr{Cint})::Cint
    return Int(grid_face_count[])
end

# Mutates input; return or not?
function BMI.get_grid_face_nodes(::ModflowModel, grid::Int, face_nodes::Vector{Cint})::Int
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_face_nodes(c_grid::Ptr{Cint}, face_nodes::Ptr{Cint})::Cint
    return
end

# Mutates input; return or not?
function BMI.get_grid_nodes_per_face(
    ::ModflowModel,
    grid::Int,
    nodes_per_face::Vector{Cint},
)::Int
    c_grid = Ref{Cint}(grid)
    @ccall libmf6.get_grid_nodes_per_face(
        c_grid::Ptr{Cint},
        nodes_per_face::Ptr{Cint},
    )::Cint
    return
end

# Strictly speaking not BMI

function get_last_bmi_error(::ModflowModel)
    buffer = zeros(UInt8, BMI_LENERRMESSAGE)
    @ccall libmf6.get_last_bmi_error(buffer::Ptr{UInt8})::Cint
    return trimmed_string(buffer)
end

function get_var_rank(::ModflowModel, name::String)
    rank = Ref(Cint(0))
    @ccall libmf6.get_var_rank(name::Ptr{UInt8}, rank::Ptr{Cint})::Cint
    return Integer(rank[])
end

function get_var_shape(m::ModflowModel, name::String)
    rank = get_var_rank(m, name)
    shape = Vector{Int32}(undef, rank)
    @ccall libmf6.get_var_shape(name::Ptr{UInt8}, shape::Ptr{Int32})::Cint
    # The BMI interface returns row major shape; Julia's memory layout is
    # column major, so we flip the shape around.
    return tuple(reverse(shape)...)
end

# XMI

function prepare_time_step(::ModflowModel, dt::Float64)
    timestep = Ref(dt)
    @ccall libmf6.prepare_time_step(timestep::Ptr{Float64})::Cint
    return
end

function do_time_step(::ModflowModel)
    @ccall libmf6.do_time_step()::Cint
    return
end

function finalize_time_step(::ModflowModel)
    @ccall libmf6.finalize_time_step()::Cint
    return
end

function get_subcomponent_count(::ModflowModel)::Int
    count = Ref{Cint}(0)
    @ccall libmf6.get_subcomponent_count(count::Ptr{Cint})::Cint
    return Int(count[])
end

function prepare_solve(::ModflowModel, component_id::Int)
    id = Ref{Cint}(component_id)
    @ccall libmf6.prepare_solve(id::Ptr{Cint})::Cint
    return
end

function solve(::ModflowModel, component_id::Int)::Bool
    id = Ref{Cint}(component_id)
    converged = Ref(false)
    @ccall libmf6.solve(id::Ptr{Cint}, converged::Ptr{Bool})::Cint
    return converged[]
end

function finalize_solve(::ModflowModel, component_id::Int)
    id = Ref{Cint}(component_id)
    @ccall libmf6.finalize_solve(id::Ptr{Cint})::Cint
    return
end

function get_var_address(
    ::ModflowModel,
    var_name,
    component_name;
    subcomponent_name = "",
)::String
    v = prepare_string(var_name)
    c = prepare_string(component_name)
    s = prepare_string(subcomponent_name)
    buffer = zeros(UInt8, BMI_LENVARADDRESS)

    @ccall libmf6.get_var_address(
        c::Ptr{UInt8},
        s::Ptr{UInt8},
        v::Ptr{UInt8},
        buffer::Ptr{UInt8},
    )::Cint
    string_end = findfirst(iszero, buffer) - 1
    return String(buffer[1:string_end])
end

end # module


text:
using SpikeNetOpt
SNO = SpikeNetOpt
using Plots
using Tests

##
# smoke tests code mostly borrowed from
# https://github.com/AStupidBear/SpikingNeuralNetworks.jl/blob/master/test/iz_neuron.jl
##

@testset "IZHI" begin

    RS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    IB = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -55, d = 4))
    CH = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -50, d = 2))
    FS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
    TC1 = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.25, c = -65, d = 0.05))
    TC2 = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.02, b = 0.25, c = -65, d = 0.05))
    RZ = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.1, b = 0.26, c = -65, d = 2))
    LTS = SNN.IZ(; N = 1, param = SNN.IZParameter(; a = 0.1, b = 0.25, c = -65, d = 2))
    P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]

    SNN.monitor(P, [:v])
    SNN.monitor(P, [:fire])

    T = 2second
    for t = 0:T
        for p in [RS, IB, CH, FS, LTS]
            p.I = [10]
        end
        TC1.I = [(t < 0.2T) ? 0mV : 2mV]
        TC2.I = [(t < 0.2T) ? -30mV : 0mV]
        RZ.I = [(0.5T < t < 0.6T) ? 10mV : 0mV]
        SNN.sim!(P, [], 0.1ms)
    end
    for p in P
        spikes = SNO.get_spikes(p)
        spikes = [s * ms for s in spikes]
        nspk = size(spikes)[1]
        v = SNN.vecplot(p, :v)
        v |> display
        @test nspk >= 1

    end
end


text:
using ATA
using ATADash
# cd("where your input files are")

# Before running the app, if you want to use a MILP solver, remember to load it
# (ex: using JuMP, Cbc; ata_app()).
ata_app()
# Navigate with the browser to localhost:8080