text:
using CSV, DataFrames
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()

all_options = Dict(
    "America" => ["New York City", "San Francisco", "Cincinnati"],
    "Canada" => ["Montreal", "Toronto", "Ottawa"],
)

app.layout = html_div() do
    html_div(
        children = [
            dcc_radioitems(
                id = "countries-radio",
                options = [(label = i, value = i) for i in keys(all_options)],
                value = "America",
            ),
            html_hr(),
            dcc_radioitems(id = "cities-radio"),
            html_hr(),
            html_div(id = "display-selected-values"),
        ],
    )
end

callback!(
    app,
    Output("cities-radio", "options"),
    Input("countries-radio", "value"),
) do selected_country
    return [(label = i, value = i) for i in all_options[selected_country]]
end

callback!(
    app,
    Output("cities-radio", "value"),
    Input("cities-radio", "options"),
) do available_options
    return available_options[1][:value]
end

callback!(
    app,
    Output("display-selected-values", "children"),
    Input("countries-radio", "value"),
    Input("cities-radio", "value"),
) do selected_country, selected_city
    return "$(selected_city) is a city in $(selected_country) "
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
#------------------------------------------------------------------------------------------------------------------------------#

abstract type AbstractRobot end

import HorizonSideRobots: move!, isborder, putmarker!, ismarker, temperature

#get(robot::AbstractRobot) = robot::Robot

move!(robot::AbstractRobot, side) = move!(get(robot), side)
isborder(robot::AbstractRobot, side) = isborder(get(robot), side)
putmarker!(robot::AbstractRobot) = putmarker!(get(robot))
ismarker(robot::AbstractRobot) = ismarker(get(robot))
temperature(robot::AbstractRobot) = temperature(get(robot))

#------------------------------------------------------------------------------------------------------------------------------#

mutable struct CountmarkersRobot <: AbstractRobot
    CmR::Robot
    counter::Integer
end

CmR = CountmarkersRobot(originRobot,0)

get(robot::CountmarkersRobot) = robot.CmR

function move!(robot::CountmarkersRobot, side::HorizonSide)
    move!(robot.CmR, side) # move!(get(CmR), side)
    if (ismarker(robot.CmR))
        robot.counter += 1
    end
end

#------------------------------------------------------------------------------------------------------------------------------#

mutable struct CoordsRobot <: AbstractRobot
    coordsrobot::Robot
    x::Int
    y::Int
end

coordsrobot = CoordsRobot(originRobot, 0, 0)

get(robot::CoordsRobot) = robot.coordsrobot

function move!(robot::CoordsRobot, side::HorizonSide)
    if side == Nord
        robot.y += 1
    elseif side == Sud
        robot.y -= 1
    elseif side == Ost
        robot.x += 1
    else
        robot.x -= 1
    end
    move!(robot.coordsrobot, side)
end

get_coords(robot::CoordsRobot) = (robot.x , robot.y)

#------------------------------------------------------------------------------------------------------------------------------#

mutable struct Beta <: AbstractRobot
    beta::Robot
    path::Vector
end

beta = Beta(originRobot, [])

get(robot::Beta) = robot.beta

function move!(robot::Beta, side::HorizonSide)
    push!(robot.path, Int(side))
    move!(robot.beta, side)
end

function get_back!(robot::Beta)
    for i in 1:length(robot.path)
        move!(robot.beta, anti_side(HorizonSide(robot.path[ length(robot.path) - i + 1 ])))
    end
    empty!(robot.path)
end

function try_move!(robot::Beta, side::HorizonSide)::Bool

    # если перед роботом есть перегородка, он начнет выполнять действия, чтобы попытаться её обойти, иначе просто пройдет в заданном направлении
    if isborder(robot.beta, side)
        steps = 0
        nextside = next_side(side)
        # робот пытается обойти перегородку со строны, которая следует за заданной по часовой стрелке
        while isborder(robot.beta, side)
            if !isborder(robot.beta, nextside)
                move!(robot, nextside)
                steps += 1
            else
                # срабатывает, когда робот уперся в угол. Он возвращается в то место, с которого начал попытку обойти по часовой стороне
                for j in 1:steps
                    move!(robot, anti_side(nextside))
                end
                return false
            end
        end

        # если у робота получилось обойти перегородку с первого раза, будет выполнено слудущее:
        move!(robot, side)
        while isborder(robot.beta, anti_side(nextside))
            move!(robot, side)
        end

        for j in 1:steps
            move!(robot, anti_side(nextside))
        end
        return true

    else
        # перегородки не было, робот прошел в заданном направлении
        move!(robot, side)
        return true
    end
end

#------------------------------------------------------------------------------------------------------------------------------#

mutable struct Gamma <: AbstractRobot
    gamma::Robot
end

gamma = Gamma(originRobot)

get(robot::Gamma) = robot.gamma

function try_move_var2!(robot::Robot, side::HorizonSide)::Bool
    steps = 0
    nextside = next_side_pr(side)
    # робот пытается обойти перегодку со стороны противоположной часовой стрелке
    while isborder(robot, side)
        if !isborder(robot, nextside)
            move!(robot, nextside)
            steps += 1
        else
            # если обойти невозможно вернется на исходное место
            for j in 1:steps
                move!(robot, anti_side(nextside))
            end
            # обойти прегородку невозможно
            return false
        end
    end

    # обход перегородки в ширину
    move!(robot, side)
    while isborder(robot, anti_side(nextside))
        move!(robot, side)
    end

    moving_defsteps!(robot, anti_side(nextside), steps)
    move!(robot,side)
    return true
end

function move!(robot::Gamma, side::HorizonSide)
    # если перед роботом есть перегородка, он начнет выполнять действия, чтобы попытаться её обойти, иначе просто пройдет в заданном направлении
    if isborder(robot.gamma, side)
        flag_nextside = true
        steps = 0
        nextside = next_side(side)
        # сначала робот пытается обойти перегородку со строны, которая следует за заданной по часовой стрелке
        while isborder(robot.gamma, side)
            if !isborder(robot.gamma, nextside)
                move!(robot.gamma, nextside)
                steps += 1
            else
                # срабатывает, когда робот уперся в угол. Он возвращается в то место, с которого начал попытку обойти по часовой стороне
                for j in 1:steps
                    move!(robot.gamma, anti_side(nextside))
                end

                flag_nextside = false

                # робот попытается обойти перегородку со сторны противоположной часовой стрелке
                return try_move_var2!(robot.gamma, side)
            end
        end

        # если у робота получилось обойти перегородку с первого раза, будет выполнено слудущее:
        if (flag_nextside)
            move!(robot.gamma, side)
            while isborder(robot.gamma, anti_side(nextside))
                move!(robot.gamma, side)
            end
            moving_defsteps!(robot.gamma, anti_side(nextside), steps)

            return true
        end
    else

        # перегородки не было, робот прошел в заданном направлении
        move!(robot.gamma,side)
        return true
    end
end

# moving_in_angle! передвинет робота в Юго-Западный угол поля и вернёт в переменную back_bath ОБРАТНЫЙ путь робота
function moving_in_angle!(robot::Gamma)::Vector{Int}
    back_path = []
    while (!(isborder(robot.gamma, West)) || !(isborder(robot.gamma, Sud)))
        for side in (Sud, West)
            while (!(isborder(robot.gamma,side)))
                move!(robot,side)
                push!(back_path, Int(side))
            end
        end
    end
    return reverse!(back_path)
end

# moving_back_to_start! вернёт робота в начальное положение через обратный путь робота
function moving_back_to_start!(robot::Gamma, back_path::Vector)::Nothing
    for i in 1:length(back_path)
        move!(robot, anti_side(HorizonSide(back_path[i])))
    end
end

#------------------------------------------------------------------------------------------------------------------------------#

struct Alpha{TypeRobot}
    alpha::TypeRobot
end

alpha = Alpha(originRobot)

get(robot::Alpha) = robot.alpha

#=
function try_move!(robot::PutmarkersRobot, side)::Bool
    result = try_move!(robot.robot, side)
    if result
        putmarker!(robot.robot) # предполагается, что для каждого возможного конкретного типа TypeRobot определен метод putmarker!
    end
    return result
end
=#

#------------------------------------------------------------------------------------------------------------------------------#

#=mutable struct CountmarkersAndTemperatureRobot <: CountmarkersRobot
    robot::CountmarkersRobot
    temperature::Real
end

getTemperature(robot::CountmarkersAndTemperatureRobot) = robot.temperature
getCountOfMarkers(robot::CountmarkersAndTemperatureRobot) = robot.CmR.counter

CmR = CountmarkersRobot(originRobot,0)
markerAndTemperature = CountmarkersAndTemperatureRobot(CmR, 0)

function move!(robot::CountmarkersAndTemperatureRobot)
    move!(markerAndTemperature.robot.CmR, side)
    if (ismarker(rmarkerAndTemperature.robot.CmR))
        markerAndTemperature.robot.counter += 1
        markerAndTemperature.temperature += temperature(markerAndTemperature.robot.CmR)
    end
