text:
using Test
using Distributions
using StatsFuns
# using Random

μ = 1
σ = 2
# using LRMoE

# Random.seed!(1234)

# tmp = LRMoESTD([1 2 ; 3 4], ["a" "b"])
# tmp1 = LRMoESTDFit(tmp, true, 20, 1.1, 1.2, 1.3)
# summary(tmp1)

using Random, LRMoE, Statistics, Clustering, Distributions

@testset "fitting" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)
    #     l = Distributions.LogNormal(μ, σ)
    #     y = rand(l, 200)

    #     # LogNormal
    #     r = LRMoE.LogNormalExpert(μ, σ)
    #     Y = hcat(fill(0, length(y)), y, y, fill(Inf, length(y)), fill(0, length(y)), 0.80.*y, 1.25.*y, fill(Inf, length(y)))

    #     model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ);
    #             LRMoE.LogNormalExpert(1.2*μ, σ) LRMoE.LogNormalExpert(0.8*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 0.5*σ)]

    #     pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf], [Inf 1.0 Inf]],
    #             [[Inf 1.0 Inf], [Inf 1.0 Inf], [Inf 1.0 Inf]]]

    #     X = rand(Uniform(-1, 1), 200, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0

    #     result = fit_main(Y, X, α_true, model, penalty = false, pen_params = pen_params)

    #     model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ);
    #             LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.80, 1.5*μ, 2.0*σ)]

    #     result = fit_main(Y, X, α_true, model, penalty = false, pen_params = pen_params)

    # end

end

@testset "fitting simulated: lognormal" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 3, 5)
        # α_true[3, :] .= 0.0
        # model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ);
        #          LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.80, 1.5*μ, 2.0*σ)]

        # pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf], [Inf 1.0 Inf]],
        #               [[Inf 1.0 Inf], [Inf 1.0 Inf], [Inf 1.0 Inf]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0
        # model_guess = [LRMoE.LogNormalExpert(0.8*μ, 1.2*σ) LRMoE.LogNormalExpert(μ, 0.9*σ) LRMoE.LogNormalExpert(1.0*μ, 2.5*σ);
        #                LRMoE.ZILogNormalExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.LogNormalExpert(0.75*μ, 0.3*σ) LRMoE.ZILogNormalExpert(0.50, 1.75*μ, 1.0*σ)]

        # # Exact observation
        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_LRMoE(Y, X, α_guess, model_guess, exact_Y = false, penalty = false, pen_α = 10.0, pen_params = pen_params)
        # result = fit_LRMoE(Y_sim, X, α_guess, model_guess, exact_Y = true, penalty = false, pen_α = 10.0, pen_params = pen_params)

        # result = fit_LRMoE(Y, X, α_guess, model_guess, exact_Y = false, penalty = false)
        # result = fit_LRMoE(Y, X, α_guess, model_guess, exact_Y = false, penalty = true)
        # result = fit_LRMoE(Y, X, α_guess, model_guess, exact_Y = false, penalty = true, pen_α = 20.0, pen_params = pen_params)
        # result = fit_LRMoE(Y_sim, X, α_guess, model_guess, exact_Y = true, penalty = true)

        # # With censoring
        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), 0.75.*Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # With truncation
        # Y = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # Y = hcat( 0.80.*Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], 2.50.*Y_sim[:,2])
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], 2.0 .*Y_sim[:,1], fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ecm_iter_max = 20)

        # gate_tmp = LogitGating(α_guess, X)
        # ll_tmp = loglik_np(Y, gate_tmp, model_guess)
        # hcat(Y, EM_E_k(ll_tmp.gate_expert_tn_bar_k))

        # model = [LRMoE.ZILogNormalExpert(0.3, μ, σ) LRMoE.ZILogNormalExpert(0.2, 0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.10, 1.5*μ, 2.0*σ);
        #          LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.ZILogNormalExpert(0.5, 0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.80, 1.5*μ, 2.0*σ)]

        # model_guess = [LRMoE.ZILogNormalExpert(0.50, 0.8*μ, 1.2*σ) LRMoE.ZILogNormalExpert(0.50, μ, 0.9*σ) LRMoE.ZILogNormalExpert(0.50, 1.0*μ, 2.5*σ);
        #                LRMoE.ZILogNormalExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.ZILogNormalExpert(0.50, 0.75*μ, 0.3*σ) LRMoE.ZILogNormalExpert(0.50, 1.75*μ, 1.0*σ)]

        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], 2.0 .*Y_sim[:,1], fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ecm_iter_max = 20)

        # model_guess = [LRMoE.ZILogNormalExpert(0.50, 0.8*μ, 1.2*σ) LRMoE.ZILogNormalExpert(0.50, μ, 0.9*σ) LRMoE.ZILogNormalExpert(0.50, 1.0*μ, 2.5*σ);
        #                LRMoE.ZILogNormalExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.ZILogNormalExpert(0.50, 0.75*μ, 0.3*σ) LRMoE.ZILogNormalExpert(0.50, 1.75*μ, 1.0*σ)]

        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], 2.0 .*Y_sim[:,1], fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ecm_iter_max = 20)

        # gate_tmp = LogitGating(α_true, X)
        # ll_tmp = loglik_np(Y, gate_tmp, model_guess)

        # hcat(Y, expm1.( - ll_tmp.gate_expert_tn ), exp.(ll_tmp.gate_expert_tn_bar_comp) ./ exp.(ll_tmp.gate_expert_tn_bar))

        # hcat(Y, expm1.( - ll_tmp.gate_expert_tn ), EM_E_z_lat(ll_tmp.gate_expert_tn_bar_comp, ll_tmp.gate_expert_tn_bar))

    # end

end

@testset "fitting simulated: poisson" begin

    # d = Distributions.Poisson(10)

    # λλ = rand(d, 1)

    # for λ in λλ

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 3, 5)
        # α_true[3, :] .= 0.0
        # model = [LRMoE.PoissonExpert(λ) LRMoE.PoissonExpert(0.5*λ) LRMoE.PoissonExpert(1.5*λ);
        #          LRMoE.ZIPoissonExpert(0.4, λ) LRMoE.PoissonExpert(1.2*λ) LRMoE.ZIPoissonExpert(0.80, 2.0*λ)]

        # pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
        #               [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

        # pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
        #               [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0
        # model_guess = [LRMoE.PoissonExpert(0.8*λ) LRMoE.PoissonExpert(λ) LRMoE.PoissonExpert(1.0*λ);
        #                LRMoE.ZIPoissonExpert(0.50, 2.0*λ) LRMoE.PoissonExpert(0.75*λ) LRMoE.ZIPoissonExpert(0.50, 1.75*λ)]

        # # Exact observation
        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # # With censoring
        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), floor.(0.75.*Y_sim[:,1]), fill(Inf, length(Y_sim[:,1])), fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # # With truncation
        # Y = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), floor.(0.30.*Y_sim[:,2]), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)
        # Y = hcat( floor.(0.80.*Y_sim[:,1]), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), floor.(0.30.*Y_sim[:,2]), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # ###########################################################
        # #! THERE IS ERROR IN THE FIRST OF THE FOLLOWING THREE CASES
        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), ceil.(2.50.*Y_sim[:,2]))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), ceil.(2.50.*Y_sim[:,2]).+1)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), floor.(0.30.*Y_sim[:,2]), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), ceil.(2.50.*Y_sim[:,2]))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

    # end

end

@testset "fitting simulated: binomial" begin

    # d = Distributions.Poisson(10)

    # λλ = rand(d, 1)

    # for λ in λλ

    #     n = 30
    #     p = 0.2

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.BinomialExpert(n, p) LRMoE.BinomialExpert(n, 0.5*p) LRMoE.BinomialExpert(n, 1.3*p);
    #              LRMoE.ZIBinomialExpert(0.4, n, p) LRMoE.BinomialExpert(n, 1.2*p) LRMoE.ZIBinomialExpert(0.70, n, 2.0*p)]

    #     pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
    #                   [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

    #     pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
    #                   [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0

    #     model_guess = [LRMoE.BinomialExpert(n, 0.8*p) LRMoE.BinomialExpert(n, p) LRMoE.BinomialExpert(n, 1.0*p);
    #              LRMoE.ZIBinomialExpert(0.5, n, 2.0*p) LRMoE.BinomialExpert(n, 0.75*p) LRMoE.ZIBinomialExpert(0.50, n, 1.75*p)]

    #     # Exact observation
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     # With censoring
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), floor.(0.75.*Y_sim[:,1]), fill(Inf, length(Y_sim[:,1])), fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     # With truncation
    #     Y = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), floor.(0.30.*Y_sim[:,2]), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)
    #     Y = hcat( floor.(0.80.*Y_sim[:,1]), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), floor.(0.30.*Y_sim[:,2]), floor.(0.80.*Y_sim[:,2]), ceil.(1.20.*Y_sim[:,2]), fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    # end

end

@testset "fitting simulated: gamma" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

        # μ = 2
        # σ = 3

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 3, 5)
        # α_true[3, :] .= 0.0
        # model = [LRMoE.GammaExpert(μ, σ) LRMoE.GammaExpert(0.5*μ, 0.6*σ) LRMoE.GammaExpert(1.5*μ, 2.0*σ);
        #          LRMoE.ZIGammaExpert(0.4, μ, σ) LRMoE.GammaExpert(0.5*μ, 0.6*σ) LRMoE.ZIGammaExpert(0.80, 1.5*μ, 2.0*σ)]

        # pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
        #               [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0
        # model_guess = [LRMoE.GammaExpert(0.8*μ, 1.2*σ) LRMoE.GammaExpert(μ, 0.9*σ) LRMoE.GammaExpert(1.0*μ, 2.5*σ);
        #                LRMoE.ZIGammaExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.GammaExpert(0.75*μ, 0.3*σ) LRMoE.ZIGammaExpert(0.50, 1.75*μ, 1.0*σ)]

    #     # Exact observation
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     # With censoring
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), 0.75.*Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

    #     # With truncation
    #     Y = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     Y = hcat( 0.80.*Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], 2.50.*Y_sim[:,2])
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], 2.0 .*Y_sim[:,1], fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ecm_iter_max = 20)

    #     gate_tmp = LogitGating(α_guess, X)
    #     ll_tmp = loglik_np(Y, gate_tmp, model_guess)
    #     hcat(Y, EM_E_k(ll_tmp.gate_expert_tn_bar_k))

    # end

end

@testset "fitting simulated: inversegaussian" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

    #     μ = 2
    #     σ = 3

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.InverseGaussianExpert(μ, σ) LRMoE.InverseGaussianExpert(0.5*μ, 0.6*σ) LRMoE.InverseGaussianExpert(1.5*μ, 2.0*σ);
    #              LRMoE.ZIInverseGaussianExpert(0.4, μ, σ) LRMoE.InverseGaussianExpert(0.5*μ, 0.6*σ) LRMoE.ZIInverseGaussianExpert(0.80, 1.5*μ, 2.0*σ)]

    #     pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
    #                   [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.InverseGaussianExpert(0.8*μ, 1.2*σ) LRMoE.InverseGaussianExpert(μ, 0.9*σ) LRMoE.InverseGaussianExpert(1.0*μ, 2.5*σ);
    #                    LRMoE.ZIInverseGaussianExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.InverseGaussianExpert(0.75*μ, 0.3*σ) LRMoE.ZIInverseGaussianExpert(0.50, 1.75*μ, 1.0*σ)]

    #     # Exact observation
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     # With censoring
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), 0.75.*Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

    #     # With truncation
    #     Y = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     Y = hcat( 0.80.*Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], 2.50.*Y_sim[:,2])
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], 2.0 .*Y_sim[:,1], fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ecm_iter_max = 20)

    #     gate_tmp = LogitGating(α_guess, X)
    #     ll_tmp = loglik_np(Y, gate_tmp, model_guess)
    #     hcat(Y, EM_E_k(ll_tmp.gate_expert_tn_bar_k))

    # end

