text:
using Test
using Distributions
using StatsFuns

using LRMoE

@testset "count with exposures" begin

    # mean(LRMoE.GammaCountExpert(1.0, 1.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1.0))
    # mean(LRMoE.GammaCountExpert(3.0, 1.0))
    # mean(LRMoE.GammaCountExpert(4.0, 1.0))
    # mean(LRMoE.GammaCountExpert(5.0, 1.0))

    # mean(LRMoE.GammaCountExpert(1.0, 2.0))
    # mean(LRMoE.GammaCountExpert(2.0, 2.0))
    # mean(LRMoE.GammaCountExpert(3.0, 2.0))
    # mean(LRMoE.GammaCountExpert(4.0, 2.0))
    # mean(LRMoE.GammaCountExpert(5.0, 2.0))

    # mean(LRMoE.GammaCountExpert(1.0, 3.0))
    # mean(LRMoE.GammaCountExpert(2.0, 3.0))
    # mean(LRMoE.GammaCountExpert(3.0, 3.0))
    # mean(LRMoE.GammaCountExpert(4.0, 3.0))
    # mean(LRMoE.GammaCountExpert(5.0, 3.0))

    # mean(LRMoE.GammaCountExpert(1.0, 1/1.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/2.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/3.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/4.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/5.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/6.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/7.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/8.0))
    # mean(LRMoE.GammaCountExpert(1.0, 1/9.0))

    # var(LRMoE.GammaCountExpert(1.0, 1/1.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/2.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/3.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/4.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/5.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/6.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/7.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/8.0))
    # var(LRMoE.GammaCountExpert(1.0, 1/9.0))

    # mean(LRMoE.GammaCountExpert(2.0, 1/1.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/2.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/3.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/4.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/5.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/6.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/7.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/8.0))
    # mean(LRMoE.GammaCountExpert(2.0, 1/9.0))

    # var(LRMoE.GammaCountExpert(2.0, 1/1.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/2.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/3.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/4.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/5.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/6.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/7.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/8.0))
    # var(LRMoE.GammaCountExpert(2.0, 1/9.0))

    # mean(LRMoE.NegativeBinomialExpert(1.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(2.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(3.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(4.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(5.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(6.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(7.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(8.0, 0.20))
    # mean(LRMoE.NegativeBinomialExpert(9.0, 0.20))

    # var(LRMoE.NegativeBinomialExpert(1.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(2.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(3.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(4.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(5.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(6.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(7.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(8.0, 0.20))
    # var(LRMoE.NegativeBinomialExpert(9.0, 0.20))

end

@testset "expert_ll: continuous, NonZI" begin

    # param1 = 3.0
    # param2 = 2.0
    # param3 = 5.0

    # l = Distributions.LogNormal(param1, param2)
    # x = rand(l, 20)

    # tmp = LRMoE.BurrExpert(param1, param2, param3)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.GammaExpert(param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.InverseGaussianExpert(param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.LogNormalExpert(param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.WeibullExpert(param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

end

@testset "expert_ll: continuous, ZI" begin

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = 3.0
    # param2 = 2.0
    # param3 = 5.0

    # l = Distributions.LogNormal(param1, param2)
    # x = rand(l, 20) .+ 0.05

    # tmp = LRMoE.ZIBurrExpert(param0, param1, param2, param3)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6)

    # tmp = LRMoE.ZIGammaExpert(param0, param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6)

    # tmp = LRMoE.ZIInverseGaussianExpert(param0, param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6)

    # tmp = LRMoE.ZILogNormalExpert(param0, param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6)

    # tmp = LRMoE.ZIWeibullExpert(param0, param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6)

end

@testset "expert_ll: discrete, NonZI" begin

    # param1 = 20
    # param2 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param3 = 5.0

    # l = Distributions.Binomial(param1, param2)
    # x = rand(l, 20)

    # tmp = LRMoE.BinomialExpert(param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.GammaCountExpert(param1, param3)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1]
    # # tmp_expo = LRMoE.GammaCountExpert(param1, param3/expo)

    # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), LRMoE.logpdf.(tmp_expo, x))

    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), LRMoE.logpdf.(tmp_expo, x))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0.0))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0))
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6)

    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x)))
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(LRMoE.logpdf.(tmp_expo, 0)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.NegativeBinomialExpert(param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1]
    # # tmp_expo = LRMoE.NegativeBinomialExpert(param1*expo, param2)

    # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), LRMoE.logpdf.(tmp_expo, x))

    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), LRMoE.logpdf.(tmp_expo, x))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0.0))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0))
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6)

    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x)))
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(LRMoE.logpdf.(tmp_expo, 0)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.PoissonExpert(param1)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1]
    # # tmp_expo = LRMoE.PoissonExpert(param1*expo)

    # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), LRMoE.logpdf.(tmp_expo, x))

    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), LRMoE.logpdf.(tmp_expo, x))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0.0))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0))
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6)

    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x)))
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(LRMoE.logpdf.(tmp_expo, 0)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6)

end

@testset "expert_ll: discrete, ZI" begin

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = 20
    # param2 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param3 = 5.0

    # l = Distributions.Binomial(param1, param2)
    # x = rand(l, 20) .+ 1

    # tmp = LRMoE.ZIBinomialExpert(param0, param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.ZIGammaCountExpert(param0, param1, param3)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1]
    # # tmp_expo = LRMoE.ZIGammaCountExpert(param0, param1, param3/expo)

    # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x))

    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0))))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0))))
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6)

    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x)))
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))), atol = 1e-06)
    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.ZINegativeBinomialExpert(param0, param1, param2)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1]
    # # tmp_expo = LRMoE.ZINegativeBinomialExpert(param0, param1*expo, param2)

    # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x))

    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0))))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0))))
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6)

    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x)))
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))), atol = 1e-06)
    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6)

    # tmp = LRMoE.ZIPoissonExpert(param0, param1)

    # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))

    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6)
    # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6)

    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06)
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0))))
    # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x)))
    # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6)

    # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1]
    # # tmp_expo = LRMoE.ZIPoissonExpert(param0, param1*expo)

    # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x))

    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0))))
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6)
    # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6)

    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06)
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0))))
    # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6)

    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x)))
    # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))), atol = 1e-06)
    # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6)

end