end=#

mutable struct MODcountmarkersRobot <: AbstractRobot
    robot::Robot
    markers::Integer
    temperature::Real
end

mod_CmR = MODcountmarkersRobot(originRobot,0,0)

getTemperature(robot::MODcountmarkersRobot) = robot.temperature
getCountOfMarkers(robot::MODcountmarkersRobot) = robot.markers

function move!(robot::MODcountmarkersRobot, side::HorizonSide)
    if (ismarker(robot.robot))
        robot.markers += 1
        robot.temperature += temperature(robot.robot)
    end
    move!(robot.robot, side)
    if (ismarker(robot.robot))
        robot.markers += 1
        robot.temperature += temperature(robot.robot)
    end
end

ext: jl

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
module Roots

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using Printf
import CommonSolve
import CommonSolve: solve, solve!, init

export fzero,
       fzeros,
       secant_method

export find_zero, find_zero!, find_zeros,
       ZeroProblem, solve, solve!, init,
       Order0, Order1, Order2, Order5, Order8, Order16

export Bisection, FalsePosition

## load in files
include("utils.jl")
include("find_zero.jl")
include("bracketing.jl")
include("derivative_free.jl")
include("order0.jl")
include("simple.jl")
include("find_zeros.jl")
include("newton.jl")
include("lith.jl")
include("fzero.jl")

end


text:
using Test
using Distributions
using StatsFuns

using LRMoE