end

@testset "fitting simulated: weibull" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

    #     μ = 2
    #     σ = 3

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.WeibullExpert(μ, σ) LRMoE.WeibullExpert(0.5*μ, 0.6*σ) LRMoE.WeibullExpert(1.5*μ, 2.0*σ);
    #              LRMoE.ZIWeibullExpert(0.4, μ, σ) LRMoE.WeibullExpert(0.5*μ, 0.6*σ) LRMoE.ZIWeibullExpert(0.80, 1.5*μ, 2.0*σ)]

    #     pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
    #                   [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.WeibullExpert(0.8*μ, 1.2*σ) LRMoE.WeibullExpert(μ, 0.9*σ) LRMoE.WeibullExpert(1.0*μ, 2.5*σ);
    #                    LRMoE.ZIWeibullExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.WeibullExpert(0.75*μ, 0.3*σ) LRMoE.ZIWeibullExpert(0.50, 1.75*μ, 1.0*σ)]

    #     # Exact observation
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     # With censoring
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), 0.75.*Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(Inf, length(Y_sim[:,1])), fill(0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

    #     # With truncation
    #     Y = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     Y = hcat( 0.80.*Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.30.*Y_sim[:,2], 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.80.*Y_sim[:,2], 1.20.*Y_sim[:,2], 2.50.*Y_sim[:,2])
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ϵ = 0.01)

    #     Y = hcat(fill(0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], 2.0 .*Y_sim[:,1], fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params, ecm_iter_max = 20)

    #     gate_tmp = LogitGating(α_guess, X)
    #     ll_tmp = loglik_np(Y, gate_tmp, model_guess)
    #     hcat(Y, EM_E_k(ll_tmp.gate_expert_tn_bar_k))

    # end

end

@testset "fitting simulated: lognormal, exact" begin

    d = Distributions.Gamma(1.0, 2.0)

    μμ = rand(d, 1)
    σσ = rand(d, 1)

    for (μ, σ) in zip(μμ, σσ)

        X = rand(Uniform(-1, 1), 20000, 5)
        α_true = rand(Uniform(-1, 1), 3, 5)
        α_true[3, :] .= 0.0
        model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ);
                 LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.80, 1.5*μ, 2.0*σ)]

        pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf], [Inf 1.0 Inf]],
                      [[Inf 1.0 Inf], [Inf 1.0 Inf], [Inf 1.0 Inf]]]

        Y_sim = sim_dataset(α_true, X, model)
        α_guess = copy(α_true)
        α_guess .= 0.0
        model_guess = [LRMoE.LogNormalExpert(0.8*μ, 1.2*σ) LRMoE.LogNormalExpert(μ, 0.9*σ) LRMoE.LogNormalExpert(1.0*μ, 2.5*σ);
                       LRMoE.ZILogNormalExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.LogNormalExpert(0.75*μ, 0.3*σ) LRMoE.ZILogNormalExpert(0.50, 1.75*μ, 1.0*σ)]

        # Exact observation
        Y = hcat(Y_sim[:,1], Y_sim[:,2])
        result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

    end

end

@testset "fitting simulated: poisson, exact" begin

    # d = Distributions.Poisson(10)

    # λλ = rand(d, 1)

    # for λ in λλ

        X = rand(Uniform(-1, 1), 20000, 5)
        α_true = rand(Uniform(-1, 1), 3, 5)
        α_true[3, :] .= 0.0
        model = [LRMoE.PoissonExpert(λ) LRMoE.PoissonExpert(0.5*λ) LRMoE.PoissonExpert(1.5*λ);
                 LRMoE.ZIPoissonExpert(0.4, λ) LRMoE.PoissonExpert(1.2*λ) LRMoE.ZIPoissonExpert(0.80, 2.0*λ)]

        pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
                      [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

        pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
                      [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

        Y_sim = sim_dataset(α_true, X, model)
        α_guess = copy(α_true)
        α_guess .= 0.0
        model_guess = [LRMoE.PoissonExpert(0.8*λ) LRMoE.PoissonExpert(λ) LRMoE.PoissonExpert(1.0*λ);
                       LRMoE.ZIPoissonExpert(0.50, 2.0*λ) LRMoE.PoissonExpert(0.75*λ) LRMoE.ZIPoissonExpert(0.50, 1.75*λ)]

        # Exact observation
        Y = hcat(Y_sim[:,1], Y_sim[:,2])
        result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.ZIPoissonExpert(0.4, λ) LRMoE.PoissonExpert(1.2*λ) LRMoE.ZIPoissonExpert(0.80, 2.0*λ)]

    #     pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.ZIPoissonExpert(0.50, 2.0*λ) LRMoE.PoissonExpert(0.75*λ) LRMoE.ZIPoissonExpert(0.50, 1.75*λ)]

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1], Y_sim[:,2])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    # end

end

@testset "fitting simulated: binomial, exact" begin

    # d = Distributions.Poisson(10)

    # λλ = rand(d, 1)

    # for λ in λλ

    #     n = 30
    #     p = 0.2

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
        # model = [LRMoE.BinomialExpert(n, p) LRMoE.BinomialExpert(n, 0.5*p) LRMoE.BinomialExpert(n, 1.3*p);
        #          LRMoE.ZIBinomialExpert(0.4, n, p) LRMoE.BinomialExpert(n, 1.2*p) LRMoE.ZIBinomialExpert(0.70, n, 2.0*p)]

    #     pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
    #                   [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

    #     pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
    #                   [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0

    #     model_guess = [LRMoE.BinomialExpert(n, 0.8*p) LRMoE.BinomialExpert(n, p) LRMoE.BinomialExpert(n, 1.0*p);
    #              LRMoE.ZIBinomialExpert(0.5, n, 2.0*p) LRMoE.BinomialExpert(n, 0.75*p) LRMoE.ZIBinomialExpert(0.50, n, 1.75*p)]

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1], Y_sim[:,2])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # n = 30
        # p = 0.2

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 2, 5)
        # α_true[2, :] .= 0.0
        # model = [LRMoE.BinomialExpert(n, p) LRMoE.BinomialExpert(n, 0.0);
        #          LRMoE.ZIBinomialExpert(0.4, n, p) LRMoE.ZIBinomialExpert(0.80, n, 2.0*p)]

        # pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
        #               [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

        # pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
        #               [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0

        # model_guess = [LRMoE.BinomialExpert(n, 0.8*p) LRMoE.BinomialExpert(n, p);
        #          LRMoE.ZIBinomialExpert(0.5, n, 2.0*p) LRMoE.ZIBinomialExpert(0.50, n, 1.75*p)]

        # # Exact observation
        # Y = hcat(Y_sim[:,1], Y_sim[:,2])
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     n = 30
    #     p = 0.2

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 2, 5)
    #     α_true[2, :] .= 0.0
    #     # model = reshape([LRMoE.BinomialExpert(n, p) LRMoE.ZIBinomialExpert(0.30, n, 0.80*p)], 1, 2)
    #     model = reshape([LRMoE.BinomialExpert(n, p) LRMoE.ZIBinomialExpert(0.30, n, 0.80*p)], 1, 2)

    #     pen_params = [[[2.0 1.0], [2.0 1.0]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0

    #     # model = reshape([LRMoE.BinomialExpert(n, 2*p) LRMoE.ZIBinomialExpert(0.50, n, 0.80*p)], 1, 2)
    #     model_guess = reshape([LRMoE.BinomialExpert(n, 2*p) LRMoE.ZIBinomialExpert(0.10, n, 0.80*p)], 1, 2)

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     gate_tmp = LogitGating(α_guess, X)
    #     ll_tmp = loglik_exact(Y, gate_tmp, model_guess)

    # end

end