@testset "exposurize_expert" begin

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(10, 20), 1)[1]
    # param2 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param3 = rand(Distributions.Uniform(5, 10), 1)[1]

    # expo = rand(Distributions.Uniform(0.5, 20), 1)[1]

    # tmp = LRMoE.BurrExpert(param1, param2, param3)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)
    # tmp = LRMoE.ZIBurrExpert(param0, param1, param2, param3)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)

    # tmp = LRMoE.GammaExpert(param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)
    # tmp = LRMoE.ZIGammaExpert(param0, param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)

    # tmp = LRMoE.InverseGaussianExpert(param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)
    # tmp = LRMoE.ZIInverseGaussianExpert(param0, param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)

    # tmp = LRMoE.LogNormalExpert(param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)
    # tmp = LRMoE.ZILogNormalExpert(param0, param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)

    # tmp = LRMoE.WeibullExpert(param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)
    # tmp = LRMoE.ZIWeibullExpert(param0, param1, param2)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)

    # tmp = LRMoE.BinomialExpert(20, param0)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)
    # tmp = LRMoE.ZIBinomialExpert(0.20, 20, param0)
    # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo)

    # tmp = LRMoE.GammaCountExpert(param1, param2)
    # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo)
    # @test isapprox([LRMoE.params(tmp)...] .* [1, 1/expo], [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...], atol = 1e-6)
    # tmp = LRMoE.ZIGammaCountExpert(param0, param1, param2)
    # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo)
    # @test isapprox([LRMoE.params(tmp)...] .* [1, 1, 1/expo], [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...], atol = 1e-6)

    # tmp = LRMoE.NegativeBinomialExpert(param1, param0)
    # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo)
    # @test [LRMoE.params(tmp)...] .* [expo, 1] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...]
    # tmp = LRMoE.ZINegativeBinomialExpert(param0, param1, param0)
    # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo)
    # @test [LRMoE.params(tmp)...] .* [1, expo, 1] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...]

    # tmp = LRMoE.PoissonExpert(param1)
    # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo)
    # @test [LRMoE.params(tmp)...] .* [expo] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...]
    # tmp = LRMoE.ZIPoissonExpert(param0, param1)
    # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo)
    # @test [LRMoE.params(tmp)...] .* [1, expo] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...]

    # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3);
    #         BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)]

    # expo = rand(Distributions.Uniform(0.5, 20), 15)

    # @test size(LRMoE.exposurize_model(model, exposure = expo)) == (2, 4, 15)

    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,1,:] == fill(mean(model[1,1]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,2,:] == fill(mean(model[1,2]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,3,:] == fill(mean(model[1,3]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,4,:] == fill(mean(model[1,4]), 15)

    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,1,:] == fill(mean(model[2,1]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,2,:] == mean(model[2,2]) .* expo

    # model_expod = LRMoE.exposurize_model(model, exposure = expo)
    # Y_tmp = hcat(rand(Distributions.LogNormal(2, 1), 15), rand(Distributions.Poisson(5), 15))

    # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[1,:], model_expod[:,:,1])) == (2, 4)
    # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[5,:], model_expod[:,:,5])) == (2, 4)
    # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[12,:], model_expod[:,:,12])) == (2, 4)
    # @test size(LRMoE.expert_ll_list_exact(Y_tmp, model_expod)) == (2, 4, 15)

    # idx = rand(1:15, 1)[1]
    # mat_tmp = LRMoE.expert_ll_ind_mat_exact(Y_tmp[idx,:], model_expod[:,:,idx])
    # @test mat_tmp[1,1] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1])
    # @test mat_tmp[1,2] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1])
    # @test mat_tmp[1,3] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1])
    # @test mat_tmp[1,4] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1])
    # @test mat_tmp[2,1] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2])
    # @test mat_tmp[2,2] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2])
    # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2])
    # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2])

    # idx = rand(1:15, 1)[1]
    # cube_tmp = LRMoE.expert_ll_list_exact(Y_tmp, model_expod)
    # @test cube_tmp[1,1,idx] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1])
    # @test cube_tmp[1,2,idx] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1])
    # @test cube_tmp[1,3,idx] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1])
    # @test cube_tmp[1,4,idx] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1])
    # @test cube_tmp[2,1,idx] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2])
    # @test cube_tmp[2,2,idx] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2])
    # @test cube_tmp[2,3,idx] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2])
    # @test cube_tmp[2,4,idx] == LRMoE.expert_ll_exact(model_expod[2,4,idx], Y_tmp[idx,2])

end

@testset "exposurize_model" begin
    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(10, 20), 1)[1]
    # param2 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param3 = rand(Distributions.Uniform(5, 10), 1)[1]

    # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3);
    #         BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)]

    # expo = rand(Distributions.Uniform(0.5, 20), 15)

    # @test size(LRMoE.exposurize_model(model, exposure = expo)) == (2, 4, 15)

    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,1,:] == fill(mean(model[1,1]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,2,:] == fill(mean(model[1,2]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,3,:] == fill(mean(model[1,3]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,4,:] == fill(mean(model[1,4]), 15)

    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,1,:] == fill(mean(model[2,1]), 15)
    # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,2,:] == mean(model[2,2]) .* expo

    # model_expod = LRMoE.exposurize_model(model, exposure = expo)
    # Y_tmp = hcat(rand(Distributions.LogNormal(2, 1), 15), rand(Distributions.Poisson(5), 15))

    # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[1,:], model_expod[:,:,1])) == (2, 4)
    # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[5,:], model_expod[:,:,5])) == (2, 4)
    # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[12,:], model_expod[:,:,12])) == (2, 4)
    # @test size(LRMoE.expert_ll_list_exact(Y_tmp, model_expod)) == (2, 4, 15)

    # idx = rand(1:15, 1)[1]
    # mat_tmp = LRMoE.expert_ll_ind_mat_exact(Y_tmp[idx,:], model_expod[:,:,idx])
    # @test mat_tmp[1,1] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1])
    # @test mat_tmp[1,2] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1])
    # @test mat_tmp[1,3] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1])
    # @test mat_tmp[1,4] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1])
    # @test mat_tmp[2,1] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2])
    # @test mat_tmp[2,2] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2])
    # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2])
    # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2])

    # idx = rand(1:15, 1)[1]
    # cube_tmp = LRMoE.expert_ll_list_exact(Y_tmp, model_expod)
    # @test cube_tmp[1,1,idx] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1])
    # @test cube_tmp[1,2,idx] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1])
    # @test cube_tmp[1,3,idx] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1])
    # @test cube_tmp[1,4,idx] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1])
    # @test cube_tmp[2,1,idx] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2])
    # @test cube_tmp[2,2,idx] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2])
    # @test cube_tmp[2,3,idx] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2])
    # @test cube_tmp[2,4,idx] == LRMoE.expert_ll_exact(model_expod[2,4,idx], Y_tmp[idx,2])

    # Y_tmp = hcat(0.2 .* Y_tmp[:,1], 0.8 .* Y_tmp[:,1], 1.2 .* Y_tmp[:,1], 5 .* Y_tmp[:,1], fill(0, length(Y_tmp[:,2])), Y_tmp[:,2], Y_tmp[:,2], fill(Inf, length(Y_tmp[:,2])))

    # @test size(LRMoE.expert_ll_ind_mat(Y_tmp[1,:], model_expod[:,:,1])) == (2, 4)
    # @test size(LRMoE.expert_ll_ind_mat(Y_tmp[5,:], model_expod[:,:,5])) == (2, 4)
    # @test size(LRMoE.expert_ll_ind_mat(Y_tmp[12,:], model_expod[:,:,12])) == (2, 4)
    # @test size(LRMoE.expert_ll_list(Y_tmp, model_expod)) == (2, 4, 15)

    # idx = rand(1:15, 1)[1]
    # mat_tmp = LRMoE.expert_ll_ind_mat(Y_tmp[idx,:], model_expod[:,:,idx])
    # @test mat_tmp[1,1] == LRMoE.expert_ll(model_expod[1,1,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test mat_tmp[1,2] == LRMoE.expert_ll(model_expod[1,2,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test mat_tmp[1,3] == LRMoE.expert_ll(model_expod[1,3,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test mat_tmp[1,4] == LRMoE.expert_ll(model_expod[1,4,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test mat_tmp[2,1] == LRMoE.expert_ll(model_expod[2,1,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])
    # @test mat_tmp[2,2] == LRMoE.expert_ll(model_expod[2,2,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])
    # @test mat_tmp[2,3] == LRMoE.expert_ll(model_expod[2,3,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])
    # @test mat_tmp[2,3] == LRMoE.expert_ll(model_expod[2,3,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])

    # idx = rand(1:15, 1)[1]
    # cube_tmp = LRMoE.expert_ll_list(Y_tmp, model_expod)
    # @test cube_tmp[1,1,idx] == LRMoE.expert_ll(model_expod[1,1,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test cube_tmp[1,2,idx] == LRMoE.expert_ll(model_expod[1,2,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test cube_tmp[1,3,idx] == LRMoE.expert_ll(model_expod[1,3,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test cube_tmp[1,4,idx] == LRMoE.expert_ll(model_expod[1,4,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4])
    # @test cube_tmp[2,1,idx] == LRMoE.expert_ll(model_expod[2,1,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])
    # @test cube_tmp[2,2,idx] == LRMoE.expert_ll(model_expod[2,2,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])
    # @test cube_tmp[2,3,idx] == LRMoE.expert_ll(model_expod[2,3,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])
    # @test cube_tmp[2,4,idx] == LRMoE.expert_ll(model_expod[2,4,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8])

    # # ll related
    # gate_tmp = rand(Distributions.Uniform(-1, 1), 15, 4)

    # dim_agg = LRMoE.loglik_aggre_dim(cube_tmp)
    # @test size(dim_agg) == (15, 4)
    # @test size(LRMoE.loglik_aggre_gate_dim(gate_tmp, dim_agg)) == (15, 4)

    # # Simulation related
    # # X_tmp = rand(Distributions.Uniform(-1, 1), 15, 7)
    # # α_tmp = rand(Distributions.Uniform(-1, 1), 4, 7)

    # # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3);
    # #         BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)]

    # # expo = rand(Distributions.Uniform(0.5, 20), 15)

    # # model_expod = LRMoE.exposurize_model(model, exposure = expo)

    # # sim_tmp = LRMoE.sim_dataset(α_tmp, X_tmp, model, exposure = expo)
    # # mean(sim_tmp, dims = 1)

    # # sim_tmp = LRMoE.sim_dataset(α_tmp, X_tmp, model)
    # # mean(sim_tmp, dims = 1)

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
## Influence coefficients for velocities
#============================================#

# Helper
function transform(x, z, x1, x2)
    s1, s2 = x - x1, x - x2
    r1, r2 = √(s1^2 + z^2), √(s2^2 + z^2)
    θ1, θ2 = atan(z, s1), atan(z, s2)

    r1, r2, θ1, θ2
end

term_1(xp, x, z, x1, x2, r1, r2, θ1, θ2) = (x2 - x1) - z * (θ2 - θ1) + (x - xp) * log(r2 / r1)
term_2(xp, x, z, r1, r2, θ1, θ2) = (x - xp) * (θ2 - θ1) + z * log(r2 / r1)

boundary_term(xp, x, x1, x2) = (1/2π, 1/2 * (x - xp) / (x2 - x1))

# Source velocities
#============================================#

function constant_source_velocity(σ, x, z, x1, x2)
    r1, r2, θ1, θ2 = transform(x, z, x1, x2)
    @. σ / 2π * SVector(-log(r2 / r1), θ2 - θ1)
end

function linear_source_velocity_a(σ1, x, z, x1, x2)
    r1, r2, θ1, θ2 = transform(x, z, x1, x2)
    if abs(z) <= 1e-12
        ua, wa = boundary_term(x2, x, x1, x2)
        σ1 .* SVector(ua, -wa)
    else
        ua = term_1(x2, x, z, x1, x2, r1, r2, θ1, θ2)
        wa = term_2(x2, x, z, r1, r2, θ1, θ2)

        @. σ1 / (2π * (x2 - x1)) * SVector(ua, -wa)
    end
end

function linear_source_velocity_b(σ2, x, z, x1, x2)
    r1, r2, θ1, θ2 = transform(x, z, x1, x2)
    if abs(z) <= 1e-12
        ub, wb = boundary_term(x1, x, x1, x2)
        σ2 .* SVector(-ub, wb)
    else
        ub = term_1(x1, x, z, x1, x2, r1, r2, θ1, θ2)
        wb = term_2(x1, x, z, r1, r2, θ1, θ2)

        @. σ2 / (2π * (x2 - x1)) * SVector(-ub, wb)
    end
end

# Vortex velocities
#============================================#

function linear_vortex_velocity_a(γ1, x, z, x1, x2)
    r1, r2, θ1, θ2 = transform(x, z, x1, x2)
    # if z == 0.
    #   ua, wa = (reverse ∘ boundary_term)(x2, x, x1, x2)
    #   γ1 .* SVector(-ua, -wa)
    # else
        ua = term_2(x2, x, z, r1, r2, θ1, θ2)
        wa = term_1(x2, x, z, x1, x2, r1, r2, θ1, θ2)

        @. γ1 / (2π * (x2 - x1)) * SVector(-ua, -wa)
    # end
end

function linear_vortex_velocity_b(γ2, x, z, x1, x2)
    r1, r2, θ1, θ2 = transform(x, z, x1, x2)
    # if z == 0.
    #   ub, wb = (reverse ∘ boundary_term)(x1, x, x1, x2)
    #   γ2 .* SVector(ub, wb)
    # else
        ub = term_2(x1, x, z, r1, r2, θ1, θ2)
        wb = term_1(x1, x, z, x1, x2, r1, r2, θ1, θ2)

        @. γ2 / (2π * (x2 - x1)) * SVector(ub, wb)
    # end
end

ext: jl


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

using HttpServer
using Swifter
import Swifter: QueryResult, ResultInfo, App

# #16107  230e767f6cecbaed426965f8badc5726bb7db72f
using Compat
if VERSION < v"0.5-dev+4194"
    const bytes_to_string = Compat.utf8
else
    const bytes_to_string = Compat.String
end

function handle(color)
    current_color = color
    HttpHandler() do req::Request, res::Response
        if ismatch(r"^/initial", req.resource)
            Response("""{"typ": "any", "value": ""}""")
        elseif ismatch(r"^/query", req.resource)
            param = JSON.parse(bytes_to_string(req.data))
            if "Setter" == param["type"]
                current_color = first(param["rhs"])
            end
            Response("""{"typ": "any", "value": "$current_color"}""")
        else
            Response("""{"typ": "symbol", "value": "Failed"}""")
        end
    end
end

param = Dict("lhs"=>Any[(:symbol,:vc),(:symbol,:view),(:symbol,:backgroundColor)], "type"=>"Getter")
result = @query vc.view.backgroundColor
@test QueryResult(ResultInfo(:symbol, Swifter.RequireToInitial), App(""),"/query",param) == result
@test Swifter.RequireToInitial == result

server_one = Server(handle("UIDeviceRGBColorSpace 0 0 0 1"))
server_two = Server(handle("UIDeviceRGBColorSpace 0 1 0 1"))
@async run(server_one, 8000)
@async run(server_two, 8001)
sleep(0.1)

vc1 = initial("http://localhost:8000")
vc2 = initial("http://localhost:8001")

@test "UIDeviceRGBColorSpace 0 0 0 1" == @query vc1.view.backgroundColor
@test "UIDeviceRGBColorSpace 0 1 0 1" == @query vc2.view.backgroundColor

@query vc1.view.backgroundColor = vc2.view.backgroundColor

@test "UIDeviceRGBColorSpace 0 1 0 1" == @query vc1.view.backgroundColor
@test "UIDeviceRGBColorSpace 0 1 0 1" == @query vc2.view.backgroundColor

try
    close(server_one.http)
    close(server_two.http)
catch e
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AliasPath <: SwaggerModel
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    apiVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiVersions
    pattern::Any # spec type: Union{ Nothing, AliasPattern } # spec name: pattern

    function AliasPath(;path=nothing, apiVersions=nothing, pattern=nothing)
        o = new()
        validate_property(AliasPath, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(AliasPath, Symbol("apiVersions"), apiVersions)
        setfield!(o, Symbol("apiVersions"), apiVersions)
        validate_property(AliasPath, Symbol("pattern"), pattern)
        setfield!(o, Symbol("pattern"), pattern)
        o
    end
end # type AliasPath

const _property_map_AliasPath = Dict{Symbol,Symbol}(Symbol("path")=>Symbol("path"), Symbol("apiVersions")=>Symbol("apiVersions"), Symbol("pattern")=>Symbol("pattern"))
const _property_types_AliasPath = Dict{Symbol,String}(Symbol("path")=>"String", Symbol("apiVersions")=>"Vector{String}", Symbol("pattern")=>"AliasPattern")
Base.propertynames(::Type{ AliasPath }) = collect(keys(_property_map_AliasPath))
Swagger.property_type(::Type{ AliasPath }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AliasPath[name]))}
Swagger.field_name(::Type{ AliasPath }, property_name::Symbol) =  _property_map_AliasPath[property_name]

function check_required(o::AliasPath)
    true
end

function validate_property(::Type{ AliasPath }, name::Symbol, val)
end


text:
(prob::DiffEqBase.AbstractODEProblem{uType,tType,isinplace},
                          Alg::FwdEulerAlg;
                          dt=(prob.tspan[2]-prob.tspan[1])/100,
                          tstops=tType[],
                          kwargs... # ignored kwargs
                          ) where {uType,tType,isinplace}

function solve(prob::FVMPDEProblem,
               U0::Matrix{Float64})

    ODE_Function(u,p,t) = HPDE_∂u∂t(prob,u)
    #ODE_Function(u,p,t) = κ_Scheme(prob,u,1,1)
    #ODE_Function(u,p,t) = FTCS(prob,u)
    ODE_Problem = ODEProblem(ODE_Function,U0,(0,prob.Grid.S[end]));
    temp_u = DifferentialEquations.solve(ODE_Problem,dt=prob.Grid.Δs,adaptive=false)
    return temp_u
end

include("HPDE_options.jl")
include("HPDE_problem.jl")

#initialize

#function solve(prob::HPDEProblem;kwargs...)
    #function ∂u∂t(u,p,t)
#q = prob.func_Q(V0)
function HPDEsolve(prob::HPDEProblem,U0::Matrix{Float64})
    U = copy(U0)

    ODE_Function(u,p,t) = HPDE_∂u∂t(prob,u,0,-1)
    #ODE_Function(u,p,t) = FTCS(prob,u)

    for i=2:length(prob.Grid.S)
        println(i)
        prob.Grid.Δy =0
        ODE_Problem = ODEProblem(ODE_Function,U,(0,prob.Grid.Δs));
        condition(u,t,integrator) = true
        affect!(integrator) = println(integrator.t)
        cb = DiscreteCallback(condition,affect!)
        temp_u = DifferentialEquations.solve(ODE_Problem)#,dt=0.001)#,callback=cb,dt=0.1,adaptive=false)
        #println(size(temp_u.t))
        U = copy(temp_u[:,:,end])
        #U = max.(U,zeros(size(U)))
        prob.Grid.U[:,:,i] = copy(U)
    end
    return prob
end

sol = HPDEsolve(prob, U0)

plot(sol.Grid.X,sol.Grid.U[:,:,end]')

    #end
using PlotlyJS
PP = copy(sol.Grid.U)
PP[2,:,:] = sol.Grid.U[2,:,:] ./ sol.Grid.U[1,:,:]
PP[3,:,:] = sol.Grid.U[3,:,:] ./ sol.Grid.U[1,:,:]

plot(PP[:,:,end]')
#end
using Plots
gr()

function animPlot(PP,i)
    animTime = 10
    fps = 10
    _plotsInterval = convert(Int64, sol.Grid.S[end] / animTime / sol.Grid.Δs / fps)
    anims = @animate for j =1:_plotsInterval:size(PP,3)
        Plots.plot(sol.Grid.X,PP[i,:,j],label="t = $((j-1)*sol.Grid.Δs)")
    end
    gif(anims, "anim_U3.gif", fps = fps)
end
animPlot(PP,1)


text:
using Documenter
# include("../src/SAT.jl")
using SAT

# index is equal to the README for the time being
cp(normpath(@__FILE__, "../../README.md"), normpath(@__FILE__, "../src/index.md"); remove_destination=true)

makedocs(modules=[SAT], doctest = false)

rm(normpath(@__FILE__, "../src/index.md"))

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    repo   = "github.com/CarloLucibello/SAT.jl.git",
    julia  = "0.5"
)


text:
# Autogenerated wrapper script for GnuTLS_jll for x86_64-linux-musl
export libgnutls

using Zlib_jll
using GMP_jll
using Nettle_jll
JLLWrappers.@generate_wrapper_header("GnuTLS")
JLLWrappers.@declare_library_product(libgnutls, "libgnutls.so.30")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, GMP_jll, Nettle_jll)
    JLLWrappers.@init_library_product(
        libgnutls,
        "lib/libgnutls.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module xsimdwrap

using VectorizationBase
using VectorizationBase: REGISTER_SIZE
using SIMDPirates
using LoopVectorization: vectorize_body
using MacroTools: @capture

export @xvectorize

const xsimdmath = joinpath(@__DIR__, "..", "deps", "libxsimdmath.so");

const float_to_float = [
    # "abs", "fabs",
    :exp, :exp2, :exp10, :expm1, :log, :log2, :log10, :log1p,
    :sin, :cos, :tan, :asin, :acos, :atan,
    :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
    :erf, :erfc, :tgamma, :lgamma
]
# const float_to_int = [
#     :ceil, :floor, :trunc, :round, :nearbyint, :rint
# ]
const floatfloat_to_float = [
    # "fmod", "remainder", "min", "max", "fmin", "fmax", "fdim",
    :pow, :hypot, :atan2
]
# const floatfloatfloat_to_float = [
#     :clip
# ]
const float_to_floatfloat = [
    :sincos
]
# const float_to_bool = [
#     :isfinite, :isinf, :isnan
# ]

for f ∈ float_to_float
    for bytecount ∈ (16, 32, 64)
        bytecount > REGISTER_SIZE && continue
        W64 = bytecount>>3
        W32 = bytecount>>2
        fname = QuoteNode(Symbol(:xsimd_,f,:_,bytecount<<3))
        dname = QuoteNode(Symbol(:xsimd_,f,:_,bytecount<<3,:d))
        @eval begin
            function $f(x::NTuple{$W64,Core.VecElement{Float64}})
                ccall(
                    ($dname, xsimdmath), NTuple{$W64,Core.VecElement{Float64}}, (NTuple{$W64,Core.VecElement{Float64}},), x
                )
            end
            function $f(x::NTuple{$W32,Core.VecElement{Float32}})
                ccall(
                    ($fname, xsimdmath), NTuple{$W32,Core.VecElement{Float32}}, (NTuple{$W32,Core.VecElement{Float32}},), x
                )
            end
        end
    end
end
for f ∈ floatfloat_to_float
    for bytecount ∈ (16, 32, 64)
        bytecount > REGISTER_SIZE && continue
        W64 = bytecount>>3
        W32 = bytecount>>2
        fname = QuoteNode(Symbol(:xsimd_,f,:_,bytecount<<3))
        dname = QuoteNode(Symbol(:xsimd_,f,:_,bytecount<<3,:d))
        @eval begin
            function $f(x::NTuple{$W64,Core.VecElement{Float64}},y::NTuple{$W64,Core.VecElement{Float64}})
                ccall(
                    ($dname, xsimdmath), NTuple{$W64,Core.VecElement{Float64}},
                    (NTuple{$W64,Core.VecElement{Float64}},NTuple{$W64,Core.VecElement{Float64}}), x, y
                )
            end
            function $f(x::NTuple{$W32,Core.VecElement{Float32}},y::NTuple{$W32,Core.VecElement{Float32}})
                ccall(
                    ($fname, xsimdmath), NTuple{$W32,Core.VecElement{Float32}},
                    (NTuple{$W32,Core.VecElement{Float32}},NTuple{$W32,Core.VecElement{Float32}}), x, y
                )
            end
        end
    end
end
for f ∈ float_to_floatfloat
    for bytecount ∈ (16, 32, 64)
        bytecount > REGISTER_SIZE && continue
        W64 = bytecount>>3
        W32 = bytecount>>2
        fname = QuoteNode(Symbol(:xsimd_,f,:_,bytecount<<3))
        dname = QuoteNode(Symbol(:xsimd_,f,:_,bytecount<<3,:d))
        @eval begin
            function $f(x::NTuple{$W64,Core.VecElement{Float64}})
                ccall(
                    ($dname, xsimdmath), Tuple{NTuple{$W64,Core.VecElement{Float64}},NTuple{$W64,Core.VecElement{Float64}}},
                    (NTuple{$W64,Core.VecElement{Float64}},), x
                )
            end
            function $f(x::NTuple{$W32,Core.VecElement{Float32}})
                ccall(
                    ($fname, xsimdmath), Tuple{NTuple{$W32,Core.VecElement{Float32}},NTuple{$W32,Core.VecElement{Float32}}},
                    (NTuple{$W32,Core.VecElement{Float32}},), x
                )
            end
        end
    end
end

include("vectorize_loops.jl")

end # module


text:
# test helpers
module CustomTypes

using ColorTypes
using ColorTypes.FixedPointNumbers

export C2, C2A, C4, AC4
export StrangeGray, Cyanotype
export AnaglyphColor, CMYK, ACMYK

struct C2{T <: Real} <: Color{T,2}
    c1::T
    c2::T
end

struct C2A{T} <: ColorAlpha{C2{T},T,3} # T should be <: Real
    c1::T
    c2::T
    alpha::T
end
ColorTypes.coloralpha(::Type{<:C2}) = C2A

struct C4{T <: Real} <: Color{T,4}
    c1::T
    c2::T
    c3::T
    c4::T
end
ColorTypes.eltype_default(::Type{<:C4}) = Int16

struct AC4{T <: Real} <: AlphaColor{C4{T},T,5}
    alpha::T
    c1::T
    c2::T
    c3::T
    c4::T
    AC4{T}(c1::T, c2::T, c3::T, c4::T, alpha::T=oneunit(T)) where {T} = new{T}(alpha, c1, c2, c3, c4)
end
ColorTypes.alphacolor(::Type{<:C4}) = AC4
ColorTypes.eltype_default(::Type{<:AC4}) = Int16
# TODO: The following should be generated automatically
AC4{T}(c1, c2, c3, c4, alpha=1) where {T} = AC4{T}(T(c1), T(c2), T(c3), T(c4), T(alpha))

struct StrangeGray{Something,T <: Integer} <: AbstractGray{Normed{T}}
    val::T
    function StrangeGray{T}(g::Normed{T,f}) where {T,f}
        f == sizeof(T) || error()
        new{:X,T}(reinterpret(g))
    end
end
ColorTypes.gray(g::StrangeGray{X,T}) where {X, T} = reinterpret(Normed{T,sizeof(T)}, g.val)

# non-gray color with a single component
struct Cyanotype{T <: Real} <: Color{T,1}
   value::T
end

function Base.convert(::Type{Cout}, c::C) where {Cout <: AbstractRGB, T, C <: Cyanotype{T}}
    r = max(1 - 1.5 * c.value, 0)
    g = 0.98 - 0.7 * c.value
    b = 0.88 - 0.5 * c.value^2
    Cout(T(r), T(g), T(b))
end

# minimal type for testing 2-component color
struct AnaglyphColor{T} <: Color{T,2} # not `TransparentGray`
    left::T
    right::T
end

# minimal type for testing 4-component color
struct CMYK{T <: Fractional} <: Color{T,4} # not `Transparent3`
    c::T
    m::T
    y::T
    k::T
    CMYK{T}(c::T, m::T, y::T, k::T) where {T} = new{T}(c, m, y, k)
end
# TODO: The following should be generated automatically
CMYK{T}(c, m, y, k) where {T} = CMYK{T}(T(c), T(m), T(y), T(k))

ColorTypes.eltype_default(::Type{<:CMYK}) = N0f8
Base.oneunit(::Type{C}) where {C <: CMYK} = (isconcretetype(C) ? C : C{N0f8})(1, 1, 1, 1)

# minimal type for testing 5-component color
struct ACMYK{T <: Fractional} <: AlphaColor{CMYK{T},T,5}
    alpha::T
    c::T
    m::T
    y::T
    k::T
    ACMYK{T}(c::T, m::T, y::T, k::T, alpha::T=oneunit(T)) where {T} = new{T}(alpha, c, m, y, k)
end
# TODO: The following should be generated automatically
ACMYK{T}(c, m, y, k, alpha=1) where {T} = ACMYK{T}(T(c), T(m), T(y), T(k), T(alpha))
ACMYK(c::T, m::T, y::T, k::T, alpha::T=oneunit(T)) where {T} = ACMYK{T}(c, m, y, k, alpha)
ACMYK{T}(col::CMYK{T}, alpha::T=oneunit(T)) where {T} = ACMYK{T}(col.c, col.m, col.y, col.k, alpha)
end # module


text:
struct Input
    a::Number
    b::Number
end

fx = function(input::Input)
    return input.a + input.b
end

max_stars_repo_name: chlins/fx

max_issues_repo_name: chlins/fx

max_forks_repo_name: chlins/fx


text:
""" IO
Input/output library for Batman.
"""
module DataStructures

import DataFrames: DataFrame
include("rootreader.jl")
include("dataframes.jl")

end # Module DataStructures


text:
using Pkg
using TOML
using BinaryBuilderBase: get_commit_sha

get_url_tree_hash(dep::Dict{String,Any}) =
    (; url=dep["repo-url"], tree_hash=dep["git-tree-sha1"])

function get_commits(dep_old::T, dep_new::T) where {T<:NamedTuple{(:url, :tree_hash), Tuple{String, String}}}
    if dep_old.url == dep_new.url && dep_old.tree_hash != dep_new.tree_hash
        # Same URL but different tree hashes
        old = get_commit_sha(dep_old.url, Base.SHA1(dep_old.tree_hash); verbose=true)
        new = get_commit_sha(dep_new.url, Base.SHA1(dep_new.tree_hash); verbose=true)
        @info "Git commit hashes for $(basename(dep_old.url))" old new
        (; old, new)
    else
        # Same tree hashes and/or maybe different URLs, but that's a weird case,
        # updating the packages shouldn't change the URL.
        nothing
    end
end

function main()
    manifest_path = joinpath(@__DIR__, "Manifest.toml")

    # Get current versions of BinaryBuilder and BinaryBuilderBase
    old_deps = TOML.parsefile(manifest_path)["deps"]
    bb_old = get_url_tree_hash(old_deps["BinaryBuilder"][1])
    bbb_old = get_url_tree_hash(old_deps["BinaryBuilderBase"][1])

    # Make sure the registry is up-to-date
    Pkg.Registry.update()
    # Update the dependencies, but no need to re-precompile the dependencies
    withenv("JULIA_PKG_PRECOMPILE_AUTO"=>0) do
        Pkg.update()
    end

    new_deps = TOML.parsefile(manifest_path)["deps"]
    if old_deps == new_deps
        # Nothing changed, bye bye
        @warn("Dependencies didn't change, exiting...")
        return nothing
    end

    # Get old and new URLs and commit hashes of BinaryBuilder and BinaryBuilderBase.
    bb_new = get_url_tree_hash(new_deps["BinaryBuilder"][1])
    bbb_new = get_url_tree_hash(new_deps["BinaryBuilderBase"][1])
    bb_commits = get_commits(bb_old, bb_new)
    bbb_commits = get_commits(bbb_old, bbb_new)

    # Commit message
    commit_title = "[CI] Update manifest"
    commit_body = ""
    # Highlight in the commit message the changes in BinaryBuilder and
    # BinaryBuilderBase.
    if bb_commits !== nothing || bbb_commits !== nothing
        commit_body *= "Notable changes:\n"
        if bb_commits !== nothing
            commit_body *= "* https://github.com/JuliaPackaging/BinaryBuilder.jl/compare/$(bb_commits.old)...$(bb_commits.new)\n"
        end
        if bbb_commits !== nothing
            commit_body *= "* https://github.com/JuliaPackaging/BinaryBuilderBase.jl/compare/$(bbb_commits.old)...$(bbb_commits.new)\n"
        end
    end

    # Save environment variables for the next steps
    write(ENV["GITHUB_ENV"],
          """
          commit_message<<EOF
          $(commit_title)

          $(commit_body)
          EOF
          commit_title<<EOF
          $(commit_title)
          EOF
          commit_body<<EOF
          $(commit_body)
          EOF
          """)
end

main()

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
## MeasurementError

## Author: Michael Creel <michael@yosemite>
## Created: 2010-03-04
## adapted to Julia 27 Aug. 2017

# this script does a Monte Carlo that shows the inconsistency of OLS estimator
# when there is measurement error of the regressors

using Plots, LinearAlgebra, Econometrics
function main()
# the function to be Monte Carlo'ed
function wrapper(n, sig)
	x = randn(n) # an exogenous regressor
	e = randn(n) # the error term
	ystar = zeros(n)
	# generate the dep var
	for t = 2:n
	  ystar[t,:] = 0.0 + 0.9*ystar[t-1] + 1.0*x[t] + e[t]
	end
    # add measurement error
	y = ystar + sig*randn(n)
	# now do OLS, using the data with meas. error in both dep. var. and regressor
	ylag = lag(y,1)
	data = [y ylag x]
	data = data[2:n,:] # drop first obs, missing due to lag
	y = data[:,1]
	ylag = data[:,2]
	x = data[:,3]
	x = [ones(size(x,1),1) ylag x] # data matrix for ols
	b = x\y
	b = b' - [0.0 0.9 1.0] # subtract true values, so mean should be approx. zero if consistent
end

# do the Monte Carlo
n = 100
sig = 0.0
reps = 1000
bs = zeros(reps,3)
for rep = 1:reps
    bs[rep,:] = wrapper(n, sig)
end
# analyze results
histogram(bs[:,2])
sig != 0.0 ? savefig("ylag_n100.png") : savefig("ylag_n100_no_error.png")
gui()
dstats(bs);
return
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
@kwdef struct FractalTriad{IT<:Integer,FT<:AbstractFloat} <: SpatialSampler
    numpoints::IT = 50
    padding::FT = 0.1
end

function _generate!(ft::FractalTriad, sdm::M) where {M<:AbstractMatrix}
    response = zeros(ft.numpoints, 2)

    return response
end


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

# Date Locales

struct DateLocale
    months::Vector{String}
    months_abbr::Vector{String}
    days_of_week::Vector{String}
    days_of_week_abbr::Vector{String}
    month_value::Dict{String, Int}
    month_abbr_value::Dict{String, Int}
    day_of_week_value::Dict{String, Int}
    day_of_week_abbr_value::Dict{String, Int}
end

function locale_dict(names::Vector{<:AbstractString})
    result = Dict{String, Int}()

    # Keep both the common case-sensitive version of the name and an all lowercase
    # version for case-insensitive matches. Storing both allows us to avoid using the
    # lowercase function during parsing.
    for i in 1:length(names)
        name = names[i]
        result[name] = i
        result[lowercase(name)] = i
    end
    return result
end

"""
    DateLocale(["January", "February",...], ["Jan", "Feb",...],
               ["Monday", "Tuesday",...], ["Mon", "Tue",...])

Create a locale for parsing or printing textual month names.

Arguments:

- `months::Vector`: 12 month names
- `months_abbr::Vector`: 12 abbreviated month names
- `days_of_week::Vector`: 7 days of week
- `days_of_week_abbr::Vector`: 7 days of week abbreviated

This object is passed as the last argument to `tryparsenext` and `format`
methods defined for each `AbstractDateToken` type.
"""
function DateLocale(months::Vector, months_abbr::Vector,
                    days_of_week::Vector, days_of_week_abbr::Vector)
    DateLocale(
        months, months_abbr, days_of_week, days_of_week_abbr,
        locale_dict(months), locale_dict(months_abbr),
        locale_dict(days_of_week), locale_dict(days_of_week_abbr),
    )
end

const ENGLISH = DateLocale(
    ["January", "February", "March", "April", "May", "June",
     "July", "August", "September", "October", "November", "December"],
    ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
     "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"],
    ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"],
    ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
)

const LOCALES = Dict{String, DateLocale}("english" => ENGLISH)

for (fn, field) in zip(
    [:dayname_to_value, :dayabbr_to_value, :monthname_to_value, :monthabbr_to_value],
    [:day_of_week_value, :day_of_week_abbr_value, :month_value, :month_abbr_value],
)
    @eval @inline function $fn(word::AbstractString, locale::DateLocale)
        # Maximize performance by attempting to avoid the use of `lowercase` and trying
        # a case-sensitive lookup first
        value = get(locale.$field, word, 0)
        if value == 0
            value = get(locale.$field, lowercase(word), 0)
        end
        value
    end
end

# Date functions

### Core query functions

# Monday = 1....Sunday = 7
dayofweek(days) = mod1(days, 7)

# Number of days in year
"""
    daysinyear(dt::TimeType) -> Int

Return 366 if the year of `dt` is a leap year, otherwise return 365.

# Examples
```jldoctest
julia> Dates.daysinyear(1999)
365

julia> Dates.daysinyear(2000)
366
```
"""
daysinyear(y) = 365 + isleapyear(y)

# Day of the year
const MONTHDAYS = (0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334)
dayofyear(y, m, d) = MONTHDAYS[m] + d + (m > 2 && isleapyear(y))

### Days of the Week
"""
    dayofweek(dt::TimeType) -> Int64

Return the day of the week as an [`Int64`](@ref) with `1 = Monday, 2 = Tuesday, etc.`.

# Examples
```jldoctest
julia> Dates.dayofweek(Date("2000-01-01"))
6
```
"""
dayofweek(dt::TimeType) = dayofweek(days(dt))

const Monday, Tuesday, Wednesday, Thursday, Friday, Saturday, Sunday = 1, 2, 3, 4, 5, 6, 7
const Mon, Tue, Wed, Thu, Fri, Sat, Sun = 1, 2, 3, 4, 5, 6, 7

dayname(day::Integer, locale::DateLocale) = locale.days_of_week[day]
dayabbr(day::Integer, locale::DateLocale) = locale.days_of_week_abbr[day]
dayname(day::Integer; locale::AbstractString="english") = dayname(day, LOCALES[locale])
dayabbr(day::Integer; locale::AbstractString="english") = dayabbr(day, LOCALES[locale])

"""
    dayname(dt::TimeType; locale="english") -> String

Return the full day name corresponding to the day of the week of the `Date` or `DateTime` in
the given `locale`.

# Examples
```jldoctest
julia> Dates.dayname(Date("2000-01-01"))
"Saturday"
```
"""
function dayname(dt::TimeType;locale::AbstractString="english")
    dayname(dayofweek(dt); locale=locale)
end

"""
    dayabbr(dt::TimeType; locale="english") -> String

Return the abbreviated name corresponding to the day of the week of the `Date` or `DateTime`
in the given `locale`.

# Examples
```jldoctest
julia> Dates.dayabbr(Date("2000-01-01"))
"Sat"
```
"""
function dayabbr(dt::TimeType;locale::AbstractString="english")
    dayabbr(dayofweek(dt); locale=locale)
end

# Convenience methods for each day
ismonday(dt::TimeType) = dayofweek(dt) == Mon
istuesday(dt::TimeType) = dayofweek(dt) == Tue
iswednesday(dt::TimeType) = dayofweek(dt) == Wed
isthursday(dt::TimeType) = dayofweek(dt) == Thu
isfriday(dt::TimeType) = dayofweek(dt) == Fri
issaturday(dt::TimeType) = dayofweek(dt) == Sat
issunday(dt::TimeType) = dayofweek(dt) == Sun

# i.e. 1st Monday? 2nd Monday? 3rd Wednesday? 5th Sunday?
"""
    dayofweekofmonth(dt::TimeType) -> Int

For the day of week of `dt`, return which number it is in `dt`'s month. So if the day of
the week of `dt` is Monday, then `1 = First Monday of the month, 2 = Second Monday of the
month, etc.` In the range 1:5.

# Examples
```jldoctest
julia> Dates.dayofweekofmonth(Date("2000-02-01"))
1

julia> Dates.dayofweekofmonth(Date("2000-02-08"))
2

julia> Dates.dayofweekofmonth(Date("2000-02-15"))
3
```
"""
function dayofweekofmonth(dt::TimeType)
    d = day(dt)
    return d < 8 ? 1 : d < 15 ? 2 : d < 22 ? 3 : d < 29 ? 4 : 5
end

# Total number of a day of week in the month
# e.g. are there 4 or 5 Mondays in this month?
const TWENTYNINE = BitSet([1, 8, 15, 22, 29])
const THIRTY = BitSet([1, 2, 8, 9, 15, 16, 22, 23, 29, 30])
const THIRTYONE = BitSet([1, 2, 3, 8, 9, 10, 15, 16, 17, 22, 23, 24, 29, 30, 31])

"""
    daysofweekinmonth(dt::TimeType) -> Int

For the day of week of `dt`, return the total number of that day of the week in `dt`'s
month. Returns 4 or 5. Useful in temporal expressions for specifying the last day of a week
in a month by including `dayofweekofmonth(dt) == daysofweekinmonth(dt)` in the adjuster
function.

# Examples
```jldoctest
julia> Dates.daysofweekinmonth(Date("2005-01-01"))
5

julia> Dates.daysofweekinmonth(Date("2005-01-04"))
4
```
"""
function daysofweekinmonth(dt::TimeType)
    y, m, d = yearmonthday(dt)
    ld = daysinmonth(y, m)
    return ld == 28 ? 4 : ld == 29 ? ((d in TWENTYNINE) ? 5 : 4) :
           ld == 30 ? ((d in THIRTY) ? 5 : 4) :
           (d in THIRTYONE) ? 5 : 4
end

### Months
const January, February, March, April, May, June = 1, 2, 3, 4, 5, 6
const July, August, September, October, November, December = 7, 8, 9, 10, 11, 12
const Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12

monthname(month::Integer, locale::DateLocale) = locale.months[month]
monthabbr(month::Integer, locale::DateLocale) = locale.months_abbr[month]
monthname(month::Integer; locale::AbstractString="english") = monthname(month, LOCALES[locale])
monthabbr(month::Integer; locale::AbstractString="english") = monthabbr(month, LOCALES[locale])

"""
    monthname(dt::TimeType; locale="english") -> String

Return the full name of the month of the `Date` or `DateTime` in the given `locale`.

# Examples
```jldoctest
julia> Dates.monthname(Date("2005-01-04"))
"January"
```
"""
function monthname(dt::TimeType; locale::AbstractString="english")
    monthname(month(dt); locale=locale)
end

"""
    monthabbr(dt::TimeType; locale="english") -> String

Return the abbreviated month name of the `Date` or `DateTime` in the given `locale`.

# Examples
```jldoctest
julia> Dates.monthabbr(Date("2005-01-04"))
"Jan"
```
"""
function monthabbr(dt::TimeType; locale::AbstractString="english")
    monthabbr(month(dt); locale=locale)
end

"""
    daysinmonth(dt::TimeType) -> Int

Return the number of days in the month of `dt`. Value will be 28, 29, 30, or 31.

# Examples
```jldoctest
julia> Dates.daysinmonth(Date("2000-01"))
31

julia> Dates.daysinmonth(Date("2001-02"))
28

julia> Dates.daysinmonth(Date("2000-02"))
29
```
"""
daysinmonth(dt::TimeType) = ((y, m) = yearmonth(dt); return daysinmonth(y, m))

### Years
"""
    isleapyear(dt::TimeType) -> Bool

Return `true` if the year of `dt` is a leap year.

# Examples
```jldoctest
julia> Dates.isleapyear(Date("2004"))
true

julia> Dates.isleapyear(Date("2005"))
false
```
"""
isleapyear(dt::TimeType) = isleapyear(year(dt))

"""
    dayofyear(dt::TimeType) -> Int

Return the day of the year for `dt` with January 1st being day 1.
"""
dayofyear(dt::TimeType) = ((y, m, d) = yearmonthday(dt); return dayofyear(y, m, d))

daysinyear(dt::TimeType) = 365 + isleapyear(dt)

### Quarters
"""
    quarterofyear(dt::TimeType) -> Int

Return the quarter that `dt` resides in. Range of value is 1:4.
"""
function quarterofyear(dt::TimeType)
    m = month(dt)
    return m < 4 ? 1 : m < 7 ? 2 : m < 10 ? 3 : 4
end
const QUARTERDAYS = (0, 90, 181, 273)

"""
    dayofquarter(dt::TimeType) -> Int

Return the day of the current quarter of `dt`. Range of value is 1:92.
"""
dayofquarter(dt::TimeType) = dayofyear(dt) - QUARTERDAYS[quarterofyear(dt)]

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
kinds = Dict("module"=>"module",
			 "variable"=>"variable",
			 "function"=>"function",
			 "datatype"=>"datatype",
			 "union"=>"union",
			 "type" =>"type"
			 )

immutable Tag
	name::AbstractString
	file::AbstractString
	address::AbstractString
	fields::Dict{AbstractString,AbstractString}
end

const mtime_line_prefix = "!_TAG_SOURCE_MTIME"
function write_header(fp::IO, file_mtimes=[])
	println(fp,"!_TAG_FILE_FORMAT   2   /extended format; --format=1 will not append ;\" to lines/")
	println(fp,"!_TAG_FILE_SORTED   0   /0=unsorted, 1=sorted, 2=foldcase/")
	println(fp,"!_TAG_PROGRAM_AUTHOR    Lyndon White aka Frames aka oxinabox /lyndon.white@reseach.uwa.edu.au")
	println(fp,"!_TAG_PROGRAM_NAME  jltags //")
	#println(fp,"!_TAG_PROGRAM_URL   http://ctags.sourceforge.net    /official site/")
	#println(fp,"!_TAG_PROGRAM_VERSION   5.9~svn20110310 //")
	for (filename, mtime) in file_mtimes
		println(fp,mtime_line_prefix,"\t",filename,"\t",string(mtime))
	end
end

"""http://ctags.sourceforge.net/FORMAT
{tagname}<Tab>{tagfile}<Tab>{tagaddress}[;"<Tab>{tagfield}..]

   {tagname}	Any identifier, not containing white space..
   <Tab>	Exactly one TAB character (although many versions of Vi can
		handle any amount of white space).
   {tagfile}	The name of the file where {tagname} is defined, relative to
   		the current directory (or location of the tags file?).
   {tagaddress}	Any Ex command.  When executed, it behaves like 'magic' was
		not set.  It may be restricted to a line number or a search
		pattern (Posix).
Optionally:
   ;"		semicolon + doublequote: Ends the tagaddress in way that looks
		like the start of a comment to Vi.
   {tagfield}	See below.

"""
function write_tag(fp::IO, tag::Tag)
	tag.fields["language"]="julia"

	print(fp, "$(tag.name)\t$(tag.file)\t$(tag.address);\"")
	for (field,value) in tag.fields
		value = escape_string(value)
		print(fp,"\t$field:$value")
	end
	println(fp)
end

###############################

function module_to_filename(mm::Module)
	name = mm==Base ? "sysimg" : mm==Core ? "coreimg": string(mm)
	path = Base.find_source_file(name*".jl")
	path!=nothing ? path : "<module $name: unkown path>"
end

function docs(func::Function, method::Method)
	try
		signiture_types = method.sig.types[2:end]
		doc = Base.Docs.doc(func,signiture_types...)
		return doc.meta[:results][1].text[1] #What a long path
	catch ee
		return "No docs found for " * string(method)
	end
end

function tags(name_sym::Symbol, mm::Module,  value::Function)
	name=string(name_sym)
	Task() do
		for func in Base.MethodList(methods(value))
			if mm!=func.module
				continue #Skip things from other modules (including submodules)
			end
			file, linenum = functionloc(func)
			address = string(linenum)
			args = map(string, func.sig.parameters[2:end])
			fields = Dict{AbstractString,AbstractString}(
				"module" => mm|>string,
				"kind" => kinds["function"],
				"arity" => string(length(args)),
				"args" => join(args,","),
				"string" => string(func),
				"doc" => docs(value, func),
				)
			produce(Tag(name,file,address,fields))
		end
	end
end

function tags(name_sym::Symbol, mm::Module,  value::Union)
	name=string(name_sym)
	Task() do
		filename=module_to_filename(mm)
			#This is a bad attempt But it is what we got right now

		fields=Dict("kind"  => kinds["union"],
					"module"=> string(mm),
					"types"  => join(",",map(string,value.types)),
					"string" => string(value)
					)
		produce(Tag(name,filename,"/name\\s*=/",fields))
	end
end

function tags(name_sym::Symbol,mm::Module, tt::DataType, record_truename=false)
	name=string(name_sym)
	Task() do
		filename = module_to_filename(mm)
		#This is a bad attempt But it is what we got right now
		fields=Dict("kind"  => kinds["datatype"],
					"module"=> string(mm),
					"other_type_parameters"=>join(",",map(string,
														  setdiff(tt.parameters,tt.types))),
					"abstract"=> tt.abstract |> string,
					"mutable"=> tt.mutable |> string,
					"size"=> tt.size |> string,
					"supertype"=>supertype(tt) |> string,
					"bitstype"=>isbits(tt) |> string,
					"string" => string(tt)
					)
		if record_truename
			fields["true_name"]=func.name.name |> string
		end

		if tt.name.names |> length >= nfields(tt) #HACK: Avoid Bug where there are not as many names as fields.
			if name_sym != :Module #HACK/BUG: Module can not have its fieldtype read.
				type_fields = ["$fname::"*string(fieldtype(tt,fname)) for fname in fieldnames(tt)]
				fields["fields"] = join(type_fields,",")
			end
		end

		produce(Tag(name,filename, "/^[typealias|type|abstract|immutable]\\s+name\\s*=/", fields))
	end
end

function tags(name_sym::Symbol, mm::Module,  value::Type)
	Task() do
		for tag in 	tags(name_sym,mm,value.body)
			#TODO: Actually undertand and  special case this fully as a type alias? idk check out `typeof(Vector)`
			tag.fields["kind"] = kinds["type"]
			tag.fields["inner_string"]=tag.fields["string"]
			tag.fields["string"] = string(value)
			produce(tag)
		end
	end
end

function tags(name_sym::Symbol,mm::Module, submodule::Module)
	name=string(name_sym)
	Task() do

		filename = module_to_filename(submodule)
		fields=Dict("kind"  => kinds["module"],
					"module"=> module_parent(submodule) |> string,
					"string"=> string(submodule)
					)
		produce(Tag(name,filename, "/^module name/",fields))
	end
end

function tags(name_sym::Symbol,mm::Module, variable::Any)
	name=string(name_sym)
	Task() do
		filename=module_to_filename(mm)
			#This is a bad attempt But it is what we got right now
			#Accessing variables from modules is a bad idea anyway
		fields=Dict("kind"  => kinds["variable"],
					"module"=> string(mm),
					"type"  => string(typeof(variable)),
					"constant" => string(isconst(mm,name_sym)),
					"string" => string(variable)
					)
		produce(Tag(name,filename, "/name\\s*=/",fields))
	end
end

###############

function tags_from_module(mm::Module)
	try
		println(STDERR, "++++++++ Generating Tags from  $mm  +++++++")
		Task() do
			for name_sym in names(mm)
				value = eval(mm,name_sym)
				map(produce, tags(name_sym, mm,  value))
			end
		end |> collect
	catch ee
		warn(string(Module)*" tagging failed. As "* string(ee))
		[]
	end
end

function name2module(name)
	eval(parse("import "*name))
    eval(parse(name))
end

##################################
#Create Module cache
const cache_location = joinpath(Pkg.Dir.path(),"../jltags_cache/")

function create_module_tagsfile(cache_name, module_name)
	mkpath(cache_location)
	mod = name2module(module_name)
	if mod==nothing
		warn("$(module_name) Not Found")
		return nothing
	end
	open(cache_name,"w") do fp
		tags = tags_from_module(mod) |> collect
		file_mtimes = Dict([tag.file=>stat(tag.file).mtime|>Dates.unix2datetime for tag in tags])
		write_header(fp, file_mtimes)
		for tag in tags
			write_tag(fp,tag)
		end
	end
	#TODO: Consider Setting mtime of the cache to the be mtime of the most recently modified file. Would that be useful?
end

function is_cache_current(cache_name)
	#It is practically impossible to add or remove a tag-able item from a module,
	#without changing one of the existing files (Julia does not do directory based submodules)
	for line in eachline(cache_name)
		if line[1:length(mtime_line_prefix)] == mtime_line_prefix
			prefix, filename, mtime_str = split(line,"\t")
			cached_mtime = DateTime(mtime_str)
			cur_mtime = stat(filename).mtime |> Dates.unix2datetime
			if cached_mtime!=cur_mtime #Forward or back we don't care. ClockSkew etc.
				return false
			end
		end
	end
	return true
end

function get_module_tag_file(module_name)
	println(STDERR, "Sourcing $module_name")
	cache_name = joinpath(cache_location,module_name*".tags")
	if !isfile(cache_name) || !is_cache_current(cache_name)
		create_module_tagsfile(cache_name, module_name)
	end
	cache_name
end

####################################################

function parseall(str)
    #Inspired by https://github.com/JuliaLang/julia/blob/9caa28d50b090cea746ce5c9ad371f5cc2c32644/test/parse.jl#L3
    Task() do
        pos = start(str)
        while !done(str, pos)
            ex, pos = parse(str, pos)
            produce(ex)
        end
    end
end

function get_modules_names(filename)
	code = readstring(filename)
    function inner(::Any)
    end
    function inner(ee::Expr)
         if ee.head == :using || ee.head==:importall
            module_name = join(map(string,ee.args),".")
            produce(module_name)
        else
            for arg in ee.args
                inner(arg)
            end
        end
    end

    Task() do
		produce("Core")
		produce("Base")
        for ee in code |> parseall
            for mod in @task inner(ee)
                produce(mod)
            end
        end
    end
end

######################

function make_standalone_tagfile(filename)
	println(STDERR,"Standalone Tagging $filename")
	mod_names = get_modules_names(filenames)

	open(joinpath(dirname(filename), "."*basename(filename)*".tags"),"w") do fp
		for mod_name in mod_names
			tag_file = get_module_tag_file(mod_name)
			println(fp, open(readstring, tag_file,"r")) #Append it all
		end
	end
end

function print_refered_tagfile_paths(filename)
	mod_names = get_modules_names(filename)
	for mod_name in mod_names
		tag_file = get_module_tag_file(mod_name)
		println(tag_file)
	end
end

###ARGS mode
if length(ARGS)>0
	mode = lowercase(ARGS[1])
	if mode == "standalone"
		make_standalone_tagfile(ARGS[2])
	elseif mode == "refer"
		print_refered_tagfile_paths(ARGS[2])
	else
		println("""Invalid mode selected.
		Usage: jltags mode target`

		jltags primary use is generating tagfiles for modules.
		It always does this, no matter the mode.
		They are stored in $cache_location
		jltags automatically checks if they need to be regenerated each time they are used.

		modes:

			- standalone:
				Generate a standalone tagfile, with the naming scheme ".<target>.tag".
				This file contained all the tags for all modules used in the target file.

			- refer:
				outputs to standard out, a linebreak seperated list of paths to tagfiles for modules that are used by the target

				"""
		)
	end

end

max_stars_repo_name: oxinabox/julia-vim-omnicomplete

max_issues_repo_name: oxinabox/julia-vim-omnicomplete

max_forks_repo_name: oxinabox/julia-vim-omnicomplete