@testset "Poisson, exact" begin

    # # Fit 20,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(0.5, 8), 1)[1]

    # X = rand(Distributions.Uniform(-5, 5), 20000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [PoissonExpert(param1) PoissonExpert(2*param1);
    #         PoissonExpert(3*param1) PoissonExpert(0.5*param1)]

    # expos = rand(Distributions.Uniform(0.1, 5), 20000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [PoissonExpert(0.5*param1) PoissonExpert(3*param1);
    #                 PoissonExpert(2*param1) PoissonExpert(param1)]

    # pen_α = 5
    # pen_params = [[[1.0 Inf], [1.0 Inf]],
    #               [[1.0 Inf], [1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIPoisson, exact" begin

    # # Fit 20,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(0.5, 8), 1)[1]

    # X = rand(Distributions.Uniform(-5, 5), 20000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIPoissonExpert(param0[1], param1) ZIPoissonExpert(param0[2], 2*param1);
    #          ZIPoissonExpert(param0[3], 3*param1) ZIPoissonExpert(param0[4], 0.5*param1)]

    # expos = rand(Distributions.Uniform(0.1, 5), 20000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIPoissonExpert(0.50, 0.5*param1) ZIPoissonExpert(0.50, 3*param1);
    #                ZIPoissonExpert(0.50, 2*param1) ZIPoissonExpert(0.50, param1)]

    # pen_α = 5
    # pen_params = [[[1.0 Inf], [1.0 Inf]],
    #               [[1.0 Inf], [1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "NegativeBinomial, exact" begin

    # Fit 50,000 observations as an example

    # using Random
    # Random.seed!(1234)

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(0.5, 8), 4)
    # param2 = rand(Distributions.Uniform(0, 1), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [NegativeBinomialExpert(param1[1], param2[1]) NegativeBinomialExpert(param1[2], param2[2]);
    #          NegativeBinomialExpert(param1[3], param2[3]) NegativeBinomialExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [NegativeBinomialExpert(0.5*param1[1], 0.5) NegativeBinomialExpert(3*param1[2], 0.5);
    #                NegativeBinomialExpert(2*param1[3], 0.5) NegativeBinomialExpert(0.2*param1[4], 0.5)]

    # pen_α = 5
    # pen_params = [[[1.0 Inf], [1.0 Inf]],
    #               [[1.0 Inf], [1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZINegativeBinomial, exact" begin

    # Fit 50,000 observations as an example

    # param0 = [0.25 0.30 0.80 0.75] # rand(Distributions.Uniform(0, 1), 4)
    # param1 = [3 2.5 8 10] # rand(Distributions.Uniform(0.5, 8), 4)
    # param2 = [0.75 0.38 0.26 0.5]# rand(Distributions.Uniform(0, 1), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZINegativeBinomialExpert(param0[1], param1[1], param2[1]) ZINegativeBinomialExpert(param0[2], param1[2], param2[2]);
    #          ZINegativeBinomialExpert(param0[3], param1[3], param2[3]) ZINegativeBinomialExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZINegativeBinomialExpert(0.5, 0.5*param1[1], 0.5) ZINegativeBinomialExpert(0.5, 3*param1[2], 0.5);
    #                ZINegativeBinomialExpert(0.5, 2*param1[3], 0.5) ZINegativeBinomialExpert(0.5, 0.2*param1[4], 0.5)]

    # pen_α = 5
    # pen_params = [[[1.0 Inf], [1.0 Inf]],
    #               [[1.0 Inf], [1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Binomial, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(5:35, 4)
    # param2 = rand(Distributions.Uniform(0, 1), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [BinomialExpert(param1[1], param2[1]) BinomialExpert(param1[2], param2[2]);
    #          BinomialExpert(param1[3], param2[3]) BinomialExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [BinomialExpert(param1[1], 0.4) BinomialExpert(param1[2], 0.5);
    #                BinomialExpert(param1[3], 0.6) BinomialExpert(param1[4], 0.5)]

    # pen_α = 5
    # pen_params = [[[1.0 Inf], [1.0 Inf]],
    #               [[1.0 Inf], [1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIBinomial, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(5:35, 4)
    # param2 = rand(Distributions.Uniform(0, 1), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIBinomialExpert(param0[1], param1[1], param2[1]) ZIBinomialExpert(param0[2], param1[2], param2[2]);
    #          ZIBinomialExpert(param0[3], param1[3], param2[3]) ZIBinomialExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIBinomialExpert(0.5, param1[1], 0.4) ZIBinomialExpert(0.5, param1[2], 0.5);
    #                ZIBinomialExpert(0.5, param1[3], 0.6) ZIBinomialExpert(0.5, param1[4], 0.5)]

    # pen_α = 5
    # pen_params = [[[1.0 Inf], [1.0 Inf]],
    #               [[1.0 Inf], [1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "GammaCount, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(0.5, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 10), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [GammaCountExpert(param1[1], param2[1]) GammaCountExpert(param1[2], param2[2]);
    #          GammaCountExpert(param1[3], param2[3]) GammaCountExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [GammaCountExpert(0.5*param1[1], 2*param2[1]) GammaCountExpert(1.5*param1[2], 1.2*param2[2]);
    #                GammaCountExpert(param1[3], 0.8*param2[3]) GammaCountExpert(1.1*param1[4], 0.75*param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIGammaCount, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(0.5, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 10), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIGammaCountExpert(param0[1], param1[1], param2[1]) ZIGammaCountExpert(param0[2], param1[2], param2[2]);
    #          ZIGammaCountExpert(param0[3], param1[3], param2[3]) ZIGammaCountExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIGammaCountExpert(0.5, 0.5*param1[1], 2*param2[1]) ZIGammaCountExpert(0.5, 1.5*param1[2], 1.2*param2[2]);
    #                ZIGammaCountExpert(0.5, param1[3], 0.8*param2[3]) ZIGammaCountExpert(0.5, 1.1*param1[4], 0.75*param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "LogNormal, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(0.5, 5), 4)
    # param2 = rand(Distributions.Uniform(0.5, 3), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [LogNormalExpert(param1[1], param2[1]) LogNormalExpert(param1[2], param2[2]);
    #          LogNormalExpert(param1[3], param2[3]) LogNormalExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [LogNormalExpert(0.5*param1[1], 1.2*param2[1]) LogNormalExpert(0.75*param1[2], 0.8*param2[2]);
    #                LogNormalExpert(2*param1[3], 1.25*param2[3]) LogNormalExpert(0.1*param1[4], param2[4])]

    # pen_α = 5
    # pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf]],
    #               [[Inf 1.0 Inf], [Inf 1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "LogNormal, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(0.5, 5), 4)
    # param2 = rand(Distributions.Uniform(0.5, 3), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [LogNormalExpert(param1[1], param2[1]) LogNormalExpert(param1[2], param2[2]);
    #          LogNormalExpert(param1[3], param2[3]) LogNormalExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [LogNormalExpert(0.8*param1[1], 1.2*param2[1]) LogNormalExpert(0.75*param1[2], 1.3*param2[2]);
    #                LogNormalExpert(1.2*param1[3], 1.25*param2[3]) LogNormalExpert(0.9*param1[4], param2[4])]

    # pen_α = 5
    # pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf]],
    #               [[Inf 1.0 Inf], [Inf 1.0 Inf]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 0.5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZILogNormal, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(0.5, 5), 4)
    # param2 = rand(Distributions.Uniform(0.5, 3), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZILogNormalExpert(param0[1], param1[1], param2[1]) ZILogNormalExpert(param0[2], param1[2], param2[2]);
    #          ZILogNormalExpert(param0[3], param1[3], param2[3]) ZILogNormalExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZILogNormalExpert(0.5, 0.5*param1[1], 1.2*param2[1]) ZILogNormalExpert(0.5, 0.75*param1[2], 0.8*param2[2]);
    #                ZILogNormalExpert(0.5, 2*param1[3], 1.25*param2[3]) ZILogNormalExpert(0.5, 0.1*param1[4], param2[4])]

    # pen_α = 5
    # pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf]],
    #               [[Inf 1.0 Inf], [Inf 1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZILogNormal, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(0.5, 5), 4)
    # param2 = rand(Distributions.Uniform(0.5, 3), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZILogNormalExpert(param0[1], param1[1], param2[1]) ZILogNormalExpert(param0[2], param1[2], param2[2]);
    #          ZILogNormalExpert(param0[3], param1[3], param2[3]) ZILogNormalExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZILogNormalExpert(0.5, 0.5*param1[1], 1.2*param2[1]) ZILogNormalExpert(0.5, 0.75*param1[2], 0.8*param2[2]);
    #                ZILogNormalExpert(0.5, 2*param1[3], 1.25*param2[3]) ZILogNormalExpert(0.5, 0.1*param1[4], param2[4])]

    # pen_α = 5
    # pen_params = [[[Inf 1.0 Inf], [Inf 1.0 Inf]],
    #               [[Inf 1.0 Inf], [Inf 1.0 Inf]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 0.5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # result = LRMoE.fit_LRMoE(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     exact_Y = false,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Gamma, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(3, 8), 4)
    # param2 = rand(Distributions.Uniform(0.5, 15), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [GammaExpert(param1[1], param2[1]) GammaExpert(param1[2], param2[2]);
    #          GammaExpert(param1[3], param2[3]) GammaExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [GammaExpert(0.5*param1[1], 1.2*param2[1]) GammaExpert(0.8*param1[2], 0.8*param2[2]);
    #                GammaExpert(2*param1[3], 1.25*param2[3]) GammaExpert(0.9*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Gamma, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(3, 8), 4)
    # param2 = rand(Distributions.Uniform(0.5, 15), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [GammaExpert(param1[1], param2[1]) GammaExpert(param1[2], param2[2]);
    #          GammaExpert(param1[3], param2[3]) GammaExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [GammaExpert(0.5*param1[1], 1.2*param2[1]) GammaExpert(0.8*param1[2], 0.8*param2[2]);
    #                GammaExpert(2*param1[3], 1.25*param2[3]) GammaExpert(0.9*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 0.5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIGamma, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(3, 8), 4)
    # param2 = rand(Distributions.Uniform(0.5, 15), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIGammaExpert(param0[1], param1[1], param2[1]) ZIGammaExpert(param0[2], param1[2], param2[2]);
    #          ZIGammaExpert(param0[3], param1[3], param2[3]) ZIGammaExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIGammaExpert(0.5, 0.5*param1[1], 1.2*param2[1]) ZIGammaExpert(0.5, 0.8*param1[2], 0.8*param2[2]);
    #                ZIGammaExpert(0.5, 2*param1[3], 1.25*param2[3]) ZIGammaExpert(0.5, 0.9*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIGamma, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(3, 8), 4)
    # param2 = rand(Distributions.Uniform(0.5, 15), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIGammaExpert(param0[1], param1[1], param2[1]) ZIGammaExpert(param0[2], param1[2], param2[2]);
    #          ZIGammaExpert(param0[3], param1[3], param2[3]) ZIGammaExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIGammaExpert(0.5, 0.5*param1[1], 1.2*param2[1]) ZIGammaExpert(0.5, 0.8*param1[2], 0.8*param2[2]);
    #                ZIGammaExpert(0.5, 2*param1[3], 1.25*param2[3]) ZIGammaExpert(0.5, 0.9*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 0.5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "InverseGaussian, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(1, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [InverseGaussianExpert(param1[1], param2[1]) InverseGaussianExpert(param1[2], param2[2]);
    #          InverseGaussianExpert(param1[3], param2[3]) InverseGaussianExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [InverseGaussianExpert(0.5*param1[1], 1.2*param2[1]) InverseGaussianExpert(0.75*param1[2], 0.8*param2[2]);
    #                InverseGaussianExpert(2*param1[3], 1.25*param2[3]) InverseGaussianExpert(0.1*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "InverseGaussian, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(1, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [InverseGaussianExpert(param1[1], param2[1]) InverseGaussianExpert(param1[2], param2[2]);
    #          InverseGaussianExpert(param1[3], param2[3]) InverseGaussianExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [InverseGaussianExpert(1.5*param1[1], 1.2*param2[1]) InverseGaussianExpert(1.75*param1[2], 0.8*param2[2]);
    #                InverseGaussianExpert(2*param1[3], 1.25*param2[3]) InverseGaussianExpert(2*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIInverseGaussian, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(1, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIInverseGaussianExpert(param0[1], param1[1], param2[1]) ZIInverseGaussianExpert(param0[2], param1[2], param2[2]);
    #          ZIInverseGaussianExpert(param0[3], param1[3], param2[3]) ZIInverseGaussianExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIInverseGaussianExpert(0.5, 0.5*param1[1], 1.2*param2[1]) ZIInverseGaussianExpert(0.5, 0.75*param1[2], 0.8*param2[2]);
    #                ZIInverseGaussianExpert(0.5, 2*param1[3], 1.25*param2[3]) ZIInverseGaussianExpert(0.5, 0.1*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIInverseGaussian, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(2, 10), 4)
    # param2 = rand(Distributions.Uniform(2, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIInverseGaussianExpert(param0[1], param1[1], param2[1]) ZIInverseGaussianExpert(param0[2], param1[2], param2[2]);
    #          ZIInverseGaussianExpert(param0[3], param1[3], param2[3]) ZIInverseGaussianExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIInverseGaussianExpert(0.5, 1.5*param1[1], 1.2*param2[1]) ZIInverseGaussianExpert(0.5, 1.75*param1[2], 0.8*param2[2]);
    #                ZIInverseGaussianExpert(0.5, 2*param1[3], 1.25*param2[3]) ZIInverseGaussianExpert(0.5, 1.1*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Weibull, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(1.5, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [WeibullExpert(param1[1], param2[1]) WeibullExpert(param1[2], param2[2]);
    #          WeibullExpert(param1[3], param2[3]) WeibullExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [WeibullExpert(1.1*param1[1], 1.2*param2[1]) WeibullExpert(2.75*param1[2], 0.8*param2[2]);
    #                WeibullExpert(2*param1[3], 1.25*param2[3]) WeibullExpert(0.9*param1[4], param2[4])]

    # pen_α = 5
    # pen_params = [[[2.0 10.0 1.0 Inf], [2.0 10.0 1.0 Inf]],
    #               [[2.0 10.0 1.0 Inf], [2.0 10.0 1.0 Inf]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Weibull, exact" begin

    # Fit 50,000 observations as an example

    # using Random
    # Random.seed!(1234)

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(1.5, 10), 4)
    # param2 = rand(Distributions.Uniform(2, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [WeibullExpert(param1[1], param2[1]) WeibullExpert(param1[2], param2[2]);
    #          WeibullExpert(param1[3], param2[3]) WeibullExpert(param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [WeibullExpert(1.1*param1[1], 1.2*param2[1]) WeibullExpert(param1[2], 0.8*param2[2]);
    #                WeibullExpert(2*param1[3], 1.25*param2[3]) WeibullExpert(param1[4], param2[4])]

    # pen_α = 5
    # pen_params = [[[2.0 10.0 1.0 Inf], [2.0 10.0 1.0 Inf]],
    #               [[2.0 10.0 1.0 Inf], [2.0 10.0 1.0 Inf]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 5,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIWeibull, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(1.5, 10), 4)
    # param2 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIWeibullExpert(param0[1], param1[1], param2[1]) ZIWeibullExpert(param0[2], param1[2], param2[2]);
    #          ZIWeibullExpert(param0[3], param1[3], param2[3]) ZIWeibullExpert(param0[4], param1[4], param2[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIWeibullExpert(0.5, 1.1*param1[1], 1.2*param2[1]) ZIWeibullExpert(0.5, 2.75*param1[2], 0.8*param2[2]);
    #                ZIWeibullExpert(0.5, 2*param1[3], 1.25*param2[3]) ZIWeibullExpert(0.5, 0.9*param1[4], param2[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Burr, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(2, 5), 4)
    # param2 = rand(Distributions.Uniform(2, 5), 4)
    # param3 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [BurrExpert(param1[1], param2[1], param3[1]) BurrExpert(param1[2], param2[2], param3[2]);
    #          BurrExpert(param1[3], param2[3], param3[3]) BurrExpert(param1[4], param2[4], param3[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [BurrExpert(1.1*param1[1], 1.2*param2[1], 2*param3[1]) BurrExpert(2.75*param1[2], 0.8*param2[2], 3*param3[2]);
    #                BurrExpert(2*param1[3], 1.25*param2[3], 1*param3[3]) BurrExpert(0.9*param1[4], param2[4], 1.2*param3[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "Burr, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 1)[1]
    # param1 = rand(Distributions.Uniform(2, 5), 4)
    # param2 = rand(Distributions.Uniform(2, 5), 4)
    # param3 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [BurrExpert(param1[1], param2[1], param3[1]) BurrExpert(param1[2], param2[2], param3[2]);
    #          BurrExpert(param1[3], param2[3], param3[3]) BurrExpert(param1[4], param2[4], param3[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [BurrExpert(1.1*param1[1], 1.2*param2[1], 2*param3[1]) BurrExpert(2.75*param1[2], 0.8*param2[2], 3*param3[2]);
    #                BurrExpert(2*param1[3], 1.25*param2[3], 1*param3[3]) BurrExpert(0.9*param1[4], param2[4], 1.2*param3[4])]

    # pen_α = 5
    # pen_params = [[[2.0 10 2.0 10 2.0 10], [2.0 10 2.0 10 2.0 10]],
    #               [[2.0 10 2.0 10 2.0 10], [2.0 10 2.0 10 2.0 10]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 5, ecm_iter_max = 20,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIBurr, exact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(2, 5), 4)
    # param2 = rand(Distributions.Uniform(2, 5), 4)
    # param3 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIBurrExpert(param0[1], param1[1], param2[1], param3[1]) ZIBurrExpert(param0[3], param1[2], param2[2], param3[2]);
    #          ZIBurrExpert(param0[3], param1[3], param2[3], param3[3]) ZIBurrExpert(param0[4], param1[4], param2[4], param3[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIBurrExpert(0.5, 1.1*param1[1], 1.2*param2[1], 2*param3[1]) ZIBurrExpert(0.5, 2.75*param1[2], 0.8*param2[2], 3*param3[2]);
    #                ZIBurrExpert(0.5, 2*param1[3], 1.25*param2[3], 1*param3[3]) ZIBurrExpert(0.5, 0.9*param1[4], param2[4], 1.2*param3[4])]

    # pen_α = 5
    # 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]]]

    # result = LRMoE.fit_exact(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