@testset "fitting simulated: negativebinomial, exact" begin

    # d = Distributions.Poisson(10)

    # λλ = rand(d, 1)

    # for λ in λλ

        # n = 30
        # p = 0.2

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 3, 5)
        # α_true[3, :] .= 0.0
        # model = [LRMoE.NegativeBinomialExpert(n, p) LRMoE.NegativeBinomialExpert(n, 0.2*p) LRMoE.NegativeBinomialExpert(n+5, 1.3*p);
        #          LRMoE.ZINegativeBinomialExpert(0.4, n, p) LRMoE.NegativeBinomialExpert(n+2, 1.2*p) LRMoE.ZINegativeBinomialExpert(0.70, n+5, 2.0*p)]

        # pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
        #               [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

        # pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
        #               [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0

        # model_guess = [LRMoE.NegativeBinomialExpert(n, 0.8*p) LRMoE.NegativeBinomialExpert(n, p) LRMoE.NegativeBinomialExpert(n+2, 1.0*p);
        #          LRMoE.ZINegativeBinomialExpert(0.5, n, 2.0*p) LRMoE.NegativeBinomialExpert(n+1, 0.75*p) LRMoE.ZINegativeBinomialExpert(0.50, n+8, 1.75*p)]

        # model_guess = [LRMoE.NegativeBinomialExpert(n, 0.8*p) LRMoE.NegativeBinomialExpert(n, p) LRMoE.NegativeBinomialExpert(n+2, 1.0*p);
        #          LRMoE.ZINegativeBinomialExpert(0.5, n, 2.0*p) LRMoE.ZINegativeBinomialExpert(0.50, n+1, 0.75*p) LRMoE.ZINegativeBinomialExpert(0.50, n+8, 1.75*p)]

        # # Exact observation
        # Y = hcat(Y_sim[:,1], Y_sim[:,2])
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # n = 10 # 30
        # p = 0.2

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 2, 5)
        # α_true[2, :] .= 0.0
        # model = [LRMoE.NegativeBinomialExpert(n, p) LRMoE.NegativeBinomialExpert(n+5, 0.80*p);
        #          LRMoE.ZINegativeBinomialExpert(0.4, n, p) LRMoE.ZINegativeBinomialExpert(0.20, n+6, 2.0*p)]

        # pen_params = [[[2.0 1.0], [2.0 1.0]],
        #               [[2.0 1.0], [2.0 1.0]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0

        # model_guess = [LRMoE.NegativeBinomialExpert(n, 0.5*p) LRMoE.NegativeBinomialExpert(n, p);
        #          LRMoE.ZINegativeBinomialExpert(0.5, n, 2.0*p) LRMoE.ZINegativeBinomialExpert(0.50, n, 1.75*p)]

        # # Exact observation
        # Y = hcat(Y_sim[:,1], Y_sim[:,2])
        # Y = hcat(fill(0, length(Y_sim[:,1])), floor.(0.60 .* Y_sim[:,1]), ceil.(1.50.*Y_sim[:,1]), fill(Inf, length(Y_sim[:,1])))
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # n = 10 # 30
        # p = 0.2

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 2, 5)
        # α_true[2, :] .= 0.0
        # model = [LRMoE.NegativeBinomialExpert(2*n, p) LRMoE.NegativeBinomialExpert(n, 0.80*p);
        #          LRMoE.ZINegativeBinomialExpert(0.4, n, p) LRMoE.ZINegativeBinomialExpert(0.60, 3*n, 2.0*p)]

        # pen_params = [[[2.0 1.0], [2.0 1.0], [2.0 1.0]],
        #               [[2.0 1.0], [2.0 1.0], [2.0 1.0]]]

        # pen_params = [[[2.0 10.0], [2.0 10.0], [2.0 10.0]],
        #               [[2.0 10.0], [2.0 10.0], [2.0 10.0]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0

        # model_guess = [LRMoE.NegativeBinomialExpert(2*n, 0.5*p) LRMoE.NegativeBinomialExpert(n, p);
        #          LRMoE.ZINegativeBinomialExpert(0.5, n, 0.8*p) LRMoE.ZINegativeBinomialExpert(0.50, 3*n, 1.75*p)]

        # # Exact observation
        # Y = hcat(Y_sim[:,1], Y_sim[:,2])
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), floor.(0.60 .* Y_sim[:,1]), ceil.(1.50.*Y_sim[:,1]), fill(Inf, length(Y_sim[:,1])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # Y = hcat( floor.(0.20 .* Y_sim[:,1]), floor.(0.60 .* Y_sim[:,1]), ceil.(1.50.*Y_sim[:,1]), ceil.(5.0 .* Y_sim[:,1]))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    # end

end

@testset "fitting simulated: gammacount, exact" begin

    # d = Distributions.Poisson(10)

    # λλ = rand(d, 1)

    # for λ in λλ

        # n = 10 # 30
        # p = 3

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 2, 5)
        # α_true[2, :] .= 0.0
        # model = [LRMoE.GammaCountExpert(n, p) LRMoE.GammaCountExpert(n, 0.80*p);
        #          LRMoE.ZIGammaCountExpert(0.4, n, p) LRMoE.ZIGammaCountExpert(0.20, n, 2.0*p)]

        # pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
        #               [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0

        # model_guess = [LRMoE.GammaCountExpert(n, 0.5*p) LRMoE.GammaCountExpert(n, p);
        #          LRMoE.ZIGammaCountExpert(0.5, n, 2.0*p) LRMoE.ZIGammaCountExpert(0.50, n, 1.75*p)]

        # # Exact observation
        # Y = hcat(Y_sim[:,1], Y_sim[:,2])
        # Y = hcat(fill(0, length(Y_sim[:,1])), floor.(0.60 .* Y_sim[:,1]), ceil.(1.50.*Y_sim[:,1]), fill(Inf, length(Y_sim[:,1])))
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # n = 10 # 30
        # p = 3

        # X = rand(Uniform(-1, 1), 20000, 5)
        # α_true = rand(Uniform(-1, 1), 2, 5)
        # α_true[2, :] .= 0.0
        # model = [LRMoE.GammaCountExpert(2*n, p) LRMoE.GammaCountExpert(n, 0.80*p);
        #          LRMoE.ZIGammaCountExpert(0.4, n, p) LRMoE.ZIGammaCountExpert(0.60, 3*n, 2.0*p)]

        # pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
        #               [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

        # Y_sim = sim_dataset(α_true, X, model)
        # α_guess = copy(α_true)
        # α_guess .= 0.0

        # model_guess = [LRMoE.GammaCountExpert(2*n, 0.5*p) LRMoE.GammaCountExpert(n, p);
        #          LRMoE.ZIGammaCountExpert(0.5, n, 0.8*p) LRMoE.ZIGammaCountExpert(0.50, 3*n, 1.75*p)]

        # gate_tmp = LogitGating(α_guess, X)
        # ll_tmp = loglik_exact(Y, gate_tmp, model_guess)

        # ll_tmp = loglik_np(Y, gate_tmp, model_guess)

        # # Exact observation
        # Y = hcat(Y_sim[:,1], Y_sim[:,2])
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        # result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

        # Y = hcat(fill(0, length(Y_sim[:,1])), floor.(0.60 .* Y_sim[:,1]), ceil.(1.50.*Y_sim[:,1]), fill(Inf, length(Y_sim[:,1])))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)

        # Y = hcat( floor.(0.20 .* Y_sim[:,1]), floor.(0.60 .* Y_sim[:,1]), ceil.(1.50.*Y_sim[:,1]), ceil.(5.0 .* Y_sim[:,1]))
        # result = fit_main(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    # end

end

@testset "fitting simulated: gamma, exact" begin

    d = Distributions.Gamma(1.0, 2.0)

    μμ = rand(d, 1)
    σσ = rand(d, 1)

    for (μ, σ) in zip(μμ, σσ)

        X = rand(Uniform(-1, 1), 20000, 5)
        α_true = rand(Uniform(-1, 1), 3, 5)
        α_true[3, :] .= 0.0
        model = [LRMoE.GammaExpert(μ, σ) LRMoE.GammaExpert(0.5*μ, 0.6*σ) LRMoE.GammaExpert(1.5*μ, 2.0*σ);
                 LRMoE.ZIGammaExpert(0.4, μ, σ) LRMoE.GammaExpert(0.5*μ, 0.6*σ) LRMoE.ZIGammaExpert(0.80, 1.5*μ, 2.0*σ)]

        pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
                      [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

        Y_sim = sim_dataset(α_true, X, model)
        α_guess = copy(α_true)
        α_guess .= 0.0
        model_guess = [LRMoE.GammaExpert(0.8*μ, 1.2*σ) LRMoE.GammaExpert(μ, 0.9*σ) LRMoE.GammaExpert(1.0*μ, 2.5*σ);
                       LRMoE.ZIGammaExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.GammaExpert(0.75*μ, 0.3*σ) LRMoE.ZIGammaExpert(0.50, 1.75*μ, 1.0*σ)]

        # Exact observation
        Y = hcat(Y_sim[:,1], Y_sim[:,2])
        result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
        result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    end

end

@testset "fitting simulated: inversegaussian, exact" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.InverseGaussianExpert(μ, σ) LRMoE.InverseGaussianExpert(0.5*μ, 0.6*σ) LRMoE.InverseGaussianExpert(1.5*μ, 2.0*σ);
    #              LRMoE.ZIInverseGaussianExpert(0.4, μ, σ) LRMoE.InverseGaussianExpert(0.5*μ, 0.6*σ) LRMoE.ZIInverseGaussianExpert(0.80, 1.5*μ, 2.0*σ)]

    #     pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
    #                   [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.InverseGaussianExpert(0.8*μ, 1.2*σ) LRMoE.InverseGaussianExpert(μ, 0.9*σ) LRMoE.InverseGaussianExpert(1.0*μ, 2.5*σ);
    #                    LRMoE.ZIInverseGaussianExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.InverseGaussianExpert(0.75*μ, 0.3*σ) LRMoE.ZIInverseGaussianExpert(0.50, 1.75*μ, 1.0*σ)]

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1], Y_sim[:,2])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    # end

end

@testset "fitting simulated: weibull, exact" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

    #     μ = 3
    #     σ = 8
    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.WeibullExpert(μ, σ) LRMoE.WeibullExpert(0.5*μ, 0.6*σ) LRMoE.WeibullExpert(1.5*μ, 2.0*σ);
    #              LRMoE.ZIWeibullExpert(0.4, μ, σ) LRMoE.WeibullExpert(0.5*μ, 0.6*σ) LRMoE.ZIWeibullExpert(0.80, 1.5*μ, 2.0*σ)]

    #     pen_params = [[[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]],
    #                   [[1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.WeibullExpert(0.8*μ, 1.2*σ) LRMoE.WeibullExpert(μ, 0.9*σ) LRMoE.WeibullExpert(1.0*μ, 2.5*σ);
    #                    LRMoE.ZIWeibullExpert(0.50, 2.0*μ, 1.2*σ) LRMoE.WeibullExpert(0.75*μ, 0.3*σ) LRMoE.ZIWeibullExpert(0.50, 1.75*μ, 1.0*σ)]

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1], Y_sim[:,2])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     gate_tmp = LogitGating(α_guess, X)
    #     ll_tmp = loglik_np(Y, gate_tmp, model_guess)
    #     hcat(Y, ll_tmp.expert_ll_comp)

    #     ll_comp = ll_tmp.expert_ll_comp
    #     pos_idx = Y_sim[:,1] .!= 0.0

    # end

end

@testset "fitting simulated: burr, exact" begin

    # d = Distributions.Gamma(1.0, 2.0)

    # μμ = rand(d, 1)
    # σσ = rand(d, 1)

    # for (μ, σ) in zip(μμ, σσ)

    #     k = 3
    #     c = 5
    #     λ = 10
    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 3, 5)
    #     α_true[3, :] .= 0.0
    #     model = [LRMoE.BurrExpert(k, c, λ) LRMoE.BurrExpert(0.5*k, 0.5*c, 2*λ) LRMoE.BurrExpert(1.5*k, 2*c, 3*λ);
    #              LRMoE.ZIBurrExpert(0.4, k, c, λ) LRMoE.BurrExpert(2*k, 0.5*c, 0.6*λ) LRMoE.ZIBurrExpert(0.80, k, 1.5*c, 2.0*λ)]

    #     pen_params = [[[1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf]],
    #                   [[1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf]]]

    #     pen_params = [[[2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0]],
    #                   [[2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.BurrExpert(k, c, λ) LRMoE.BurrExpert(0.5*k, 0.5*c, 2*λ) LRMoE.BurrExpert(1.5*k, 2*c, 3*λ);
    #              LRMoE.ZIBurrExpert(0.4, k, 2*c, λ) LRMoE.BurrExpert(k, c, 0.6*λ) LRMoE.ZIBurrExpert(0.50, k, 2*c, 5*λ)]

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1], Y_sim[:,2])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     gate_tmp = LogitGating(α_guess, X)
    #     ll_tmp = loglik_np(Y, gate_tmp, model_guess)
    #     hcat(Y, ll_tmp.expert_ll_comp)

    #     ll_comp = ll_tmp.expert_ll_comp
    #     pos_idx = Y_sim[:,1] .!= 0.0

    #     k = 3
    #     c = 5
    #     λ = 10
    #     X = rand(Uniform(-1, 1), 20000, 5)
    #     α_true = rand(Uniform(-1, 1), 2, 5)
    #     α_true[2, :] .= 0.0
    #     model = [LRMoE.BurrExpert(k, c, λ) LRMoE.BurrExpert(0.5*k, 0.5*c, 2*λ)]

    #     pen_params = [[[1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf]],
    #                   [[1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf], [1.0 Inf 1.0 Inf 1.0 Inf]]]

    #     pen_params = [[[2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0]],
    #                   [[2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0], [2.0 1.0 2.0 1.0 2.0 1.0]]]

    #     Y_sim = sim_dataset(α_true, X, model)
    #     α_guess = copy(α_true)
    #     α_guess .= 0.0
    #     model_guess = [LRMoE.BurrExpert(k, 2*c, 0.8*λ) LRMoE.BurrExpert(1.5*k, c, 3*λ)]

    #     # Exact observation
    #     Y = hcat(Y_sim[:,1])
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = false, pen_params = pen_params)
    #     result = fit_exact(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)

    #     Y = hcat(fill(0.0, length(Y_sim[:,1])), 0.80 .* Y_sim[:,1], 1.20 .* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])))
    #     result = fit_main(Y, X, α_guess, model_guess, penalty = true, pen_α = 10, pen_params = pen_params)
    # end

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
"""
    matte_example(example, path; force = false, activate = true)

Create a new Matte app, based on the example specified by `example`. See the manual for a
list of examples that you can use here.

`path` should be a folder that doesn't exist. If it doesnt, you can use `force = true` to
overwrite. `activate` controls whether the helper will activate the project for the newly
created app.
"""
function matte_example(example, path; force = false, activate = true)
    src = joinpath(dirname(pathof(Matte)), "..", "examples", example)
    cp(src, path)
    println("App created.\nChanging working directory to `$path`...")
    cd(path)
    if activate
        println("And activating the project...")
        Pkg.activate(".")
    end
    println("Done!")