@testset "ZIBurr, inexact" begin

    # Fit 50,000 observations as an example

    # param0 = rand(Distributions.Uniform(0, 1), 4)
    # param1 = rand(Distributions.Uniform(2, 5), 4)
    # param2 = rand(Distributions.Uniform(2, 5), 4)
    # param3 = rand(Distributions.Uniform(0.5, 30), 4)

    # X = rand(Distributions.Uniform(-5, 5), 50000, 7)
    # α = rand(Distributions.Uniform(-1, 1), 2, 7)
    # α[2,:] .= 0

    # model = [ZIBurrExpert(param0[1], param1[1], param2[1], param3[1]) ZIBurrExpert(param0[3], param1[2], param2[2], param3[2]);
    #          ZIBurrExpert(param0[3], param1[3], param2[3], param3[3]) ZIBurrExpert(param0[4], param1[4], param2[4], param3[4])]

    # expos = rand(Distributions.Uniform(0.1, 5), 50000)
    # Y_sim = LRMoE.sim_dataset(α, X, model, exposure = expos)

    # # All exact observations
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small censoring
    # Y_sim = hcat(fill(0.0, length(Y_sim[:,1])), 0.975.* Y_sim[:,1], 1.025.* Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), fill(0.0, length(Y_sim[:,2])), 0.975.* Y_sim[:,2], 1.025.* Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Small truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], Y_sim[:,1], Y_sim[:,1], fill(Inf, length(Y_sim[:,1])), 0.1 .* Y_sim[:,2], Y_sim[:,2], Y_sim[:,2], fill(Inf, length(Y_sim[:,2])))
    # # Censoring + Truncation
    # Y_sim = hcat(0.1 .* Y_sim[:,1], 0.95.* Y_sim[:,1], 1.25.* Y_sim[:,1], 10 .* Y_sim[:,1], 0.1 .* Y_sim[:,2], 0.95.* Y_sim[:,2], 1.25.* Y_sim[:,2], 10 .* Y_sim[:,2])

    # α_guess = fill(0.0, 2, 7)
    # model_guess = [ZIBurrExpert(0.5, 1.1*param1[1], 1.2*param2[1], 2*param3[1]) ZIBurrExpert(0.5, 1.75*param1[2], 0.8*param2[2], 3*param3[2]);
    #                ZIBurrExpert(0.5, 2*param1[3], 1.25*param2[3], 1*param3[3]) ZIBurrExpert(0.5, 0.9*param1[4], param2[4], 1.2*param3[4])]

    # pen_α = 5
    # pen_params = [[[2.0 10 2.0 10 2.0 10], [2.0 10 2.0 10 2.0 10]],
    #               [[2.0 10 2.0 10 2.0 10], [2.0 10 2.0 10 2.0 10]]]

    # result = LRMoE.fit_main(Y_sim, X, α_guess, model_guess,
    #                     exposure = expos,
    #                     penalty = true,
    #                     pen_α = pen_α, pen_params = pen_params,
    #                     ϵ = 5, ecm_iter_max = 20,
    #                     # ϵ = 1e-03, α_iter_max = 5.0, ecm_iter_max = 200,
    #                     # grad_jump = true, grad_seq = nothing,
    #                     print_steps = true)

    # α
    # result.α_fit

    # model
    # result.model_fit

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
module Custom

using PyCall
export load_model, unstructured_predict

function to_bytes(n::Integer; bigendian=true, len=sizeof(n))
    bytes = Array{UInt8}(undef, len)
    for byte in (bigendian ? (1:len) : reverse(1:len))
        bytes[byte] = n & 0xff
        n >>= 8
    end
    return bytes
end

function load_model(code_dir)
    return nothing
end

function score_unstructured(data; model = nothing, kwargs = nothing)
    println("Model: $model")
    println("Incoming content type params: $kwargs")
    println("Incoming data type $(typeof(data))")
    println("Incoming data: $data")
    println("Incoming query params: $(kwargs.query)")

    if haskey(kwargs.query, "ret_mode")
        ret_mode = kwargs.query["ret_mode"]
    else
        ret_mode = nothing
    end

    word_count = split( replace( data, "\n" => " "), " ")
    word_count = length(word_count)

    if ret_mode == "binary"
        ret_data = to_bytes(word_count)
        ret_kwargs = PyCall.PyDict( Dict("mimetype" => "application/octet-stream"))
        ret = PyCall.pybytes(ret_data), ret_kwargs
    else
        ret_kwargs = PyCall.PyDict( Dict("mimetype" => kwargs.mimetype))
        ret = string(word_count)
        #, ret_kwargs
    end
    println("outgoing")
    println(ret)
    return ret
end

end

max_stars_repo_name: andreakropp/datarobot-user-models

max_issues_repo_name: andreakropp/datarobot-user-models

max_forks_repo_name: andreakropp/datarobot-user-models


text:
abstract type SSM end

mutable struct StateSpaceModel <: SSM
    y::REALMAT # T x p (response matrix)
    F::REALVEC # d x 1 (design vector)
    G::REALMAT # d x d (evolution matrix)

    m::REALMAT # d x p (parameter means)
    P::REALMAT # d x d (parameter covariance)
    S::REALMAT # p x p ()

    β::Real # Discount factor for covariance matrix
    δ::Real # Discount factor for parameters
    Δ::REALMAT
    ν::Real # Degrees of freedom of student-t distribution (calculated from β)
    n::Real #
    k::Real # Ensures that Σ is a random walk

    function StateSpaceModel(y, F, G, m, P, S, β, δ)
        p = size(y, 2)
        d = size(m, 1)

        # Check dimensions of the output equation
        !(p == size(m, 2) == size(S, 1) == size(S, 2)) && throw(DimensionMismatch("y has $p columns, m has $(size(m, 2)) columns, and S is a symmetric $(size(S)) matrix."))

        # Check dimensions of the state equation
        !(d == length(F) == size(P, 1) == size(P, 2) == size(G, 1) == size(G, 2)) && throw(DimensionMismatch("F has length $(length(F)), m has $(size(m, 1)) rows and P is a symmetric $(size(P)) matrix."))

        (δ  > 0 && δ <= 1) || throw(ArgumentError("0 < δ ≤ 1 required (currently $δ)."))
        (β > 2/3 && β < 1) || throw(ArgumentError("$(2//3) < β < 1 required (currently $β)."))
        n = 1/(1 - β)
        ν = β*n
        k = (β - p*β + p)/(2β - p*β + p - 1)
        Δ = diagm(fill(1/δ, d))
        return new(y, F, G, m, P, S, β, δ, Δ, ν, n, k)
    end
end

# Constructors for local level
function LocalLevel(y::REALMAT, m::REALMAT, P::REALMAT, S::REALMAT, β::Real, δ::Real)
    F = [1]
    G = ones(1, 1)
    return StateSpaceModel(y, F, G, m, P, S, β, δ)
end

# Constructors for local level trend
function LocalLevelTrend(y::REALMAT, m::REALMAT, P::REALMAT, S::REALMAT, β::Real, δ::Real)
    F = [1, 0]
    G = [1 1;
         0 1]
    return StateSpaceModel(y, F, G, m, P, S, β, δ)
end

# Constructors for local level cycle
function LocalLevelCycle(y::REALMAT, m::REALMAT, P::REALMAT, S::REALMAT, β::Real, δ::Real)
    F = [1, 1, 0]
    G = [1 0 0;
         0 0 0;
         0 0 0]
    return StateSpaceModel(y, F, G, m, P, S, β, δ)
end


text:
N, faces = 10^6, 1:6

numSol = sum([iseven(i+j) for i in faces, j in faces]) / length(faces)^2
mcEst  = sum([iseven(rand(faces) + rand(faces)) for i in 1:N]) / N

println("Numerical solution = $numSol \nMonte Carlo estimate = $mcEst")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Test
using ProxyInterfaces

# @forward
# --------

# taken from https://github.com/MikeInnes/Lazy.jl/blob/0de1643f37555396d615e825fb17b73edb0b5939/test/runtests.jl

# define structs for @forward macro testing below
struct Foo112 end
struct Bar112 f::Foo112 end

@testset "Forward macro" begin
  play(x::Foo112; y) = y                    # uses keyword arg
  play(x::Foo112, z) = z                    # uses regular arg
  play(x::Foo112, z1, z2; y) = y + z1 + z2  # uses both

  @forward Bar112.f play                    # forward `play` function to field `f`

  let f = Foo112(), b = Bar112(f)
    @test play(f, y = 1) === play(b, y = 1)
    @test play(f, 2) === play(b, 2)
    @test play(f, 2, 3, y = 1) === play(b, 2, 3, y = 1)
  end
end

# ProxyInterfaces
# ---------------

@testset "iterator" begin
  include("iterator.jl")
end

@testset "indexing" begin
  include("indexing.jl")
end

@testset "array" begin
  include("array.jl")
end

@testset "dict" begin
  include("dict.jl")
end


text:
using RNG
using Base.Test: @test


text:
export CuPtr

type CuPtr{T}
    ptr::CUdeviceptr
    dev::Int
    id::Int

    function CuPtr(ptr, dev, id)
        p = new(ptr, dev, id)
        finalizer(p, free)
        p
    end
end

Base.convert(::Type{CUdeviceptr}, p::CuPtr) = p.ptr
Base.convert{T}(::Type{Ptr{T}}, p::CuPtr) = Ptr{T}(p.ptr)
Base.unsafe_convert(::Type{CUdeviceptr}, p::CuPtr) = p.ptr
Base.unsafe_convert{T}(::Type{Ptr{T}}, p::CuPtr) = Ptr{T}(p.ptr)

const memstep = 1024
#const alloc_ptrs = CUdeviceptr[]
const free_ptrs = Dict{Int,Vector{CUdeviceptr}}[]

function malloc(T::Type, n::Int)
    dev = getdevice()
    while length(free_ptrs) < dev+1
        push!(free_ptrs, Dict{Int,Vector{CUdeviceptr}}())
    end

    id = sizeof(T) * n ÷ memstep + 1
    ptrs = get!(free_ptrs[dev+1], id, CUdeviceptr[])
    if length(ptrs) == 0
        p = CUdeviceptr[0]
        cuMemAlloc(p, id*memstep)
        #push!(alloc_ptrs, p[1])
        CuPtr{T}(p[1], dev, id)
    else
        CuPtr{T}(pop!(ptrs), dev, id)
    end
end

free(p::CuPtr) = push!(free_ptrs[p.dev+1][p.id], p.ptr)

function devicegc()
    gc()
    _dev = getdevice()
    for dev = 0:ndevices()-1
        setdevice(dev)
        for (id,ptrs) in free_ptrs[dev+1]
            for p in ptrs
                cuMemFree(p)
            end
            empty!(ptrs)
        end
    end
    setdevice(_dev)
    run(`nvidia-smi`)
end


text:
# This file was generated, do not modify it. # hide
using MLJ
using PrettyPrinting
MLJ.color_off() # hide
X, y = @load_iris
@load DecisionTreeClassifier


text:
module test_millboard_dicts

using Millboard
using Test

Millboard.set_table_mode(:markdown)

board = Dict()
@test sprint(show, table(board)) == ""

board = Dict(:a=>:b, "a"=>"b")
@test sprint(show, table(board)) == """
| KEY | VALUE |
|-----|-------|
| :a  | :b    |
| a   | b     |"""

board = Dict("1x3"=>[1 2 3], "2x3"=>[1 2 3; 4 5 6], "3x1"=> [1; 2; 3])
@test sprint(show, table(board)) == """
| KEY | VALUE |
|-----|-------|
| 1x3 | 1 2 3 |
|     | 1 2 3 |
| 2x3 | 4 5 6 |
|     | 1     |
|     | 2     |
| 3x1 | 3     |"""

end # module test_millboard_dicts


text:
using Weave
using Literate

"""
    function hfun_menubar()

Return the html code to display the menu bar, with the table of contents and other book info.
"""
@delay function hfun_menubar()
    menu = globvar(:menu)
    isnothing(menu) && return ""
    page_numbering = globvar(:page_numbering) === true
    toc = build_toc(menu, page_numbering)
    getfield.(last.(toc), :title)
    io = IOBuffer()

    write(
        io,
        """

        """,
    )
    for entry in last.(toc)
        item =
            entry.filename === nothing ? entry.title :
            "$(entry.title)"
        write(
            io,
            """

$item

            """,
        )
    end
    write(
        io,
        """

        """,
    )

    return String(take!(io))
end

"""
    function hfun_get_title(section = nothing)

Return the title of the section prepended with the section number.

The section referred to is that in the path given by the argument `section` or that of the current page if `section === nothing`.
"""
@delay function hfun_get_title(section = nothing)
    menu = globvar(:menu)
    filename = section === nothing ? locvar(:fd_rpath) : "$(section[1]).md"
    title = pagevar(filename, :title)
    isnothing(title) && return filename
    isnothing(menu) && return pagevar(filename, :title)

    filename_noext = first(splitext(filename))

    page_numbering = globvar(:page_numbering) === true
    toc = build_toc(menu, page_numbering)
    entries = toc[map(x -> x.filename, last.(toc)).==filename_noext]
    return length(entries) == 1 ? only(last.(entries)).title : "Title not found"
end

"""
    function hfun_link_section(section)

Return a link to the given section, with the proper numbered title for the section.
"""
@delay function hfun_link_section(section)
    return """$(hfun_get_title(section))"""
end

"""
    function hfun_navigation()

Return the html code to display the navigation buttons on the top and/or bottom of each page.
"""
@delay function hfun_navigation()
    io = IOBuffer()

    menu = globvar(:menu)
    isnothing(menu) && return String(take!(io))

    filename = locvar(:fd_rpath)

    filename_noext = first(splitext(filename))

    page_numbering = globvar(:page_numbering) === true
    toc = build_toc(menu, page_numbering)

    prevnext = build_prevnext(toc)
    entry = prevnext[first.(prevnext).==filename_noext]
    length(entry) == 1 || return String(take!(io))
    prev, next = last(only(entry))

    if prev !== nothing || next !== nothing
        write(
            io,
            """

            """,
        )
    end

    if prev !== nothing
        prev_link = "/$prev"
        prev_title = only(last.(toc[getfield.(last.(toc), :filename).==prev])).title
        write(
            io,
            """

            $prev_title ←

            """,
        )
    end

    if next !== nothing
        next_link = "/$next"
        next_title = only(last.(toc[getfield.(last.(toc), :filename).==next])).title
        write(
            io,
            """

                → $next_title

            """,
        )
    end

    if prev !== nothing || next !== nothing
        write(
            io,
            """

            """,
        )
    end

    return String(take!(io))
end

"""
    function hfun_github_repo_link()

Return the html code to display the link to the github repo.
"""
function hfun_github_repo_link()
    github_repo = globvar(:github_repo)
    (github_repo === nothing || github_repo == "") && return ""

    io = IOBuffer()

    write(
        io,
        """

        """,
    )

    return String(take!(io))
end