end

"""
    new_matte_app(path; force = false, activate = true)

Create a new skeleton Matte app in the folder `path`. The folder should be empty. Matte will
create the final folder if it doesn't exist, but will fail if the rest of the path doesn't
exist.

`path` should be a folder that doesn't exist. If it doesnt, you can use `force = true` to
overwrite. `activate` controls whether the helper will activate the project for the newly
created app.
"""
function new_matte_app(path; force = false, activate = true)
    matte_example("blank", path, force = force, activate = activate)
end

ext: jl


text:
"""
    SeisWriteHeaders(filename,h;<keyword arguments>)

Write seismic headers in seis format

# Arguments

- `itrace=1` : First trace number to write
- `update_tracenum=true`

*Credits: AS, 2015*

"""
function SeisWriteHeaders(filename,h;itrace=1,update_tracenum=true)

    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_d = join([DATAPATH filename "@data@"])
    filename_h = join([DATAPATH filename "@headers@"])
    if (itrace==1)
	stream_hout = open(filename_h,"w")
    else
	stream_hout = open(filename_h,"a")
    end
    nx = length(h)
    h1 = Header32Bits[]
    for j = itrace : itrace + nx - 1
	if update_tracenum == true
	    h[j - itrace + 1].tracenum = j
	end
	h2 = HeaderToBits(h[j - itrace + 1])
	append!(h1,h2)
    end
    write(stream_hout,h1)
    close(stream_hout)

end


text:
###### NORMAL SHOCK RELATIONS #######

function MachNumber(M1::MachNumber,::Type{NormalShock};gas::PerfectGas=DefaultPerfectGas)
    value(M1) >= 1.0 || error("M1 must be at least 1")
    γ = SpecificHeatRatio(gas)
    M2sq = (1 + 0.5*(γ-1)*M1^2)/(γ*M1^2-0.5*(γ-1))
    return MachNumber(sqrt(M2sq))
end

function TemperatureRatio(M1::MachNumber,::Type{NormalShock};gas::PerfectGas=DefaultPerfectGas)
    value(M1) >= 1.0 || error("M1 must be at least 1")
    γ = SpecificHeatRatio(gas)
    return TemperatureRatio(1 + 2*(γ-1)/(γ+1)^2*((1+γ*M1^2)/M1^2)*(M1^2-1))
end

function PressureRatio(M1::MachNumber,::Type{NormalShock};gas::PerfectGas=DefaultPerfectGas)
    value(M1) >= 1.0 || error("M1 must be at least 1")
    γ = SpecificHeatRatio(gas)
    return PressureRatio(1 + 2*γ/(γ+1)*(M1^2-1))
end

function DensityRatio(M1::MachNumber,::Type{NormalShock};gas::PerfectGas=DefaultPerfectGas)
    value(M1) >= 1.0 || error("M1 must be at least 1")
    γ = SpecificHeatRatio(gas)
    return DensityRatio((γ+1)*M1^2/(2+(γ-1)*M1^2))
end

function StagnationPressureRatio(M1::MachNumber,::Type{NormalShock};gas::PerfectGas=DefaultPerfectGas)
    value(M1) >= 1.0 || error("M1 must be at least 1")
    M2 = MachNumber(M1,NormalShock,gas=gas)
    pratio = PressureRatio(M1,NormalShock,gas=gas)
    p01_unit = StagnationPressure(Pressure(1),M1,Isentropic,gas=gas)
    p02_unit = StagnationPressure(Pressure(1),M2,Isentropic,gas=gas)
    return StagnationPressureRatio(pratio*p02_unit/p01_unit)
end

function Entropy(s1::Entropy,M1::MachNumber,::Type{NormalShock};gas::PerfectGas=DefaultPerfectGas)
    value(M1) >= 1.0 || error("M1 must be at least 1")
    γ = SpecificHeatRatio(gas)
    cv = SpecificHeatVolume(gas)
    ds = cv*(log(1+2*γ/(γ+1)*(M1^2-1)) - γ*log((γ+1)*M1^2/(2+(γ-1)*M1^2)))
    return Entropy(s1+ds)

end


text:
module DiffEqFlux

using GalacticOptim, DataInterpolations, DiffEqBase, DiffResults, DiffEqSensitivity,
      Distributions, ForwardDiff, Flux, Requires, Adapt, LinearAlgebra, RecursiveArrayTools,
      StaticArrays, Base.Iterators, Printf, Zygote

using DistributionsAD
import ProgressLogging, ZygoteRules
import ConsoleProgressMonitor, TerminalLoggers, LoggingExtras

import Logging

gpu_or_cpu(x) = Array

function diffeq_fd(p,f,n,prob,solver=nothing,args...;u0=prob.u0,kwargs...)
  @warn("diffeq_fd has been deprecated in the update of DiffEqFlux to Zygote support. Use the solve function with sensealg=ForwardDiffSensitivity() to recover the same functionality. See https://docs.juliadiffeq.org/latest/analysis/sensitivity/ for more details")
  f(concrete_solve(prob,solver,u0,p,args...;kwargs...))
end

function diffeq_rd(p,prob,solver=nothing,args...;u0=prob.u0,kwargs...)
  @warn("diffeq_rd has been deprecated in the update of DiffEqFlux to Zygote support. Use the concrete_solve function with sensealg=TrackerAdjoint() to recover the same functionality. See https://docs.juliadiffeq.org/latest/analysis/sensitivity/ for more details")
  concrete_solve(prob,solver,u0,p,args...;kwargs...)
end

function diffeq_adjoint(p,prob,solver=nothing,args...;u0=prob.u0,kwargs...)
  @warn("diffeq_adjoint has been deprecated in the update of DiffEqFlux to Zygote support. Use the concrete_solve function to recover the same functionality. See https://docs.juliadiffeq.org/latest/analysis/sensitivity/ for more details")
  concrete_solve(prob,solver,u0,p,args...;kwargs...)
end

function neural_ode(args...;kwargs...)
  @error("neural_ode has be removed and replaced by NeuralODE. Please consult the README for more details.")
end

function neural_ode_rd(args...;kwargs...)
  @error("neural_ode_rd has be removed and replaced by NeuralODE. Please consult the README for more details.")
end

function neural_dmsde(args...;kwargs...)
  @error("neural_dmsde has be removed and replaced by NeuralDMSDE. Please consult the README for more details.")
end

function neural_dmsde(model,x,mp,tspan,
                      args...;kwargs...)
    error("neural_dmsde has been deprecated with the change to Zygote. Please see the documentation on the new NeuralDSDE layer.")
end

function neural_ode_rd(model,x,tspan,
                       args...;
                       kwargs...)
    error("neural_ode_rd has been deprecated with the change to Zygote. Please see the documentation on the new NeuralODE layer.")
end

function neural_ode(model,x,tspan,args...;kwargs...)
    error("neural_ode has been deprecated with the change to Zygote. Please see the documentation on the new NeuralODE layer.")
end

# ForwardDiff integration

ZygoteRules.@adjoint function ForwardDiff.Dual{T}(x, ẋ::Tuple) where T
  @assert length(ẋ) == 1
  ForwardDiff.Dual{T}(x, ẋ), ḋ -> (ḋ.partials[1], (ḋ.value,))
end

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(d::ForwardDiff.Dual{T}, ::Val{:partials}) where T =
  d.partials, ṗ -> (ForwardDiff.Dual{T}(ṗ[1], 0),)

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(d::ForwardDiff.Dual{T}, ::Val{:value}) where T =
  d.value, ẋ -> (ForwardDiff.Dual{T}(0, ẋ),)

ZygoteRules.@adjoint ZygoteRules.literal_getproperty(A::Tridiagonal, ::Val{:dl}) = A.dl,y -> Tridiagonal(dl,zeros(length(d)),zeros(length(du)),)
ZygoteRules.@adjoint ZygoteRules.literal_getproperty(A::Tridiagonal, ::Val{:d}) = A.d,y -> Tridiagonal(zeros(length(dl)),d,zeros(length(du)),)
ZygoteRules.@adjoint ZygoteRules.literal_getproperty(A::Tridiagonal, ::Val{:du}) = A.dl,y -> Tridiagonal(zeros(length(dl)),zeros(length(d),du),)
ZygoteRules.@adjoint Tridiagonal(dl, d, du) = Tridiagonal(dl, d, du), p̄ -> (diag(p̄[2:end,1:end-1]),diag(p̄),diag(p̄[1:end-1,2:end]))

include("ffjord.jl")
include("train.jl")
include("fast_layers.jl")
include("neural_de.jl")
include("require.jl")
include("spline_layer.jl")
include("tensor_product_basis.jl")
include("tensor_product_layer.jl")
include("collocation.jl")
include("hnn.jl")
include("multiple_shooting.jl")

export diffeq_fd, diffeq_rd, diffeq_adjoint
export DeterministicCNF, FFJORD, NeuralODE, NeuralDSDE, NeuralSDE, NeuralCDDE, NeuralDAE, NeuralODEMM, TensorLayer, AugmentedNDELayer, SplineLayer, NeuralHamiltonianDE
export HamiltonianNN
export ChebyshevBasis, SinBasis, CosBasis, FourierBasis, LegendreBasis, PolynomialBasis
export neural_ode, neural_ode_rd
export neural_dmsde
export FastDense, StaticDense, FastChain, initial_params

export EpanechnikovKernel, UniformKernel, TriangularKernel, QuarticKernel
export TriweightKernel, TricubeKernel, GaussianKernel, CosineKernel
export LogisticKernel, SigmoidKernel, SilvermanKernel
export collocate_data

export multiple_shoot

end


text:
import Remark
using Literate
ENV["GKSwstype"]="100"
using Plots

files =  filter( f -> startswith(f, "0"), readdir("src")) |> collect

for file in files
    Literate.notebook("src/$file", "notebooks",  execute=true)
    slides_path = joinpath("docs",file[1:2])
    mkpath(slides_path)
    s = Remark.slideshow("src/$file", slides_path)
end


text:
using BlackBoxOptim

if length(ARGS) >= 1
  maxtime = parse(Float64, ARGS[1])
else
  maxtime = 0.3 # Time used in orignal (flaky) tests...
end

if length(ARGS) >= 2
  N = parse(Int, ARGS[2])
else
  N = 30
end

schaffer1(x) = (sumabs2(x), sumabs2(x .- 2.0))

results1 = Float64[]
results2 = Float64[]

totaltime = 0.0