"""
    function build_toc(menu, page_numbering = true, level = 1, pre = "")

Build the table of contents out of the given `menu` variable, in a recursive way.
"""
function build_toc(menu, page_numbering = true, level = 1, pre = "")
    toc = []
    isnothing(menu) && return toc
    i = 0
    for m in menu
        sec, subsecs = m isa Pair ? m : (m, nothing)
        this_page_numbering = page_numbering * !startswith(sec, '*')
        sec = lstrip(sec, '*')

        new_pre = this_page_numbering === false ? "" : "$pre$(i += 1)."

        filename = process_it(sec)
        filename_noext = filename !== nothing ? first(splitext(filename)) : nothing

        title = filename === nothing ? sec : pagevar(filename, :title)
        if this_page_numbering !== false
            title = "$new_pre $title"
        end

        push!(toc, sec => (filename = filename_noext, title = title, level = level))
        if subsecs !== nothing
            append!(toc, build_toc(subsecs, this_page_numbering, level + 1, new_pre))
        end
    end
    return toc
end

"""
    function build_prevnext(toc)

Build the list with the information about the previous and next pages of each page.
"""
function build_prevnext(toc)
    prevnext = []
    isnothing(toc) && return prevnext
    prev = nothing
    ftoc = filter(x -> x.filename !== nothing, last.(toc))
    for i = 1:length(ftoc)-1
        push!(prevnext, ftoc[i].filename => (prev = prev, next = ftoc[i+1].filename))
        prev = ftoc[i].filename
    end
    push!(prevnext, ftoc[end].filename => (prev = prev, next = nothing))
    return prevnext
end

"""
    function process_it(filename)

Process a given julia script or Juno's style markdown file `.jmd` or jupyter notebook.

Processing is done via either Weave or Literate, depending on the folder it belongs to
(whether `_weave` or `_jupyter`, for Weave, or `_literate`, for Literate)

Processing generates a Markdown file for Franklin and a Jupyter notebook for the reader, if desired.
"""
function process_it(filename)
    startswith(filename, "pages/") && return "$filename.md"
    startswith(filename, "src/") || return nothing

    out_path = replace(
        dirname(filename),
        r"^src/weave" => "pages/weaved",
        r"^src/literate" => "pages/literated",
        r"^src/jupyter" => "pages/jupytered",
    )
    processed_filename = first(splitext("$out_path/$(basename(filename))")) * ".md"

    link_download_notebook = globvar(:link_download_notebook)
    link_nbview_notebook = globvar(:link_nbview_notebook)
    link_binder_notebook = globvar(:link_binder_notebook)
    exec_notebook = globvar(:exec_notebook)

    if mtime(filename) > mtime(processed_filename)
        if startswith(filename, "src/literate/")
            Literate.markdown(
                filename,
                out_path,
                execute = true,
                credit = false,
                flavor = Literate.FranklinFlavor(),
            )
        else
            Weave.weave(filename; out_path, fig_path = "images", doctype = "github")
        end
        postprocess_md(processed_filename, out_path)
    end

    notebook_output_dir = "generated/$(replace(dirname(filename), r"^src/weave" => "weaved", r"^src/literate" => "literated", r"^src/jupyter" => "jupytered"))"
    notebook_path = first(splitext("$notebook_output_dir/$(basename(filename))")) * ".ipynb"

    if any(
        ==(true),
        (link_download_notebook, link_nbview_notebook, link_binder_notebook),
    ) && mtime(filename) > mtime(notebook_path)

        if startswith(filename, "src/weave/")
            mkpath(notebook_output_dir)
            Weave.notebook(
                filename,
                out_path = notebook_output_dir,
                nbconvert_options = "--allow-errors",
            )
        elseif startswith(filename, "src/literate/")
            Literate.notebook(
                filename,
                notebook_output_dir,
                execute = exec_notebook,
                credit = false,
            )
        elseif startswith(filename, "src/jupyter/")
            mkpath(notebook_output_dir)
            cp(filename, "$notebook_output_dir/$(basename(filename))", force = true)
            #=
            # The following would have the advantage of forcing a clean
            # execution of the notebook, but Weave currently has some issues
            # with that
            Weave.notebook(
                filename,
                out_path = notebook_output_dir,
                nbconvert_options = "--allow-errors"
            ) =#
        end
        postprocess_ipynb(notebook_path)
    end

    return first(splitext(processed_filename)) # remove extension ".md"
end

"""
    function hfun_linkbadges()

Return the html code for displaying the badges for accessing the notebooks and the source
code associated with each page that has been processed via Weave or Literate.
"""
@delay function hfun_linkbadges()
    filename = locvar(:fd_rpath)

    link_view_source = globvar(:link_view_source)
    link_download_notebook = globvar(:link_download_notebook)
    link_nbview_notebook = globvar(:link_nbview_notebook)
    link_binder_notebook = globvar(:link_binder_notebook)

    github_repo = globvar(:github_repo)

    notebook_path = "generated/$(replace(dirname(filename), r"^pages/?" => ""))/$(replace(basename(filename), r".md$" => ".ipynb"))"

    io = IOBuffer()
    if (
        link_view_source == true &&
        startswith(filename, r"pages/(?:weaved|literated|jupytered)/")
    ) || (
        any(
            ==(true),
            (link_download_notebook, link_nbview_notebook, link_binder_notebook),
        ) && (isfile("$notebook_path"))
    )
        write(
            io,
            """

            """,
        )
        if link_nbview_notebook == true && isfile("$notebook_path")
            website = globvar(:website)
            write(
                io,
                """

                """,
            )
        end
        if link_binder_notebook == true && isfile("$notebook_path")
            nbgitpuller_repo = globvar(:nbgitpuller_repo)
            nbgitpuller_branch = globvar(:nbgitpuller_branch)
            binder_application = globvar(:binder_application)
            prepath = globvar(:prepath)
            if all(
                !isnothing,
                (nbgitpuller_repo, nbgitpuller_branch, binder_application, prepath),
            )
                nbgitpuller_repo = replace(nbgitpuller_repo, "https://github.com/" => "")

                binder_application =
                    binder_application == "" ? "tree" :
                    binder_application in ("lab", "retro") ? "$binder_application/tree" :
                    binder_application # not sure how the others work

                prepath == "" || (prepath = "$prepath/")
                write(
                    io,
                    """

                    """,
                )
            end
        end

        if link_download_notebook == true && isfile("$notebook_path")
            write(
                io,
                """

                """,
            )
        end

        if link_view_source == true &&
           startswith(filename, r"pages/(?:weaved|literated|jupytered)/")
            menu = globvar(:menu)
            page_numbering = globvar(:page_numbering) === true
            toc = build_toc(menu, page_numbering)
            ftoc = filter(x -> last(x).filename == first(splitext(filename)), toc)
            source_path = length(ftoc) > 0 ? "/$(first(first.(ftoc)))" : nothing
            write(
                io,
                """

                """,
            )
        end

        write(
            io,
            """

            """,
        )
    end
    return String(take!(io))
end

"""
    function postprocess_md(filename, out_path)

Post-process a markdown file generated by Weave or Literate for compatibility with
Franklin, by replacing `![...](...)`` with either a full path or with `\figalt`/`\fig`
depending on the situation; by moving any output image to the proper place; and by
properly handling the title, adding the section number, if needed.
"""
function postprocess_md(filename, out_path)
    tmppath, tmpio = mktemp()
    no_match_so_far = true
    open(filename, "r") do io
        for line in eachline(io, keep = true)
            if (m = match(r"^#\s+(.*)$", line)) !== nothing && no_match_so_far == true
                line = "@def title = \"$(m.captures[1])\"\n\n# {{ get_title }}\n"
                no_match_so_far = false
            else
                line = replace(
                    line,
                    r"!\[([^\]]*)\]\((?:\.\./)*_assets/([^\)]*)\)" =>
                        s"![\1](/assets/\2)",
                    r"^!\[\]\(([^/)][^\)]*)\)" => s"\\fig{\1}",
                    r"^!\[([^\]]*)\]\(([^/)][^\)]*)\)" => s"\\figalt{\1}{\2}",
                    r"(?<!\!)\[[^\]]*\]\(/pages/([^\)]*)\)" =>
                        s"{{link_section pages/\1}}",
                )
            end
            write(tmpio, line)
        end
    end
    close(tmpio)
    mv(tmppath, filename, force = true)
    if startswith(filename, r"pages/(?:weaved|jupytered)/") && isdir("$out_path/images/")
        destination = "_assets/$(first(splitext(filename)))/code/images"
        mkpath(destination)
        mv("$out_path/images", destination, force = true)
    elseif startswith(filename, "pages/literated/")
        destination = "_assets/$(first(splitext(filename)))/code"
        mkpath(destination)
        for file in filter(f -> !endswith(f, r"\.md"), readdir(out_path))
            mv("$out_path/$file", "$destination/$file", force = true)
        end
    end
end

function postprocess_ipynb(filename)
    tmppath, tmpio = mktemp()
    prepend = "../"^count(==('/'), filename)
    open(filename, "r") do io
        for line in eachline(io, keep = true)
            line = replace(
                line,
                r"(?<!')!\[([^\]]*)\]\(((?:\.\./)*)_assets/([^\)]*)\)" =>
                    s"![\1](\2assets/\3)",
                r"(?<!')\"(?:\.\./)*_assets/([^\)]*)\"" => SubstitutionString("\"$(prepend)assets/\\1\""),
                r"(?<!')(?<!\!)\[([^\]]*)\]\(/pages/([^\)]*)\)" =>
                    SubstitutionString("[\\1]($(prepend)pages/\\1)"),
            )
            write(tmpio, line)
        end
    end
    close(tmpio)
    mv(tmppath, filename, force = true)
end


text:
gt = PlotlyBase.GenericTrace("scatter"; x=1:10, y=(1:10).^2)
pplot = PlotlyBase.Plot(gt)
layout = PlotlyBase.Layout()

using PlotlyBase.JSON, JSON3

@testset "Convert to json" begin
    @test JSON.json(gt) == "{\"y\":[1,4,9,16,25,36,49,64,81,100],\"type\":\"scatter\",\"x\":[1,2,3,4,5,6,7,8,9,10]}"
    @test JSON.json(gt) == JSON3.write(gt)
    @test JSON.json(pplot) == JSON3.write(pplot)
    @test JSON.json(layout) == JSON3.write(layout)
end


text:
################################################################################
# Assign Model and Target Parameter
_modelProb = ModelWrapper(ProbModel(), val_dist)
_syms = [
    keys(_modelProb.val),
    keys(_modelProb.val)[1],
    keys(_modelProb.val)[end],
    (:d_a1, :d_a2, :d_a3, :d_a5, :d_a6, :d_a7, :d_a8),
]
_targets = [Tagged(_modelProb, _syms[iter]) for iter in eachindex(_syms)]
_params = [sample(_modelProb, _targets[iter]) for iter in eachindex(_syms)]

@testset "Tagged - Model parameter" begin
        ## Assign Sub Model
        for iter in eachindex(_syms)
        _sym = _syms[iter]
        _target = _targets[iter]
        _param = _params[iter]
        _model_temp = ModelWrapper(subset(val_dist, _sym))
        ## Compute logdensities and check length
        @test length(_model_temp) == length(_target)
        @test abs(log_prior(_model_temp, _target) - log_prior(_model_temp)) <= _TOL
        @test abs(log_abs_det_jac(_model_temp, _target) - log_abs_det_jac(_model_temp)) <=
              _TOL
        ## Type Check - Unflatten/Constrain
        theta_unconstrained_vec = randn(length(_target))
        theta_constrained = unflatten_constrain(_model_temp, theta_unconstrained_vec)
        theta_constrained2 = unflatten_constrain(
            _model_temp, _target, theta_unconstrained_vec
        )
        @test typeof(theta_constrained) == typeof(_model_temp.val)
        @test typeof(theta_constrained2) == typeof(_model_temp.val)
        _θ1  = flatten(_model_temp.info.reconstruct, theta_constrained)
        _θ2  = flatten(_target.info.reconstruct, theta_constrained2)
        @test sum(abs.(_θ1 - _θ2)) ≈ 0 atol = _TOL
        ## Utility functions
        log_prior(_target, _model_temp.val)
        θ_flat = flatten(_model_temp, _target)
        unflatten(_model_temp, _target, θ_flat)
        unflatten!(_model_temp, _target, θ_flat)
        θ_flat_unconstrained = unconstrain_flatten(_model_temp, _target)
        unflatten_constrain!(_model_temp, _target, θ_flat_unconstrained)
        log_prior_with_transform(_model_temp, _target)

        subset(_model_temp, _target)
        ModelWrappers.length(_target)
        ModelWrappers.paramnames(_target)
        fill(_model_temp, _target, _model_temp.val)
        fill!(_model_temp, _target, _model_temp.val)
        _model_temp.val
        sample(_RNG, _model_temp, _target)
        sample(_model_temp, _target)
        sample!(_RNG, _model_temp, _target)
        sample!(_model_temp, _target)
    end
end


text:
if !isdefined(Main, :Groebner)
    import Groebner
end

include("MQ/parser.jl")

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

BenchmarkTools.DEFAULT_PARAMETERS.samples = 5

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    @btime gb = Groebner.groebner($system)
    # gb = Groebner.groebner(system, reduced=false)
end

function run_f4_MQ_benchmarks()
    systems = [
        ("mq_10_20_2_0", read_MQ_GF("mq_n10_m20_p2_s0")),
        ("mq_10_20_2_1", read_MQ_GF("mq_n10_m20_p2_s1"))
    ]

    println()
    for (name, system) in systems
        println("$name")
        # benchmark_system_my(system)
        benchmark_system_my(system)
    end
end

run_f4_MQ_benchmarks()

#=

=#


text:
# entities immediately available as private members in IIF, but require Flux for actual use

struct FluxModelsDistribution{ID,OD,P,D<:AbstractArray}
  # shape of the input data
  inputDim::NTuple{ID,Int}
  # shape of the output data
  outputDim::NTuple{OD,Int}
  # actual Flux models
  models::Vector{P}
  # the data used for prediction, must be <: AbstractArray
  data::D
  # shuffle model predictions relative to particle index at each sampling
  shuffle::Base.RefValue{Bool}
  # false for default serialization with model info, set true for separate storage of models
  serializeHollow::Base.RefValue{Bool}
    # # TODO remove requirement and standardize sampler API
    # specialSampler::Function
end


text:
using Lazy

include("utils.jl")

fibs = @lazy BigInt(0):BigInt(1):(fibs + drop(BigInt(1), fibs));

fibonacci(n) = (n, fibs) |> take |> collect

fibonaccinum(n) = fibonacci(n)[end]

factor(n, base) = begin
    counter = 0
    if base == 0 || base == 1 || n == 0 return 0 end
    val = n
    while val % base == 0
        counter += 1
        val = val ÷ base
    end
    return counter
end

max_stars_repo_name: manuelelucchi/Fibonacci

max_issues_repo_name: manuelelucchi/Fibonacci

max_forks_repo_name: manuelelucchi/Fibonacci


text:
using Revise # reduce recompilation time
using Plots
using LinearAlgebra
using StaticArrays
using SparseArrays
using BenchmarkTools
using UnPack
using ForwardDiff

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DTri
using UniformTriMesh

using SetupDG
using ExplicitJacobians

"Approximation parameters"
N = 2 # The order of approximation
K1D = 8
CFL = 1
T = 1 # endtime

"Mesh related variables"
VX, VY, EToV = uniform_tri_mesh(K1D)
# VX = @. VX - .3*sin(pi*VX)

# initialize ref element and mesh
rd = init_reference_tri(N)
md = init_mesh((VX,VY),EToV,rd)

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,FToF,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps!(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB,FToF)
mapP[mapB] = mapPB

println("Done initializing mesh + ref elements")

## construct hybridized SBP operators

function build_SBP_ops(rd::RefElemData)

        @unpack M,Dr,Ds,Vq,Pq,Vf,wf,nrJ,nsJ = rd
        Qr = Pq'*M*Dr*Pq
        Qs = Pq'*M*Ds*Pq
        Ef = Vf*Pq
        Br = diagm(wf.*nrJ)
        Bs = diagm(wf.*nsJ)
        Qrh = .5*[Qr-Qr' Ef'*Br;
                -Br*Ef  Br]
        Qsh = .5*[Qs-Qs' Ef'*Bs;
                -Bs*Ef  Bs]

        Vh = [Vq;Vf]
        Ph = M\transpose(Vh)
        VhP = Vh*Pq

        # make skew symmetric versions of the operators"
        Qrhskew = .5*(Qrh-transpose(Qrh))
        Qshskew = .5*(Qsh-transpose(Qsh))

        return Qrhskew,Qshskew,Vh,Ph,VhP,M
end

# interpolate geofacs to both vol/surf nodes
@unpack Vq,Vf = rd
@unpack rxJ, sxJ, ryJ, syJ = md
rxJ, sxJ, ryJ, syJ = (x->[Vq;Vf]*x).((rxJ, sxJ, ryJ, syJ)) # interp to hybridized points

## global matrices

function build_global_ops(rd::RefElemData,md::MeshData)

        Qrhskew,Qshskew,Vh,Ph,VhP,M = build_SBP_ops(rd)

        Ax,Ay,Bx,By,B = assemble_global_SBP_matrices_2D(rd,md,Qrhskew,Qshskew)

        # add off-diagonal couplings
        Ax += Bx
        Ay += By

        Ax *= 2 # for flux differencing
        Ay *= 2

        AxTr = sparse(transpose(Ax))
        AyTr = sparse(transpose(Ay))
        Bx   = abs.(Bx) # create LF penalization term

        # globalize operators and nodes
        @unpack x,y,J = md
        VhTr = kron(speye(K),sparse(transpose(Vh)))
        Vh   = kron(speye(K),sparse(Vh))
        M    = kron(spdiagm(0 => J[1,:]),sparse(M))
        Ph   = kron(spdiagm(0 => 1 ./ J[1,:]), sparse(Ph))
        x,y = (a->a[:]).((x,y))

        return Ax,Ay,AxTr,AyTr,Bx,B,VhTr,Vh,M,Ph,x,y