info("MaxTime=$maxtime, N_repeats=$N")
while length(results1) < N
  print("Optimization run $(length(results1)+1) of $N: ")
  ctrl = bbsetup(schaffer1; Method = :borg_moea,
    SearchRange = [(-10.0, 10.0), (-10.0, 10.0)], TraceMode = :silent,
    FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true), ϵ=0.01)
  tic()
  @time res = bboptimize(ctrl, MaxTime = maxtime)
  t = toq()
  if isapprox(t, maxtime; rtol=0.01, atol=0.01)
    totaltime += t
    f = best_fitness(res)
    push!(results1, f[1])
    push!(results2, f[2])
  else
    # Ignore runs with large MaxTime deviation:
    # usually initial runs have compilation hiccups (both BBO and GC)
    @warn("Run ignored: large MaxTime deviation")
  end
end

function report(vs, label, t = nothing)
  r(x) = round(x, digits=3)
  println(label, ": ", r(mean(vs)), " ± ", r(std(vs)),
    " (", r(minimum(vs)), "-", r(maximum(vs)), ")")
  if t != nothing
    println("Time: ", r(t))
  end
end

report(results1, "Fitness 1")
report(results2, "Fitness 2", totaltime)

# Results below from running on Robert's MacBook Pro 13" from 2015, on 2016-07-29:

# julia04 -L src/BlackBoxOptim.jl examples/stability_test.jl 0.3 50
# Fitness 1: 1.999 +/- 0.007 (1.974-2.016)
# Fitness 2: 2.001 +/- 0.007 (1.985-2.026)
# Time: 15.028

# julia05 -L src/BlackBoxOptim.jl examples/stability_test.jl 0.3 50
# Fitness 1: 2.031 +/- 0.234 (1.655-2.85)
# Fitness 2: 2.059 +/- 0.423 (1.48-4.653)
# Time: 15.112

# julia04 -L src/BlackBoxOptim.jl examples/stability_test.jl 3.0 50
# Fitness 1: 2.0 +/- 0.005 (1.995-2.005)
# Fitness 2: 2.0 +/- 0.005 (1.995-2.005)
# Time: 147.778

# julia05 -L src/BlackBoxOptim.jl examples/stability_test.jl 3.0 50
# Fitness 1: 2.004 +/- 0.019 (1.955-2.048)
# Fitness 2: 1.996 +/- 0.019 (1.952-2.046)
# Time: 151.186


text:
export kernelN2D, kernelN3D # haven't figured out how best to dispath for kernelN
export transmission_angle_wiener, transmission_angle, transmission_direction

function kernelN2D(n::Int,x::Union{T,Complex{T}},y::Union{T,Complex{T}}) where T<:AbstractFloat
        h = hankelh1(n,x); dh = diffhankelh1(n,x)
        j = besselj(n,y);  dj = diffbesselj(n,y)

    return x * dh * j - y * h * dj
end

function kernelN3D(n::Int,x::Union{T,Complex{T}},y::Union{T,Complex{T}}) where T<:AbstractFloat
    h = shankelh1(n,x); dh = diffshankelh1(n,x)
    j = sbesselj(n,y);  dj = diffsbesselj(n,y)

    return x * dh * j - y * h * dj
end

"""
    transmission_direction(k_eff::Complex, incident_wavevector::AbstractArray, surface_normal::AbstractArray)

Gives the effective direction `direction` where `sum(direction .^ 2) = 1.0` and the components of `k_eff .* direction` and `incident_wavevector` which are orthogonal to the surface are the same. Note `surface_normal` is the outward pointing normal and the incident medium's wavenumber `k = sqrt(sum(incident_wavevector .^2))`. Below we deduce the result.

Assume that `dot(v,w) = conj(v[i])w[i]` and `surface_normal[i]*surface_normal[i] = 1`. Let `wnp` be orthogonal to `surface_normal`

    wnp =  incident_wavevector -  dot(surface_normal,incident_wavevector) .* surface_normal

Then we know that the effective transmission wavevector has to equal

    k_eff_vec = wnp + α .* surface_normal

where α is determined so that

    sum(x^2 for x in wnp) + α^2  = sum(x^2 for x in k_eff_vec) = k_eff^2
"""
function transmission_direction(k_eff::Complex{T}, incident_wavevector::AbstractArray{CT} where CT <: Union{T,Complex{T}}, surface_normal::AbstractArray{T};
        tol::T = sqrt(eps(T))) where {T<:AbstractFloat,Dim}

    # Let us first make the surface normal a unit vector which is outward pointing
    surface_normal = surface_normal / norm(surface_normal)
    if real(dot(surface_normal,incident_wavevector)) > 0
        surface_normal = -surface_normal
    end

    wnp = incident_wavevector -  dot(surface_normal,incident_wavevector) .* surface_normal
    α = sqrt(k_eff^2 - sum(x^2 for x in wnp))

    # The conditions below guarantee (in order of priority) that either: 1) the wave attenuates when travelling into the material (imag(α) < 0), or 2) the wave does not attenuate, but travels into the material.
    if tol * real(α) > abs(imag(α))
        α = -α # leads to real(α) < 0
    elseif imag(α) > 0
        α = -α # leads to imag(α) < 0
    elseif isnan(α)
        α = -one(T) # covers the cases k_eff = Inf and k_eff = 0.0
        return wnp + α .* surface_normal
    end

    return (wnp + α .* surface_normal) ./ k_eff
end

function transmission_direction(k_eff::Complex{T},  psource::PlaneSource{T,Dim}, material::Material{Dim}; tol::T = sqrt(eps(T))) where {T,Dim}
    transmission_direction(k_eff, psource.direction, material.shape.normal; tol = tol)
end

transmission_angle(pwave::Union{PlaneSource,EffectivePlaneWaveMode}, material::Material) = transmission_angle(pwave, material.shape)

transmission_angle(pwave::Union{PlaneSource,EffectivePlaneWaveMode}, shape::Union{Plate,Halfspace}) = transmission_angle(pwave.direction, shape.normal)

"""
    transmission_angle(wavevector, surface_normal)

Calculates the transmission angle for a wave propagting inside a material which has the outward normal `surface_normal`. The wave propagates in the direction of `wavevector` with a wavenumber given by `k = sqrt(sum(wavevector .^2))`, where `k` should have a positive imaginary part.
"""
transmission_angle(wavevector::AbstractVector,surface_normal::AbstractVector) = transmission_angle(SVector(wavevector...),SVector(surface_normal...))

"""
    transmission_angle(wavevector::StaticVector{3}, surface_normal::StaticVector{3})

Some care is needed when `wavevector` is a complex vector in 3 dimensions.

## The maths

Let's define the normal projection `vn`:

    n = - surface_normal
    vn = dot(n,wavevector) .* n

remembering that in Julia, and mathematics, `dot(v,w) = sum(conj(v[i])*w[i] for i in eachindex(v))`.

We now define the orthogonal component `vo` and its direction `no`:

    vo = wavevector - vn
    no = vo / sqrt(sum(vo.^2))

Note this assumes that `vo` is a real vector, which it should be due to symmetry.

To determine the transmission angle `θT` we calculate the wavenumber

    k = ± sqrt(sum(wavevector .^2))

where we choose the sign in `±` such that `imag(k) >= 0`.

From the above we can write the `wavevector` in the form

    wavevector = k .* (n .* cos(θ) + no .* sin(θ))

which implies that

    k * cos(θ) = dot(conj(n), wavevector)
    k * sin(θ) = dot(conj(no),wavevector) = dot(conj(no),vo) = sqrt(sum(vo.^2))

From these two we can robustly calculate `θT` as

    θT = atan(sqrt(sum(vo.^2)),dot(n,wavevector))

"""
function transmission_angle(wavevector::SVector{3,CT} where CT <: Union{T,Complex{T}}, surface_normal::SVector{3,T}) where {T<:AbstractFloat}
    n = - surface_normal / norm(surface_normal)
    if real(dot(n,wavevector)) < 0
        n = -n
    end

    kcosθ = dot(n,wavevector)

    vo = wavevector - kcosθ .* n
    ksinθ = sqrt(sum(vo.^2))

    return θ = atan(ksinθ,kcosθ)
end

"""
    transmission_angle(wavevector::StaticVector{2}, surface_normal::StaticVector{2})

Specialised to 2 spatial dimensions.
"""
function transmission_angle(wavevector::SVector{2,CT} where CT <: Union{T,Complex{T}}, surface_normal::SVector{2,T}) where {T<:AbstractFloat}
    n = - surface_normal
    no = [-n[2], n[1]] # guarantee θ grows anti-clockwise

    θ = atan(dot(conj(no),wavevector),dot(n,wavevector))

    return θ
end

"""
calculate effective transmission angle θ_eff. We restrict -pi/2 < Re θ_eff < pi/2 for 2 dimensions.
"""
function transmission_angle_wiener(k::Union{T,Complex{T}}, k_eff::Union{T,Complex{T}}, θin) where T<:Number
    # snell(θ::Array{T}) = abs(k*sin(θin) - k_eff*sin(θ[1] + im*θ[2]))
    # result = optimize(snell, [θin,0.]; x_tol= tol, g_tol= tol^2.0)

    θ_eff = asin(k * sin(θin) / k_eff)

    # if !(abs(real(θ_eff - θin)) <= pi/T(2))
    if sign(real(k_eff * cos(θ_eff))) != sign(real(k_eff))
        θ_eff = pi - θ_eff
    end

    return θ_eff
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
# must be run in Julia 0.6: PerceptualColourMaps is an out of date package with
# some great functionality
import PerceptualColourMaps: cmap
import Colors: RGB
import Colors

open("reds.txt","w") do stream
  for col in hex.(Colors.colormap("Reds"))
    println(stream,col)
  end
end

open("diverging_colors.txt","w") do stream
  for col in hex.(RGB.(cmap("D1")))
    println(stream,col)
  end
end

open("circular_colors.txt","w") do stream
  for col in hex.(RGB.(cmap("C6")))
    println(stream,col)
  end
end


text: #First try: Just going to allow users to edit cell parameters

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
@expectation function banana()
    x1 ~ Normal(0, 4)
    x2 ~ Normal(0, 4)
    Turing.acclogp!(_varinfo, banana_density(x1, x2))
    return banana_f(x1, x2)
end

function banana_f(x1, x2)
    cond = 1 / (1 + exp(50 * (x2 + 5)))
    return cond * (x1 - 2)^3
end

function banana_density(x1, x2)
    return -0.5 * (0.03 * x1^2 + (x2/2 + 0.03 * (x1^2 - 100))^2 )
end


text:
export SeisChannel

@doc (@doc SeisData)
mutable struct SeisChannel <: GphysChannel
  id    ::String
  name  ::String
  loc   ::InstrumentPosition
  fs    ::Float64
  gain  ::Float64
  resp  ::InstrumentResponse
  units ::String
  src   ::String
  misc  ::Dict{String,Any}
  notes ::Array{String,1}
  t     ::Array{Int64,2}
  x     ::FloatArray

  function SeisChannel(
      id    ::String,
      name  ::String,
      loc   ::InstrumentPosition,
      fs    ::Float64,
      gain  ::Float64,
      resp  ::InstrumentResponse,
      units ::String,
      src   ::String,
      misc  ::Dict{String,Any},
      notes ::Array{String,1},
      t     ::Array{Int64,2},
      x     ::FloatArray
      )

      return new(id, name, loc, fs, gain, resp, units, src, misc, notes, t, x)
    end
end