end

Ax,Ay,AxTr,AyTr,Bx,B,VhTr,Vh,M,Ph,x,y =
        build_global_ops(rd::RefElemData,md::MeshData)
println("Done building global ops")

# jacspy = I+Ax+Ay
# jac = droptol!(jacspy,1e-12)
# jac = droptol!(VhTr*jacspy*Vh,1e-12)
# display(spy(jac .!= 0,ms=2.25))
# # title="nnz = $(nnz(jac))",
# error("d")

## define Burgers fluxes
function F(uL,uR)
        Fx = @. (uL^2 + uL*uR + uR^2)/6
        Fy = @. 0*uL
        return Fx,Fy
end

# extract coordinate fluxes
Fx = (uL,uR)->F(uL,uR)[1]
Fy = (uL,uR)->F(uL,uR)[2]

# AD for jacobians
dFx(uL,uR) = ForwardDiff.jacobian(uR->F(uL,uR)[1],uR)
dFy(uL,uR) = ForwardDiff.jacobian(uR->F(uL,uR)[2],uR)

## nonlinear solver stuff
function init_newton_fxn(Q,ops,rd::RefElemData,md::MeshData,funs,dt)

        Ax,Ay,AxTr,AyTr,Bx,B,M,Vh,Ph = ops
        Fx,Fy,dFx,dFy,LF,dLF = funs

        # set up normals for use in penalty term
        @unpack Vq,Vf = rd
        @unpack nxJ,nyJ,sJ = md
        Nq,Np = size(Vq)
        Nf    = size(Vf,1)
        Nh    = Nq + Nf
        fids  = (Nq+1):Nh
        nxhmat,nyhmat = ntuple(x->zeros(Nh,K),2)
        nxhmat[fids,:] = nxJ[:]./sJ[:]
        nyhmat[fids,:] = nyJ[:]./sJ[:]
        nxh,nyh = (x->x[:]).((nxhmat,nyhmat))

        # get lengths of arrays
        Nfields = length(Q)
        Id_fields = speye(Nfields) # for Kronecker expansion to large matrices - fix later with lazy evals
        Vh_fields = droptol!(kron(Id_fields,Vh),1e-12)
        M_fields  = droptol!(kron(Id_fields,M),1e-12)

        # init jacobian matrix
        dFdU_h = repeat(I+Ax+Ay,Nfields,Nfields)

        function midpt_newton_iter!(Qnew, Qprev) # for Burgers' eqn specifically

                Qh    = SVector((x->Vh*x).(Qnew)) # tuples are faster, but need arrays for ForwardDiff

                ftmp  = hadamard_sum(AxTr,Fx,Qh) + hadamard_sum(AyTr,Fy,Qh) + hadamard_sum(B,LF,Qh,nxh,nyh)
                f     = vcat((x->Ph*x).(ftmp)...)
                res   = vcat(Qnew...) + .5*dt*f - vcat(Qprev...)

                fill!(dFdU_h.nzval,0.0)
                accum_hadamard_jacobian!(dFdU_h, Ax, dFx, Qh)
                accum_hadamard_jacobian!(dFdU_h, Ay, dFy, Qh)
                accum_hadamard_jacobian!(dFdU_h, B, dLF, Qh,nxh,nyh) # flux term involving normals
                dFdU = droptol!(transpose(Vh_fields)*(dFdU_h*Vh_fields),1e-12)

                # solve and update
                dQ   = (M_fields + .5*dt*dFdU)\(M_fields*res)
                Qtmp = reshape(vcat(Qnew...) - dQ,length(Q[1]),Nfields)       # convert Qnew to column vector for update
                for fld = 1:Nfields
                        Qnew[fld] .= Qtmp[:,fld]
                end
                return norm(dQ) # convert back to array of arrays
        end
        return midpt_newton_iter!
end

## init condition, rhs

u = @. -sin(pi*x)
# u = randn(size(x))
Q = [u]
Qnew = deepcopy(Q)

# convert to tuple
Q = tuple(Q...)
Qnew = tuple(Qnew...)

# set time-stepping constants
CN = (N+1)*(N+2)/2  # estimated trace constant
h = minimum(md.J[1,:]./md.sJ[1,:]) # ratio of J/Jf = O(h^d/h^d-1)
dt = CFL * 2 * h / CN
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

function LF(uL,uR,nxL,nyL,nxR,nyR)
        absnx = @. (abs(nxL) + abs(nxR))/2
        return (@. max(abs(uL),abs(uR))*(uL-uR)*absnx)
end
dLF(uL,uR,args...) = ForwardDiff.jacobian(uR->LF(uL,uR,args...),uR)

# initialize Newton vars
ops = (Ax,Ay,copy(transpose(Ax)),copy(transpose(Ay)),Bx,B,M,Vh,Ph) # pack inputs together
fxns = (Fx,Fy,dFx,dFy,LF,dLF)
midpt_newton_iter! = init_newton_fxn(Q,ops,rd,md,fxns,dt)

## newton time iteration

bcopy!(utarget,u) = utarget .= u

it_count = zeros(Nsteps)
energy   = zeros(Nsteps)
for i = 1:Nsteps
        global Q,Qnew

        #Qnew = copy(Q)  # copy / over-write at each timestep
        bcopy!.(Qnew,Q)
        iter = 0
        dQnorm = 1
        while dQnorm > 1e-12
                dQnorm = midpt_newton_iter!(Qnew,Q)
                iter += 1
                if iter > 10
                        println("iter = $iter")
                end
        end
        it_count[i] = iter
        bcopy!.(Q, (@. 2*Qnew - Q)) # implicit midpoint rule

        u = Q[1]
        energy[i] = u'*M*u

        if i%10==0 || i==Nsteps
                println("Number of time steps $i out of $Nsteps")
                # display(scatter(x,Q[1]))
        end
end

@unpack VDM = rd
rp, sp = Basis2DTri.equi_nodes_2D(25)
Vp = Basis2DTri.vandermonde_2D(N,rp,sp)/VDM

xp,yp,up = (x->Vp*reshape(x,size(Vp,2),K)).((x,y,Q[1]))
gr(aspect_ratio=1, legend=false,
   markerstrokewidth=0, markersize=2)
# display(scatter(xp,yp,up,zcolor=up,cam=(3,25),axis=false))
scatter(xp,yp,up,zcolor=up,cam=(0,90),border=:none,axis=false)
# contourf(xp,yp,up,lw=1,fill=true,levels=20)
# png("sol_unif_mesh.png")

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
# # "Point exponential decay" time-stepper convergence test

using PyPlot

# Define a few utilities for running tests and unpacking and plotting results

include("ConvergenceTests/ConvergenceTests.jl") # we use the GaussianAdvectionDiffusion module here.

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

run_single_test = ConvergenceTests.PointExponentialDecay.run_test

""" Run advection-diffusion test for all Nx in resolutions. """
function run_convergence_test(stop_time, proposal_Δt...)

    # Adjust time-step
    stop_iterations = [round(Int, stop_time / dt) for dt in proposal_Δt]
                 Δt = [stop_time / stop_iter for stop_iter in stop_iterations]

    # Run the tests
    results = [run_single_test(Δt=dt, stop_iteration=stop_iter)
               for (dt, stop_iter) in zip(Δt, stop_iterations)]

    return results, Δt
end

unpack_errors(results) = map(r -> r.L₁, results)

stop_time = 3
Δt = vcat(1 ./ collect(1000:-100:100), 1 ./ collect(90:-10:10), 1 ./ collect(9:-1:1))
results, Δt = run_convergence_test(stop_time, Δt...)
L₁ = unpack_errors(results)

fig, axs = subplots()

loglog(Δt, L₁, "o--", alpha=0.8, linewidth=1, label="Measured error")

# Guide line to confirm second-order scaling
loglog(Δt, L₁[1] .* Δt / Δt[1], "k-", alpha=0.4, linewidth=2, label=L"\sim \Delta t")

xlabel(L"\Delta t")
ylabel("\$L_1\$-norm of \$ | c_\\mathrm{simulated}(t) - \\mathrm{e}^{-t} | \$")
title("Oceananigans time-stepper convergence for \$ c(t) = \\mathrm{e}^{-t} \$")
removespines("top", "right")
legend()

savefig("figs/point_exponential_decay_time_stepper_convergence.png", dpi=480)


text:
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]
A = m["A"]
C = m["C"]
Q = m["Q"]
R = m["R"]
x10 = m["m1"]
P10 = m["P1"]

###########################
print("Testing ACQR_ksmoother: ")
xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = ACQR_kfilter(y,A,C,Q,R,x10,P10)
xtN,PtN,Pt1tN = ACQR_ksmoother(A,xtt,Ptt,xtt1,Ptt1)
if sum(xtN.^2) > 1e-6
    println("OK")
end