# Are keywords type-stable now?
SeisChannel(;
            id    ::String              = "",
            name  ::String              = "",
            loc   ::InstrumentPosition  = GeoLoc(),
            fs    ::Float64             = default_fs,
            gain  ::Float64             = default_gain,
            resp  ::InstrumentResponse  = PZResp(),
            units ::String              = "",
            src   ::String              = "",
            misc  ::Dict{String,Any}    = Dict{String,Any}(),
            notes ::Array{String,1}     = Array{String,1}(undef, 0),
            t     ::Array{Int64,2}      = Array{Int64,2}(undef, 0, 2),
            x     ::FloatArray          = Array{Float32,1}(undef, 0)
            ) = SeisChannel(id, name, loc, fs, gain, resp, units, src, misc, notes, t, x)

function getindex(S::SeisData, j::Int)
  C = SeisChannel()
  for f in datafields
    setfield!(C, f, getindex(getfield(S,f), j))
  end
  return C
end
setindex!(S::SeisData, C::SeisChannel, j::Int) = (
  [(getfield(S, f))[j] = getfield(C, f) for f in datafields];
  return S)

function isempty(Ch::SeisChannel)
  q::Bool = min(Ch.gain == default_gain, Ch.fs == default_fs)
  if q == true
    for f in (:id, :loc, :misc, :name, :notes, :resp, :src, :t, :units, :x)
      q = min(q, isempty(getfield(Ch, f)))
    end
  end
  return q
end

# ============================================================================
# Conversion and push to SeisData
function SeisData(C::SeisChannel)
  S = SeisData(1)
  for f in datafields
    setindex!(getfield(S, f), getfield(C, f), 1)
  end
  return S
end

function push!(S::SeisData, C::SeisChannel)
  for i in datafields
    push!(getfield(S,i), getfield(C,i))
  end
  S.n += 1
  return nothing
end

# This intentionally undercounts exotic objects in :misc (e.g. a nested Dict)
# because those objects aren't written to disk or created by SeisIO
function sizeof(C::SeisChannel)
  s = 96
  for f in datafields
    v = getfield(C,f)
    s += sizeof(v)
    if f == :notes
      if !isempty(v)
        s += sum([sizeof(j) for j in v])
      end
    elseif f == :misc
      k = collect(keys(v))
      s += sizeof(k) + 64 + sum([sizeof(j) for j in k])
      for p in values(v)
        s += sizeof(p)
        if typeof(p) == Array{String,1}
          s += sum([sizeof(j) for j in p])
        end
      end
    end
  end
  return s
end

function write(io::IO, S::SeisChannel)
  write(io, Int64(sizeof(S.id)))
  write(io, S.id)                                                     # id
  write(io, Int64(sizeof(S.name)))
  write(io, S.name)                                                   # name
  write(io, loctyp2code(S.loc))
  write(io, S.loc)                                                    # loc
  write(io, S.fs)                                                     # fs
  write(io, S.gain)                                                   # gain
  write(io, resptyp2code(S.resp))
  write(io, S.resp)                                                   # resp
  write(io, Int64(sizeof(S.units)))
  write(io, S.units)                                                  # units
  write(io, Int64(sizeof(S.src)))
  write(io, S.src)                                                    # src
  write_misc(io, S.misc)                                              # misc
  write_string_vec(io, S.notes)                                       # notes
  write(io, Int64(size(S.t,1)))
  write(io, S.t)                                                      # t
  write(io, typ2code(eltype(S.x)))
  write(io, Int64(length(S.x)))
  write(io, S.x)                                                      # x
  return nothing
end

read(io::IO, ::Type{SeisChannel}) = SeisChannel(
  String(fastread(io, fastread(io, Int64))),                          # id
  String(fastread(io, fastread(io, Int64))),                          # name
  read(io, code2loctyp(fastread(io))),                                # loc
  fastread(io, Float64),                                              # fs
  fastread(io, Float64),                                              # gain
  read(io, code2resptyp(fastread(io))),                               # resp
  String(fastread(io, fastread(io, Int64))),                          # units
  String(fastread(io, fastread(io, Int64))),                          # src
  read_misc(io, getfield(BUF, :buf)),                                 # misc
  read_string_vec(io, getfield(BUF, :buf)),                           # notes
  read!(io, Array{Int64, 2}(undef, fastread(io, Int64), 2)),          # t
  read!(io, Array{code2typ(read(io,UInt8)), 1}(undef, fastread(io, Int64))), # x
  )

convert(::Type{SeisData}, C::SeisChannel) = SeisData(C)


text:
using HTTP, JSON, Base64

###############################################################################
# Module Variables & Constants
###############################################################################

MINIMUM_TIME = -(16 << 56)
MAXIMUM_TIME = (48 << 56) - 1

###############################################################################
# UUID Encode/Decode
###############################################################################

function encodeUUID(txt::String)
    parsed = join(split(txt, "-"))
    parts = [parsed[i:i+1] for i in 1:2:32]
    return base64encode(IOBuffer(hex2bytes(parsed)))
end

function decodeUUID(txt::String)
    function zeropad(val::String)
        if length(val) == 2
            return val
        elseif length(val) == 1
            return "0" * val
        elseif length(val) == 0
            return "00"
        else
            # throw(Exception("invalid String"))
        end
    end
    parts = string.(base64decode(txt), base=16)
    parts = [zeropad(val) for val in parts]
    parts = join([zeropad(val) for val in parts])

    return "$(parts[1:8])-$(parts[9:12])-$(parts[13:16])-$(parts[17:20])-$(parts[21:end])"
end

###############################################################################
# APIFrontEnd Related
###############################################################################

function stitch(parts::Array{String,1}, label::String)
    pattern = Regex("{\"result\": ")
    results = nothing
    for item in parts
        data = JSON.parse(item)["result"]
        if results === nothing
            results = data[label]
        else
            results = vcat(results, data[label])
        end
    end
    return results
end

function parse_api_results(content, label)
    pattern = Regex("{\"result\": ")
    numtokens = length(collect(eachmatch(Regex("{\"result"), content)))

    if numtokens == 1
        data = JSON.parse(content)
        return data["result"][label]
    elseif numtokens > 1
        parts = string.(split(content, "\n"))
        parts = [item for item in parts if item != ""]
        return stitch(parts, label)
    else

    end
end

function apicall(url::String, body::String)
    headers = [
        "Content-Type" => "application/json",
        "Authorization" => "Bearer $APIKEY"
    ]

    io = Base.BufferStream()
    response = HTTP.request("POST", url, headers, body, response_stream=io; verbose=0)
    close(io)

    # TODO
    # do not check for stat error here, as it could be invalid
    # JSON in the body.  We CAN check for HTTP status error
    # though

    return String(read(io))
end

function checkstat(data::Dict)
    if data["stat"] != nothing
        if haskey(data["stat"], "code")
            throw(BTrDBException(data["stat"]["msg"]))
        end
        throw(BTrDBException("unknown error occurred"))
    end
end


text:
@doc raw"""
# Hermite's Interpolating Polynomials

Implements a two point Hermite inter-/extrapolation function which passes
through the function and its first derivative for the interval ``[0,1]``.
The polynomial is determined by four constraint equations, matching the
function and its derivative at the points ``0`` and ``1``.

Call with one of the following methods
```julia
_hermite_extrapolation!(t₀, t₁, x₀, x₁, ẋ₀, ẋ₁, t, x)
_hermite_extrapolation!(t₀, t₁, x₀, x₁, ẋ₀, ẋ₁, t, x, ẋ)
_hermite_extrapolation!(v, t₀, t₁, x₀, x₁, t, x)
_hermite_extrapolation!(v, t₀, t₁, x₀, x₁, t, x, ẋ)
```

where

* `t₀`: first  sample time $t_0$
* `t₁`: second sample time $t_1$
* `x₀`: first  solution value $x_0 = x(t_0)$
* `x₁`: second solution value $x_1 = x(t_1)$
* `ẋ₀`: first  vector field value $ẋ_0 = v(t_0, x(t_0))$
* `ẋ₁`: second vector field value $ẋ_1 = v(t_1, x(t_1))$
* `v`:  function to compute vector field with signature `v(t,x,ẋ)`
* `t`:  time $t$ to extrapolate
* `x`:  extrapolated solution value $x(t)$
* `v`:  extrapolated vector field value $ẋ(t)$

#### Derivation

The interpolation works as follows:
Start by defining the 3rd degree polynomial and its derivative by
```math
\begin{aligned}
g(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 , \\
g'(x) &= a_1 + 2 a_2 x + 3 a_3 x^2 ,
\end{aligned}
```
and apply the constraints
```math
\begin{aligned}
g(0) &= f_0 & & \Rightarrow & a_0 &= f_0 , \\
g(1) &= f_1 & & \Rightarrow & a_0 + a_1 + a_2 + a_3 &= f_1 , \\
g'(0) &= f'_0 & & \Rightarrow & a_1 &= f'_0 , \\
g'(1) &= f'_1 & & \Rightarrow & a_1 + 2 a_2 + 3 a_3 &= f'_1 . \\
\end{aligned}
```
Solving for ``a_0, a_1, a_2, a_3`` leads to
```math
\begin{aligned}
a_0 &= f_0 , &
a_1 &= f'_0 , &
a_2 &= - 3 f_0 + 3 f_1 - 2 f'_0 - f'_1 , &
a_3 &= 2 f_0 - 2 f_1 + f'_0 + f'_1 ,
\end{aligned}
```
so that the polynomial ``g(x)`` reads
```math
g(x) = f_0 + f'_0 x + (- 3 f_0 + 3 f_1 - 2 f'_0 - f'_1) x^2 + (2 f_0 - 2 f_1 + f'_0 + f'_1) x^3 .
```
The function and derivative values can be factored out, so that ``g(x)`` can be rewritten as
```math
g(x) = f_0 (1 - 3 x^2 + 2 x^3) + f_1 (3 x^2 - 2 x^3) + f'_0 (x - 2 x^2 + x^3) + f'_1 (- x^2 + x^3) ,
```
or in generic form as
```math
g(x) = f_0 a_0(x) + f_1 a_1(x) + f'_0 b_0(x) + f'_1 b_1(x) ,
```
with basis functions
```math
\begin{aligned}
a_0 (x) &= 1 - 3 x^2 + 2 x^3 , &
b_0 (x) &= x - 2 x^2 + x^3 , \\
a_1 (x) &= 3 x^2 - 2 x^3 , &
b_1 (x) &= - x^2 + x^3 .
\end{aligned}
```
The derivative ``g'(x)`` accordingly reads
```math
g'(x) = f_0 a'_0(x) + f_1 a'_1(x) + f'_0 b'_0(x) + f'_1 b'_1(x) ,
```
with
```math
\begin{aligned}
a'_0 (x) &= - 6 x + 6 x^2 , &
b'_0 (x) &= 1 - 4 x + 3 x^2 , \\
a'_1 (x) &= 6 x - 6 x^2 , &
b'_1 (x) &= - 2 x + 3 x^2 .
\end{aligned}
```
The basis functions ``a_0``and ``a_1`` are associated with the function
values at ``x_0`` and ``x_1``, respectively, while the basis functions
``b_0`` and ``b_1`` are associated with the derivative values at
``x_0`` and ``x_1``.
The basis functions satisfy the following relations,
```math
\begin{aligned}
a_i (x_j) &= \delta_{ij} , &
b_i (x_j) &= 0 , &
a'_i (x_j) &= 0 , &
b'_i (x_j) &= \delta_{ij} , &
i,j &= 0, 1 ,
\end{aligned}
```
where ``\delta_{ij}`` denotes the Kronecker-delta, so that
```math
\begin{aligned}
g(0) &= f_0 , &
g(1) &= f_1 , &
g'(0) &= f'_0 , &
g'(1) &= f'_1 .
\end{aligned}
```
"""
function _hermite_extrapolation! end

function _hermite_extrapolation!(t₀::TT, t₁::TT, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}, ẋ₀::AbstractArray{DT}, ẋ₁::AbstractArray{DT}, t::TT, x::AbstractArray{DT}) where {DT,TT}
    local a₀::TT
    local a₁::TT
    local b₀::TT
    local b₁::TT
    local Δt = t₁ - t₀
    local s = (t - t₀) / Δt

    # Interpolate x at t
    if t == t₀
        x .= x₀
    elseif t == t₁
        x .= x₁
    else
        a₁ = 3s^2 - 2s^3
        a₀ = 1 - a₁
        b₁ = s^2*(s-1)
        b₀ = s*(1-s)+b₁
        x .= a₀ .* x₀ .+ a₁ .* x₁ .+ b₀ .* Δt .* ẋ₀ .+ b₁ .* Δt .* ẋ₁
    end

    return x
end

function _hermite_extrapolation!(t₀::TT, t₁::TT, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}, ẋ₀::AbstractArray{DT}, ẋ₁::AbstractArray{DT}, t::TT, x::AbstractArray{DT}, ẋ::AbstractArray{DT}) where {DT,TT}
    local a₀::TT
    local a₁::TT
    local b₀::TT
    local b₁::TT
    local Δt = t₁ - t₀
    local s = (t - t₀) / Δt

    _hermite_extrapolation!(t₀, t₁, x₀, x₁, ẋ₀, ẋ₁, t, x)

    # Interpolate ẋ at t
    if t == t₀
        ẋ .= ẋ₀
    elseif t == t₁
        ẋ .= ẋ₁
    else
        a₁ = (6s - 6s^2) / Δt
        a₀ = - a₁
        b₁ = s*(3s-2)
        b₀ = 1-2s+b₁
        ẋ .= a₀ .* x₀ .+ a₁ .* x₁ .+ b₀ .* ẋ₀ .+ b₁ .* ẋ₁
    end

    return (x, ẋ)
end

function _get_velocities(v::Function, t₀::TT, t₁::TT, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}) where {DT,TT}
    ẋ₀ = zero(x₀)
    ẋ₁ = zero(x₁)
    v(t₀, x₀, ẋ₀)
    v(t₁, x₁, ẋ₁)
    return (ẋ₀, ẋ₁)
end

function _hermite_extrapolation!(v::Function, t₀::TT, t₁::TT, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}, t::TT, x::AbstractArray{DT}) where {DT,TT}
    _hermite_extrapolation!(t₀, t₁, x₀, x₁, _get_velocities(v, t₀, t₁, x₀, x₁)..., t, x)
end

function _hermite_extrapolation!(v::Function, t₀::TT, t₁::TT, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}, t::TT, x::AbstractArray{DT}, ẋ::AbstractArray{DT}) where {DT,TT}
    _hermite_extrapolation!(t₀, t₁, x₀, x₁, _get_velocities(v, t₀, t₁, x₀, x₁)..., t, x, ẋ)
end

struct HermiteExtrapolation{T} <: Extrapolation
    t₀::T
    t₁::T
    Δt::T

    function HermiteExtrapolation{T}(t₀, t₁) where {T}
        new(t₀, t₁, t₁-t₀)
    end
end

function HermiteExtrapolation(t₀::T, t₁::T) where {T}
    HermiteExtrapolation{T}(t₀, t₁)
end

function GeometricBase.evaluate!(int::HermiteExtrapolation{TT}, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}, ẋ₀::AbstractArray{DT}, ẋ₁::AbstractArray{DT}, t::TT, x::AbstractArray{DT}) where {DT,TT}
    _hermite_extrapolation!(int.t₀, int.t₁, x₀, x₁, ẋ₀, ẋ₁, t, x)
end

function GeometricBase.evaluate!(int::HermiteExtrapolation{TT}, x₀::AbstractArray{DT}, x₁::AbstractArray{DT}, ẋ₀::AbstractArray{DT}, ẋ₁::AbstractArray{DT}, t::TT, x::AbstractArray{DT}, ẋ::AbstractArray{DT}) where {DT,TT}
    _hermite_extrapolation!(int.t₀, int.t₁, x₀, x₁, ẋ₀, ẋ₁, t, x, ẋ)
end


text:
include("helpers.jl")

function laff_copy(x::T, y::T) where T <: Matrix{Number}

    my, ny = verifydimensions(x, y)

    # Converter x para o formato de y
    x = reshape(x, my, ny)

    # Copiar x para y
    y = copy(x)

    y

end

## Testar ##

## Usaremos matrizes, pois não há propriamente vetores linha e coluna
## em Julia
x = Matrix([1,2,3], 3, 1)
y = Matrix([4,5,6], 1, 3)

y_out = laff_copy(x, y)

max_stars_repo_path: LAFF/Semana 1/laff_copy.jl

max_issues_repo_path: LAFF/Semana 1/laff_copy.jl

max_forks_repo_path: LAFF/Semana 1/laff_copy.jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LearningProblem(sdata, tdata, task)

A geostatistical learning problem with source data `sdata`,
target data `tdata`, and learning `task`.

## References

* Hoffimann et al. 2021. [Geostatistical Learning: Challenges and Opportunities]
  (https://www.frontiersin.org/articles/10.3389/fams.2021.689393/full)
```
"""
struct LearningProblem{Dₛ,Dₜ,T} <: Problem
  sdata::Dₛ
  tdata::Dₜ
  task::T

  function LearningProblem{Dₛ,Dₜ,T}(sdata, tdata, task) where {Dₛ,Dₜ,T}
    sourcevars = name.(variables(sdata))
    targetvars = name.(variables(tdata))

    # assert task is compatible with the data
    @assert features(task) ⊆ sourcevars "features must be present in source data"
    @assert features(task) ⊆ targetvars "features must be present in target data"
    if issupervised(task)
      @assert label(task) ∈ sourcevars "label must be present in source data"
    end

    new(sdata, tdata, task)
  end
end

LearningProblem(sdata::Dₛ, tdata::Dₜ, task::T) where {Dₛ,Dₜ,T} =
  LearningProblem{Dₛ,Dₜ,T}(sdata, tdata, task)

"""
    sourcedata(problem)

Return the source data of the learning `problem`.
"""
sourcedata(problem::LearningProblem) = problem.sdata

"""
    targetdata(problem)

Return the target data of the learning `problem`.
"""
targetdata(problem::LearningProblem) = problem.tdata

"""
    task(problem)

Return the learning task of the learning `problem`.
"""
task(problem::LearningProblem) = problem.task

# ------------
# IO methods
# ------------
function Base.show(io::IO, problem::LearningProblem)
  N = embeddim(problem.sdata)
  print(io, "$(N)D LearningProblem")
end

function Base.show(io::IO, ::MIME"text/plain", problem::LearningProblem)
  println(io, problem)
  println(io, "  source: ", problem.sdata)
  println(io, "  target: ", problem.tdata)
  print(  io, "  task:   ", problem.task)
end


text:
abstract type AbstractSpace{T} end

"""
    EffSpace{T}  <: AbstractSpace{T}

`EffSpace{T}` is a type for efficient design spaces.
This means that they implement the functions  `size`, `length`, `sample`,
`getindex` can be used without the construction of all possible constructs in the
space. Explicit calculation of all constructs is still possible.
"""
abstract type EffSpace{T}  <: AbstractSpace{T} end

"""
  Base.summary(io::IO, a::EffSpace)

provoides a summary of `a`
"""
function Base.summary(io::IO, a::EffSpace)
    space_type,n = typeof(a),length(a)
    print(io,"$n-element $space_type")
end

function Base.show(io::IO,a::EffSpace)
    space_type,construct_type,n = typeof(a),eltype(a),length(a)
    print(io," spacetype | $space_type \n generated constructs | $construct_type \n n_consturcts | $n")
end

"""
    getindex(space::EffSpace,i::Int)

Returns the construct on the i-th possition in the design space.
"""

Base.getindex(space::EffSpace,i::Int) = space.space[i]
"""
    getindex(space::EffSpace,A::Array)

Returns a vector of constructs out of `space` corresponding to the indices in `A`.
"""
Base.getindex(space::EffSpace,A::Array) = map(i -> space.space[i], A)

"""
    length(space::EffSpace)

Returns the number of constructs in the `space`.
"""
Base.length(space::EffSpace) = length(space.space)

"""
    iterate(space::EffSpace, state = 1)

Creates an iterater object for an `EffSpace`.
"""
function Base.iterate(space::EffSpace, state = 1 )
    if  state <= length(space)
        construct = space[state]
        state += 1
        return (construct,state)
    else
        return
    end
end

"""
    FullOrderedSpace{T}

`EffSpace{T}` space type for space filled with `OrderedConstruct` and without any constraints.
"""
struct FullOrderedSpace{T} <: EffSpace{T}
    space::KroneckerPower{T}
end

"""
    eltype(::Type{FullOrderedSpace)

Returns the type of `FullOrderedSpace{T}` if they are collected.
"""
Base.eltype(::Type{FullOrderedSpace{T}}) where {T} = OrderedConstruct{T}

"""
    FullUnorderedSpace{T}

`EffSpace{T}` space type for space filled with `UnorderedConstruct` and without any constraints
"""
struct FullUnorderedSpace{T}<: EffSpace{T}
    space::Combination{T}
end

"""
    eltype(::Type{FullUnorderedSpace{T}}) where {T}

Returns the type of `FullUnorderedSpace{T}` if they are collected.
"""
Base.eltype(::Type{FullUnorderedSpace{T}}) where {T} = UnorderedConstruct{T}

"""
    ComputedSpace{T}

`EffSpace{T}` space type for space where all constructs are explicitly generated.
"""
struct ComputedSpace{T} <: EffSpace{T}
    space::T
end

"""
    eltype(::Type{ComputedSpace)

Return the type of `ComputedSpace` when they are collected.
"""
Base.eltype(::Type{ComputedSpace{T}}) where {T} = eltype(T)

"""
    FrameSpace{T, Tc <: ConstructConstraints}

Structure generated for a space where no efficient alterative is implented currently.
The closed efficient design space is stored together with the uncheck Constraints.
"""
struct FrameSpace{Ts <: EffSpace{T} where T, Tc <: ConstructConstraints} <: AbstractSpace{Ts}
    space::Ts
    con::Tc
end

function Base.show(io::IO,a::FrameSpace)
    print(io,"Effspace:")
    print(io,a.space, "\n")
    print(io,"constraint:", a.con)
end

# getindex seem not very usefull in this case
getindex(space::FrameSpace ,i::Int) = @error "no efficient indexing, first use getspace function to generated the full space"

"""
    iterate(space::FrameSpace, state = 1)

Creates iterater object for an `FrameSpace`. Filters the not allowed constructs
"""
function Base.iterate(space::FrameSpace, state = 1 )
    # make costruct

    if  state > length(space.space)
        return
    end

    construct = space.space[state]

    #evaluated Constraints
    while filterconstraint(construct, space.con) == true
        state += 1
        if state > length(space.space)
            return
        end

        construct = space.space[state]
    end
    #update state
    state += 1

    return (construct,state)

end

"""
    length(space::FrameSpace)

Fallback to calculate the length FrameSpace.
Inefficient calculation has tot iterated overall allowed constructs.
"""

function Base.length(space::FrameSpace)
    #NOTE: not sure if you use `@warn` correctly, please check!
    @info "No efficient calculation iterate to all allowed options to get the length"
    len = 0
    for i in space
        len +=1
    end

    return len
end

"""
    eltype(K::FrameSpace)

Returns the type of `FrameSpace{T}` if they are collected.
"""
Base.eltype(K::FrameSpace{T}) where {T} = eltype(K.space)

struct MultiSpace{T} <: AbstractSpace{T}
    space::Array{T}
    MultiSpace(space::Array) = promote_type([typeof(i) for i in space]...) |> x -> new{x}(space)
end

Base.eltype(::MultiSpace{T}) where {T} = eltype(T)

"""
    length(space::MultiSpace)

Returns the number of constructs in the `space`.
"""
Base.length(space::MultiSpace) = map(x -> length(x),space.space) |> sum

"""
    nspace(space::MultiSpac)

Return the number of single space that are used to construct the whole `space`
"""

_nspace(space::MultiSpace) = length(space.space)

"""
    Base.size(space::MultiSpace)

Returns a tuple of two numbers, the first position contains the number of
constructs in the whole `space`. The second number is the number of single space
that are construct the entire design space
"""
Base.size(space::MultiSpace) = (length(space), _nspace(space))

# chained the signal iterators, can be done with Base.Iterators.flatten(space.space) but less consistent with other spaces.

"""
    Base.iterate(space::MultiSpace)

Returns a tuple of two numbers, the first position contains the number of
constructs in the `space`. The second number is the number of `sing` that are
used to construct the entire design space
"""
function Base.iterate(space::MultiSpace, state = [1 1])

    temp = iterate(space.space[state[1]],state[2])

    if temp == nothing
        state[1] += 1
        if (state[1] <= _nspace(space))
            temp = iterate(space.space[state[1]],1)
            state[2] = temp[2]
            return (temp[1],state)
        else
            return
        end

    else
        state[2] = temp[2]
        return (temp[1],state )
    end
end

# special length, prevent printing of the warning

function getindex(space::MultiSpace,i::Int)
    @assert typeintersect(eltype(space), FrameSpace) !=
                    eltype(space) "Framespace don't have efficient indexing"
    @assert length(space) >= i  "BoundsError index is to high "
    return  _getindex(space,i, 1)
end

function _getindex(space::MultiSpace,index::Int, state = 1)
    controle = length(space.space[state])
    if controle >= index
        @inbounds return space.space[state][index]
    else
        _getindex(space,index-controle,state+=1)
    end
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
@testset "emcee.jl" begin
    @testset "gdemo" begin
        Random.seed!(9876)

        n_samples = 1000
        n_walkers = 250

        spl = Emcee(n_walkers, 2.0)
        @test DynamicPPL.alg_str(Sampler(spl, gdemo_default)) == "Emcee"

        chain = sample(gdemo_default, spl, n_samples)
        check_gdemo(chain)
    end

    @testset "initial parameters" begin
        nwalkers = 250
        spl = Emcee(nwalkers, 2.0)

        # No initial parameters, with im- and explicit `init_params=nothing`
        Random.seed!(1234)
        chain1 = sample(gdemo_default, spl, 1)
        Random.seed!(1234)
        chain2 = sample(gdemo_default, spl, 1; init_params=nothing)
        @test Array(chain1) == Array(chain2)

        # Initial parameters have to be specified for every walker
        @test_throws ArgumentError sample(gdemo_default, spl, 1; init_params=[2.0, 1.0])

        # Initial parameters
        chain = sample(gdemo_default, spl, 1; init_params=fill([2.0, 1.0], nwalkers))
        @test chain[:s] == fill(2.0, nwalkers, 1)
        @test chain[:m] == fill(1.0, nwalkers, 1)
    end
end


text:
function pagerank(
    g::AbstractGraph{U},
    α=0.85,
    n=100::Integer,
    ϵ=1.0e-6
    ) where U <: Integer

    # indegree(g, v) is estimated run-time to iterate over inneighbors(g, v)
    partitions = LightGraphs.optimal_contiguous_partition(indegree(g), nthreads(), nv(g))

    # collect dangling nodes
    dangling_nodes = [v for v in vertices(g) if outdegree(g, v) == 0]
    N = Int(nv(g))
    # solution vector and temporary vector
    x = fill(1.0 / N, N)
    xlast = copy(x)
    # personalization vector
    p = fill(1.0 / N, N)
    # adjustment for leaf nodes in digraph
    dangling_weights = p
    for _ in 1:n
        dangling_sum = 0.0
        for v in dangling_nodes
            dangling_sum += x[v]
        end
        # flow from teleprotation
        for v in vertices(g)
            xlast[v] = (1 - α + α * dangling_sum) * p[v]
        end
        # flow from edges
        @threads for v_set in partitions
            for v in v_set
                for u in inneighbors(g, v)
                    xlast[v] += α * x[u] / outdegree(g, u)
                end
            end
        end

        # l1 change in solution convergence criterion
        err = 0.0
        for v in vertices(g)
            err += abs(xlast[v] - x[v])
            x[v] = xlast[v]
        end
        if (err < N * ϵ)
            return x
        end
    end
    error("Pagerank did not converge after $n iterations.") # TODO 0.7: change to InexactError with appropriate msg.
end


text:
include(string(bin_dir, "/RNXCMP/", wrapperjlfile))
include(string(bin_dir, "/GFZRNX/", wrapperjlfile))
include(string(bin_dir, "/rinex3_lht/", wrapperjlfile))


text:
function besselj0(x::BigFloat)
    x = abs(x)
    T = eltype(x)
    if iszero(x)
        return one(x)
    elseif isinf(x)
        return zero(x)
    end

    if x <= 2.0
        z = x * x
        p = z * z * evalpoly(z, J0_2N) / evalpoly(z, J0_2D)
        p -= 0.25 * z
        p += one(T)
        return p
    else
        xinv = inv(x)
        z  = xinv * xinv
        if xinv <= 0.25
            if xinv <= 0.125
                if xinv <= 0.0625
                    p = evalpoly(z, P16_IN) / evalpoly(z, P16_ID)
                    q = evalpoly(z, Q16_IN) / evalpoly(z, Q16_ID)
                else
                    p = evalpoly(z, P8_16N) / evalpoly(z, P8_16D)
                    q = evalpoly(z, Q8_16N) /  evalpoly(z, Q8_16D)
                end
            elseif xinv <= 0.1875
                p = evalpoly(z, P5_8N) / evalpoly(z, P5_8D)
	            q = evalpoly(z, Q5_8N) / evalpoly(z, Q5_8D)
            else
                p = evalpoly(z, P4_5N) / evalpoly(z, P4_5D)
	            q = evalpoly(z, Q4_5N) / evalpoly(z, Q4_5D)
            end
        elseif xinv <= 0.5
            if xinv <= 0.375
                if xinv <= 0.3125
                    p = evalpoly(z, P3r2_4N) / evalpoly(z, P3r2_4D)
	                q = evalpoly(z, Q3r2_4N) / evalpoly(z, Q3r2_4D)
                else
                    p = evalpoly(z, P2r7_3r2N) / evalpoly(z, P2r7_3r2D)
                    q = evalpoly(z, Q2r7_3r2N) / evalpoly(z, Q2r7_3r2D)
                end
            elseif xinv <= 0.4375
                p = evalpoly(z, P2r3_2r7N) / evalpoly(z, P2r3_2r7D)
	            q = evalpoly(z, Q2r3_2r7N) / evalpoly(z, Q2r3_2r7D)
            else
                p = evalpoly(z, P2_2r3N) / evalpoly(z, P2_2r3D)
	            q = evalpoly(z, Q2_2r3N) / evalpoly(z, Q2_2r3D)
            end
        end

        p = 1.0 + z * p
        q = z * xinv * q
        q = q - 0.125 * xinv
        c = cos(x)
        s = sin(x)
        ss = s - c
        cc = s + c
        z = -cos(x + x)
        if (s * c) < 0
            cc = z / ss;
        else
            ss = z / cc;
        end
        z = ONEOSQPI(T) * (p * cc - q * ss) / sqrt(x)
        return z
    end
end


text:
# This file is a part of BitOperations.jl, licensed under the MIT License (MIT).

__precompile__(true)

module BitOperations

include("bitops.jl")
include("endianess.jl")
include("zigzagzenc.jl")

end # module


text:
function liu_fit(xdata, ydata)
    # each term in Liu's time-domain impedance is denoted by m# (i.e. m1, m2, m3)
    m1(L) = 2*L/xdata[end]
    # m2 is complicated because of the sqrt(t). The units don't strictly work out
    # because there's supposed to be a 1/sqrt(Hz) in the denominator. Now, L is to
    # be computed nH and t should be provided as ns. So, scaling everything would
    # result in something like (Rs*1e-6)*sqrt.(t*1e-9)/(pi*L*1e-9). But, we don't
    # want Julia to have to work with big numbers if she doesn't have to. So, we'll
    # combine all of these exponentials into one: 10^-6*10^-4.5*10^9 = 10^-1.5
    m2(t,Rs,L) = 10^-1.5*Rs*sqrt.(t)/(pi*L) # Rs - microOhms, t - ns, L - nH
    m3(t,Rdc,L) = .5*Rdc*t/L # Rdc- ohms, t - ns, L - nH
    m(t,p) = m1(p[1])*(1+m2(t,p[2],p[1])+m3(t,p[3],p[1])) # p[1:3] -> L,Rs,Rdc
    fit = LsqFit.curve_fit(m,xdata,ydata,rand(3))
    errors = LsqFit.estimate_errors(fit)
    L = fit.param[1]
    C = xdata[end]^2/(4*L)
    Rs = fit.param[2]
    Rdc = fit.param[3]
    #@printf("L [nH]: %2.2f +/- %2.2f (%2.2f,%2.2f)\n", L, errors[1],
            #fit.param[1]-errors[1],fit.param[1]+errors[1])
    #@printf("C [μF]: %2.2f\n", C)
    #@printf("Rs [μΩ]: %2.2f +/- %2.2f (%2.2f,%2.2f)\n", Rs, errors[2],
            #fit.param[2]-errors[2],fit.param[2]+errors[2])
    #@printf("Rdc [Ω]: %2.2f +/- %2.2f (%2.2f,%2.2f)\n", Rdc, errors[3],
            #fit.param[3]-errors[3],fit.param[3]+errors[3])
    return Dict(
                "L" => L,
                "C" => C,
                "Rs" => Rs,
                "Rdc" => Rdc,
               )
end


text:
function lowrank(target, rank)
    f = svd(target)
    f.S[(rank+1):end] .= 0.0
    return f.U * Diagonal(f.S) * f.Vt
end

n1(Y) = fill(mean(Y), size(Y))
n2(Y) = (mean(Y; dims=2).+mean(Y; dims=1))./2.0
n3(Y) = (mean(Y; dims=2).+mean(Y; dims=1).+mean(Y))./3.0

function impute(query, template, position; rank=4)
    initial_value = query[position]
    query[position] = template[position]
    for iteration in 1:20
        _tmp = lowrank(query, rank)
        query[position] = _tmp[position]
    end
    prediction = (position=position, initial=template[position], updated=query[position])
    query[position] = initial_value
    return prediction
end

max_stars_repo_name: viralemergence/trefle

max_issues_repo_name: viralemergence/trefle

max_forks_repo_name: viralemergence/trefle