text:
# Wrapper function(s) to simulate from a Count Data Model
# Only check if number of parameters matches length of θ

"""
    simulate(T, model, θ; burnin, pinfirst)
Function to generate time series from Count Data model.

* `T`: Length of time series
* `model`: Model specification
* `θ`: Parameters (vector or struct)
* `burnin`: Number of burnin observations, default 500
* `pinfirst`: Set values for first observations (instead of burnin)

# Example
```julia-repl
model = CountModel(pastObs = 1)
simulate(100, model, [10, 0.5])
```
"""
function simulate(T::U where U<:Integer, model::T where T<:INGARCH, θ::parameter;
    burnin::U where U<:Integer = 500 ,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0) )

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    y = fill(0, T + burnin)
    X = model.X

    p = length(model.pastObs)
    q = length(model.pastMean)
    if p == 0
        P = 0
    else
        P = maximum(model.pastObs)
    end
    if q == 0
        Q = 0
    else
        Q = maximum(model.pastMean)
    end

    r = length(model.external)
    M = maximum([P, Q])

    if (size(X)[2] > 0) & (size(X)[2] != T)
        error("Number of regressors does not match T.")
    end

    logl = (model.link == "Log")
    lin = !logl

    nb = model.distr == "NegativeBinomial"
    zi = model.zi

    β0 = θ.β0
    α = θ.α
    β = θ.β
    η = θ.η
    ϕ = θ.ϕ
    ϕ = θ.ϕ
    ω = θ.ω

    rI = 0
    rE = 0
    if r > 0
        iE = findall(model.external)
        iI = setdiff(1:r, rE)
        rE = sum(model.external)
        rI = r - rE
    end

    ν = fill(β0, T + burnin)
    λ = fill(ifelse(lin, β0, exp(β0)), T + burnin)

    if p == 0
        if (q != 0) & (r == 0)
            error("INGARCH(0, q) without regressors not allowed.")
        end

        for i = 1:r
            ν[burnin + 1:end] .+= η[i]*X[i, :]
        end
        if logl
            λ = exp.(ν)
        end

        if nb
            y = rand.(NegativeBinomial.(ϕ[1], ϕ[1]./(ϕ[1] .+ λ)))
        else
            y = rand.(Poisson.(λ))
        end

        y = y.*(rand(T + burnin) .>= ω)
        λzi = λ.*(1 - ω)

        return y[burnin + 1:end], ν[burnin + 1:end], λ[burnin + 1:end], λzi[burnin + 1:end]
    end

    if p + q > 0
        νinit = β0/(1 - sum([α; β]))
    else
        vinit = β0
    end
    λinit = ifelse(logl, exp(νinit), νinit)

    ν[1:M] = fill(νinit, M)
    λ[1:M] = fill(λinit, M)

    if nb
        y[1:M] = rand.(NegativeBinomial.(ϕ[1], ϕ[1]./(ϕ[1] .+ λ[1:M])))
    else
        y[1:M] = rand.(Poisson.(λ[1:M]))
    end

    if zi
        y[1:M] = y[1:M].*(rand(M) .>= ω)
    end

    for t = (M+1):(T + burnin)
        ν[t] = β0
        if logl
            ν[t] += dot(α, log.(y[t .- model.pastObs] .+ 1))
        else
            ν[t] += dot(α, y[t .- model.pastObs])
        end

        if q > 0
            ν[t] += dot(β, ν[t .- model.pastMean])
        end

        if rI > 0
            ν[t] += dot(η[iI], X[iI, maximum([1, t - burnin])])
        end

        λ[t] = ν[t]

        if rE > 0
            λ[t] += dot(η[iE], X[iE, maximum([1, t - burnin])])
        end

        if logl
            λ[t] = exp(λ[t])
        end

        if nb
            y[t] = rand(NegativeBinomial(ϕ[1], ϕ[1]/(ϕ[1] + λ[t])))
        else
            y[t] = rand(Poisson(λ[t]))
        end

        if zi
            y[t] = y[t]*(rand() >= ω)
        end

        if length(pinfirst) > 0
            if burnin + 1 <= t <= burnin + length(pinfirst)
                y[t] = pinfirst[t - burnin]
                λ[t] = pinfirst[t - burnin]
                ν[t] = ifelse(logl, log(pinfirst[t - burnin] + 1), pinfirst[t - burnin])
            end
        end
    end

    λzi = λ.*(1 - ω)

    return y[burnin + 1:end], ν[burnin + 1:end], λ[burnin + 1:end], λzi[burnin+1:end]
end

function simulate(T::U where U<:Integer, model::INARCHModel, θ::parameter;
    burnin::U where U<:Integer = 500 ,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0) )

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    y = fill(0, T + burnin)
    X = model.X

    p = length(model.pastObs)
    if p == 0
        P = 0
    else
        P = maximum(model.pastObs)
    end

    r = length(model.external)
    M = P

    if (size(X)[2] > 0) & (size(X)[2] != T)
        error("Number of regressors does not match T.")
    end

    logl = (model.link == "Log")
    lin = !logl

    nb = model.distr == "NegativeBinomial"
    zi = model.zi

    β0 = θ.β0
    α = θ.α
    η = θ.η
    ϕ = θ.ϕ
    ϕ = θ.ϕ
    ω = θ.ω

    ν = fill(β0, T + burnin)
    λ = fill(ifelse(lin, β0, exp(β0)), T + burnin)

    if p == 0
        for i = 1:r
            ν[burnin + 1:end] .+= η[i]*X[i, :]
        end
        if logl
            λ = exp.(ν)
        end

        if nb
            y = rand.(NegativeBinomial.(ϕ[1], ϕ[1]./(ϕ[1] .+ λ)))
        else
            y = rand.(Poisson.(λ))
        end

        y = y.*(rand(T + burnin) .>= ω)
        λzi = λ.*(1 - ω)

        return y[burnin + 1:end], ν[burnin + 1:end], λ[burnin + 1:end], λzi[burnin + 1:end]
    end

    νinit = β0/(1 - sum(α))
    λinit = ifelse(logl, exp(νinit), νinit)

    ν[1:M] = fill(νinit, M)
    λ[1:M] = fill(λinit, M)

    if nb
        y[1:M] = rand.(NegativeBinomial.(ϕ[1], ϕ[1]./(ϕ[1] .+ λ[1:M])))
    else
        y[1:M] = rand.(Poisson.(λ[1:M]))
    end

    if zi
        y[1:M] = y[1:M].*(rand(M) .>= ω)
    end

    for t = (M+1):(T + burnin)
        ν[t] = β0
        if logl
            ν[t] += dot(α, log.(y[t .- model.pastObs] .+ 1))
        else
            ν[t] += dot(α, y[t .- model.pastObs])
        end

        if r > 0
            ν[t] += dot(η, X[:, maximum([1, t - burnin])])
        end

        λ[t] = ν[t]

        if logl
            λ[t] = exp(λ[t])
        end

        if nb
            y[t] = rand(NegativeBinomial(ϕ[1], ϕ[1]/(ϕ[1] + λ[t])))
        else
            y[t] = rand(Poisson(λ[t]))
        end

        if zi
            y[t] = y[t]*(rand() >= ω)
        end

        if length(pinfirst) > 0
            if burnin + 1 <= t <= burnin + length(pinfirst)
                y[t] = pinfirst[t - burnin]
                λ[t] = pinfirst[t - burnin]
                ν[t] = ifelse(logl, log(pinfirst[t - burnin] + 1), pinfirst[t - burnin])
            end
        end
    end

    λzi = λ.*(1 - ω)

    return y[burnin + 1:end], ν[burnin + 1:end], λ[burnin + 1:end], λzi[burnin+1:end]
end

function simulate(T::U where U<:Integer, model::IIDModel, θ::parameter;
    burnin::U where U<:Integer = 500 ,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0))

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    y = fill(0, T + burnin)
    X = model.X
    r = length(model.external)

    if (size(X)[2] > 0) & (size(X)[2] != T)
        error("Number of regressors does not match T.")
    end

    logl = (model.link == "Log")
    lin = !logl

    nb = model.distr == "NegativeBinomial"
    zi = model.zi

    β0 = θ.β0
    η = θ.η
    ϕ = θ.ϕ
    ϕ = θ.ϕ
    ω = θ.ω

    ν = fill(β0, T + burnin)
    λ = fill(ifelse(lin, β0, exp(β0)), T + burnin)

    for i = 1:r
        ν[burnin + 1:end] .+= η[i]*X[i, :]
    end
    if logl
        λ = exp.(ν)
    end

    if nb
        y = rand.(NegativeBinomial.(ϕ[1], ϕ[1]./(ϕ[1] .+ λ)))
    else
        y = rand.(Poisson.(λ))
    end

    y = y.*(rand(T + burnin) .>= ω)
    λzi = λ.*(1 - ω)

    return y[burnin + 1:end], ν[burnin + 1:end], λ[burnin + 1:end], λzi[burnin + 1:end]
end

function simulate(T::U where U<:Integer, model::T where T<:INARMA, θ::parameter;
    burnin::U where U<:Integer = 500 ,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0))

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    Y = fill(0, T + burnin)
    X = model.X

    p = length(model.pastObs)
    q = length(model.pastMean)
    if p == 0
        P = 0
    else
        P = maximum(model.pastObs)
    end
    if q == 0
        Q = 0
    else
        Q = maximum(model.pastMean)
    end

    r = length(model.external)

    if (size(X)[2] > 0) & (size(X)[2] != T)
        error("Number of regressors does not match T.")
    end

    M = maximum([P, Q])

    logl1 = (model.link[1] == "Log")
    logl2 = (model.link[2] == "Log")

    lin1 = !logl1
    lin2 = !logl2

    nb1 = model.distr[1] == "NegativeBinomial"
    nb2 = model.distr[2] == "NegativeBinomial"

    β0 = θ.β0
    α = θ.α
    β = θ.β
    η = θ.η
    ϕ = θ.ϕ
    ω = θ.ω

    rI = 0
    rE = 0
    iI = []
    iE = []
    if r > 0
        iE = findall(model.external)
        iI = setdiff(1:r, rE)
        rE = sum(iE)
        rI = r - rE
    end

    R = fill(0, T + burnin)
    Z = fill(0, T+ burnin)

    μ = zeros(T + burnin)
    λ = fill(β0, T + burnin)

    if rI > 0
        for i = 1:rI
            λ[1:burnin] .+= fill(η[iI[i]]*X[iI[i], 1], burnin)
            λ[burnin + 1:end] .+= η[iI[i]].*X[iI[i], :]
        end
    end

    if logl1
        λ = exp.(λ)
    end

    if rE > 0
        for i = 1:rE
            μ[1:burnin] .+= fill(η[iE[i]]*X[iE[i], 1], burnin)
            μ[burnin + 1:end] .+= η[iE[i]].*X[iE[i], :]
        end
        if logl2
            μ = exp.(μ)
        end
    end

    if nb1
        p1 = ϕ[1]./(ϕ[1] .+ λ)
        R = rand.(NegativeBinomial.(ϕ[1], p1))
    else
        R = rand.(Poisson.(λ))
    end

    if model.zi
        R = R.*(rand(length(R)) .> ω)
    end

    if nb2
        p2 = ϕ[1 + nb1]./(ϕ[1 + nb1] .+ μ)
        Z = rand.(NegativeBinomial.(ϕ[1 + nb1], p2))
    else
        Z = rand.(Poisson.(μ))
    end

    Y[1:M] = R[1:M]

    pf = length(pinfirst)
    pfBool = pf > 0

    for t = M+1:T+burnin
        Y[t] = R[t] + Z[t]
        if p > 0
            Y[t] += sum(α.∘Y[t .- model.pastObs])
        end

        if q > 0
            Y[t] += sum(β.∘R[t .- model.pastMean])
        end

        if pfBool
            if 1 <= (t - burnin) <= pf
                Y[t] = pinfirst[t - burnin]
            end
        end
    end

    return Y[burnin + 1:end], R[burnin+1:end]
end

function simulate(T::U where U<:Integer, model::INARModel, θ::parameter;
    burnin::U where U<:Integer = 500 ,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0))

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    Y = fill(0, T + burnin)
    X = model.X

    p = length(model.pastObs)
    if p == 0
        P = 0
    else
        P = maximum(model.pastObs)
    end

    r = length(model.external)

    if (size(X)[2] > 0) & (size(X)[2] != T)
        error("Number of regressors does not match T.")
    end

    M = P

    logl1 = (model.link[1] == "Log")
    logl2 = (model.link[2] == "Log")

    lin1 = !logl1
    lin2 = !logl2

    nb1 = model.distr[1] == "NegativeBinomial"
    nb2 = model.distr[2] == "NegativeBinomial"

    β0 = θ.β0
    α = θ.α
    η = θ.η
    ϕ = θ.ϕ
    ω = θ.ω

    rI = 0
    rE = 0
    iI = []
    iE = []
    if r > 0
        iE = findall(model.external)
        iI = setdiff(1:r, rE)
        rE = sum(iE)
        rI = r - rE
    end

    R = fill(0, T + burnin)
    Z = fill(0, T+ burnin)

    μ = zeros(T + burnin)
    λ = fill(β0, T + burnin)

    if rI > 0
        for i = 1:rI
            λ[1:burnin] .+= fill(η[iI[i]]*X[iI[i], 1], burnin)
            λ[burnin + 1:end] .+= η[iI[i]].*X[iI[i], :]
        end
    end

    if logl1
        λ = exp.(λ)
    end

    if rE > 0
        for i = 1:rE
            μ[1:burnin] .+= fill(η[iE[i]]*X[iE[i], 1], burnin)
            μ[burnin + 1:end] .+= η[iE[i]].*X[iE[i], :]
        end
        if logl2
            μ = exp.(μ)
        end
    end

    if nb1
        p1 = ϕ[1]./(ϕ[1] .+ λ)
        R = rand.(NegativeBinomial.(ϕ[1], p1))
    else
        R = rand.(Poisson.(λ))
    end

    if model.zi
        R = R.*(rand(length(R)) .> ω)
    end

    if nb2
        p2 = ϕ[1 + nb1]./(ϕ[1 + nb1] .+ μ)
        Z = rand.(NegativeBinomial.(ϕ[1 + nb1], p2))
    else
        Z = rand.(Poisson.(μ))
    end

    Y[1:M] = R[1:M]

    pf = length(pinfirst)
    pfBool = pf > 0

    for t = M+1:T+burnin
        Y[t] = R[t] + Z[t]
        if p > 0
            Y[t] += sum(α.∘Y[t .- model.pastObs])
        end

        if pfBool
            if 1 <= (t - burnin) <= pf
                Y[t] = pinfirst[t - burnin]
            end
        end
    end

    return Y[burnin + 1:end], R[burnin+1:end]
end

function simulate(T::U where U<:Integer, model::INMAModel, θ::parameter;
    burnin::U where U<:Integer = 500 ,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0))

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    Y = fill(0, T + burnin)
    X = model.X

    q = length(model.pastMean)
    if q == 0
        Q = 0
    else
        Q = maximum(model.pastMean)
    end

    r = length(model.external)

    if (size(X)[2] > 0) & (size(X)[2] != T)
        error("Number of regressors does not match T.")
    end

    M = Q

    logl1 = (model.link[1] == "Log")
    logl2 = (model.link[2] == "Log")

    lin1 = !logl1
    lin2 = !logl2

    nb1 = model.distr[1] == "NegativeBinomial"
    nb2 = model.distr[2] == "NegativeBinomial"

    β0 = θ.β0
    β = θ.β
    η = θ.η
    ϕ = θ.ϕ
    ω = θ.ω

    rI = 0
    rE = 0
    iI = []
    iE = []
    if r > 0
        iE = findall(model.external)
        iI = setdiff(1:r, rE)
        rE = sum(iE)
        rI = r - rE
    end

    R = fill(0, T + burnin)
    Z = fill(0, T+ burnin)

    μ = zeros(T + burnin)
    λ = fill(β0, T + burnin)

    if rI > 0
        for i = 1:rI
            λ[1:burnin] .+= fill(η[iI[i]]*X[iI[i], 1], burnin)
            λ[burnin + 1:end] .+= η[iI[i]].*X[iI[i], :]
        end
    end

    if logl1
        λ = exp.(λ)
    end

    if rE > 0
        for i = 1:rE
            μ[1:burnin] .+= fill(η[iE[i]]*X[iE[i], 1], burnin)
            μ[burnin + 1:end] .+= η[iE[i]].*X[iE[i], :]
        end
        if logl2
            μ = exp.(μ)
        end
    end

    if nb1
        p1 = ϕ[1]./(ϕ[1] .+ λ)
        R = rand.(NegativeBinomial.(ϕ[1], p1))
    else
        R = rand.(Poisson.(λ))
    end

    if model.zi
        R = R.*(rand(length(R)) .> ω)
    end

    if nb2
        p2 = ϕ[1 + nb1]./(ϕ[1 + nb1] .+ μ)
        Z = rand.(NegativeBinomial.(ϕ[1 + nb1], p2))
    else
        Z = rand.(Poisson.(μ))
    end

    Y[1:M] = R[1:M]

    pf = length(pinfirst)
    pfBool = pf > 0

    for t = M+1:T+burnin
        Y[t] = R[t] + Z[t]

        if q > 0
            Y[t] += sum(β.∘R[t .- model.pastMean])
        end

        if pfBool
            if 1 <= (t - burnin) <= pf
                Y[t] = pinfirst[t - burnin]
            end
        end
    end

    return Y[burnin + 1:end], R[burnin+1:end]
end

function simulate(T::U where U<:Integer, model::T where T<:CountModel, θ::Array{U, 1} where U<:AbstractFloat;
    burnin::U where U<:Integer = 500,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0))

    if !parametercheck(θ, model)
        error("Parameters not valid.")
    end

    simulate(T, model, θ2par(θ, model), burnin = burnin, pinfirst = pinfirst)
end

function simulate(T::U where U<:Integer, model::T where T<:CountModel, θ::Array{U, 1} where U<:Integer;
    burnin::U where U<:Integer = 500,
    pinfirst::Array{U, 1} where U<:AbstractFloat = Array{Float64, 1}(undef, 0))

    simulate(T, model, θ2par(θ, model), burnin = burnin, pinfirst = pinfirst)
end


text:
using BifurcationKitDocs
using Test

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

ext: jl


text:
#
# Metropolis-Hastings for logistic regression
# with missing input values
#

using PyPlot, PyCall
using Distributions
PyDict(matplotlib["rcParams"])["mathtext.fontset"] = "cm"
PyDict(matplotlib["rcParams"])["mathtext.rm"] = "serif"
PyDict(matplotlib["rcParams"])["lines.linewidth"] = 1.5
PyDict(matplotlib["rcParams"])["font.family"] = "TakaoPGothic"

function visualize_contour(samples_W, X, Y, text)
    R = 100
    xmin = minimum(X[1,:])
    xmax = maximum(X[1,:])
    ymin = minimum(X[2,:])
    ymax = maximum(X[2,:])
    lx = xmax - xmin
    ly = ymax - ymin
    xmin = xmin - 0.25 * lx
    xmax = xmax + 0.25 * lx
    ymin = ymin - 0.25 * ly
    ymax = ymax + 0.25 * ly

    x1 = linspace(xmin,xmax,R)
    x2 = linspace(ymin,ymax,R)
    x1grid = repmat(x1, 1, R)
    x2grid = repmat(x2', R, 1)
    val = [x1grid[:] x2grid[:]]'

    z_list = []
    for n in 1 : size(samples_W, 2)
        W = samples_W[:,n]
        z_tmp = [sigmoid(W'*val[:,i]) for i in 1 : size(val, 2)]
        push!(z_list, z_tmp)
    end
    z = mean(z_list)

    zgrid = reshape(z, R, R)

    # precition
    contour(x1grid, x2grid, zgrid, alpha=0.5, cmap=get_cmap("bwr"))
    scatter(X[1,Y.==1], X[2,Y.==1], c="r")
    scatter(X[1,Y.==0], X[2,Y.==0], c="b")
    xlim([xmin, xmax])
    ylim([ymin, ymax])
    title("$text")
end

function draw_line(W, xmin, xmax)
    y1 = - xmin*W[1]/W[2]
    y2 = - xmax*W[1]/W[2]
    plot([xmin, xmax], [y1, y2], c="k")
end

function sigmoid(x)
    return 1.0 / (1.0 + exp(-x[1]))
end

function bern_sample(mu)
    i = rand(Bernoulli(mu))
    val = zeros(2)
    val[i+1] = 1
    return val
end

function insert_X(X, X_miss)
    X_copy = deepcopy(X)
    X_copy[isnan(X)] = X_miss
    return X_copy
end

function MH(Y, X, M, Sigma_w, sigma_x, sigma_mh, max_iter, burnin)
    function ln_p_tilde(W, Y, X, Sigma_w, sigma_x)
        lkh = [Y[n]*log(sigmoid(W'*X[:, n])) + (1 - Y[n])*log(1 - sigmoid(W'*X[:, n])) for n in 1 : size(X, 2)]
        pri_w = -0.5 * (W'*inv(Sigma_w)*W)[1]
        pri_x = -0.5 * (sum(X.^2) * inv(sigma_x^2)) # Warning: Observation is not needed.
        return sum(lkh) + pri_w + pri_x
    end

    M = size(X, 1)
    N_miss = sum(isnan(X))

    ## sampling
    W = randn(M)
    X_miss = randn(N_miss)

    samples_W = Array{Float64, 2}(M, max_iter)
    samples_X = Array{Float64, 2}(N_miss, max_iter)
    for i in 1 : max_iter
        # sample W
        W_new = rand(MvNormal(W, sigma_mh^2*eye(M)))
        X_tmp = insert_X(X, X_miss)
        r = exp(minimum([ln_p_tilde(W_new, Y, X_tmp, Sigma_w, sigma_x) - ln_p_tilde(W, Y, X_tmp, Sigma_w, sigma_x), 0]))
        W = (rand() < r) ? W_new : W
        samples_W[:, i] = W

        # sample X_miss
        X_miss_new = rand(MvNormal(X_miss, sigma_mh^2*eye(N_miss)))
        X_tmp1 = insert_X(X, X_miss_new)
        X_tmp2 =  insert_X(X, X_miss)
        r = exp(minimum([ln_p_tilde(W, Y, X_tmp1, Sigma_w, sigma_x) - ln_p_tilde(W, Y, X_tmp2, Sigma_w, sigma_x), 0]))
        X_miss = (rand() < r) ? X_miss_new : X_miss
        samples_X[:, i] = X_miss
    end

    return samples_W[:,burnin+1:end], samples_X[:,burnin+1:end]
end

function test()
    ## compare results by discarding, interpolation and full observatino

    ########################
    # create model

    # hyperparameters
    M = 2
    Sigma_w = 20.0 * eye(M)

    ########################
    # create toy-data

    # sample parameters
    W = rand(MvNormal(zeros(M), Sigma_w))
    sigma_x = 1.0

    # sample data
    N = 30
    X_raw = 2 * rand(M, N) - 1.0
    Y = [rand(Bernoulli(sigmoid(W'*X_raw[:, n]))) for n in 1 : N]

    # missing values
    rate = 0.50
    X = deepcopy(X_raw)
    for n in 1 : N
        for m in 1 : M
            X[m,n] = (rand() < rate) ? NaN : X[m,n]
        end
    end

    idx = vec(!any(isnan(X), 1))
    X_red = X[:,idx]
    Y_red = Y[idx]

    ########################
    # inference
    sigma_mh = 0.10
    max_iter = 4000
    burnin = 1000

    samples_W1, samples_X1 = MH(Y_red, X_red, M, Sigma_w, sigma_x, sigma_mh, max_iter, burnin)
    samples_W2, samples_X2 = MH(Y, X, M, Sigma_w, sigma_x, sigma_mh, max_iter, burnin)
    samples_W3, samples_X3 = MH(Y, X_raw, M, Sigma_w, sigma_x, sigma_mh, max_iter, burnin)

    ########################
    # summarize
    idx_nan = vec(any(isnan(X), 1))
    X_est = insert_X(X, mean(samples_X2, 2))
    X_std = 0*deepcopy(X)
    X_std = insert_X(X_std, std(samples_X2, 2))

    ########################
    # visualize
    figure("result", figsize=(12,3))
    clf()
    subplot(1,3,1)
    visualize_contour(samples_W1, X_red, Y_red, "Reduced data (N=$(length(Y_red)))")

    subplot(1,3,2)
    Y_tmp = Y[idx_nan]
    X_tmp = X_est[:, idx_nan]
    scatter(X_tmp[1,Y_tmp.==1], X_tmp[2,Y_tmp.==1], c="r", marker="x")
    scatter(X_tmp[1,Y_tmp.==0], X_tmp[2,Y_tmp.==0], c="b", marker="x")
    for i in find(idx_nan)
        if Y[i] == 1
            plot([X_est[1,i] + X_std[1,i], X_est[1,i] - X_std[1,i]], [X_est[2,i], X_est[2,i]], "r", alpha=0.20, linewidth=5)
            plot([X_est[1,i], X_est[1,i]], [X_est[2,i] + X_std[2,i], X_est[2,i] - X_std[2,i]], "r", alpha=0.20, linewidth=5)
        else
            plot([X_est[1,i] + X_std[1,i], X_est[1,i] - X_std[1,i]], [X_est[2,i], X_est[2,i]], "b", alpha=0.20, linewidth=5)
            plot([X_est[1,i], X_est[1,i]], [X_est[2,i] + X_std[2,i], X_est[2,i] - X_std[2,i]], "b", alpha=0.20, linewidth=5)
        end
    end
    visualize_contour(samples_W2, X, Y, "Incomplete data (N=$N)")

    subplot(1,3,3)
    visualize_contour(samples_W3, X_raw, Y, "Full data (N=$N)")

    figure("traces", figsize=(12, 12))
    clf()
    subplot(4,1,1)
    plot(samples_W1')
    subplot(4,1,2)
    plot(samples_W2')
    subplot(4,1,3)
    plot(samples_X2')
    subplot(4,1,4)
    plot(samples_W3')

end

test()

max_stars_repo_name: sammy-suyama/mlblog

max_issues_repo_name: sammy-suyama/mlblog

max_forks_repo_name: sammy-suyama/mlblog


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using MocosSimLauncher

project_dir = joinpath(@__DIR__, "..")
cd(project_dir) do
  launch(["--help"])

  outdir = mktempdir()
  args = String[
    joinpath("example", "example.json"),
    "--output-summary", joinpath(outdir, "summary.jld2"),
    "--output-run-dump-prefix", joinpath(outdir, "run_dump"),
    "--output-daily", joinpath(outdir, "daily.jld2"),
    "--output-params-dump", joinpath(outdir, "params_dump.jld2"),
  ]
  launch(args)
end


text:
using Plots, PyPlot, FFTW
using Oceananigans

function rayleigh_benard_convection(Ra_desired, Nx, Ny, Nz, Lx, Ly, Lz, Nt, Δt)
    model = Model((Nx, Ny, Nz), (Lx, Ly, Lz))

    α = 207e-6  # Volumetric expansion coefficient [K⁻¹] of water at 20°C.
    ΔT = 1      # Temperature difference [K] between top and bottom.
    Pr = 0.7    # Prandtl number Pr = 𝜈/κ.

    # Calculate viscosity needed to get flow with desired Rayleigh number using
    # Eq. (3.5) of Kerr (1996).
    𝜈 = √((α*model.constants.g *Lz^3 / Ra_desired) * Pr * ΔT)
    Ra = (α*model.constants.g*Lz^3 / 𝜈^2) * Pr * ΔT

    println("Aspect ratio:    Γxz=$(Lx/Lz), Γyz=$(Ly/Lz)")
    println("Rayleigh number: Ra=$Ra")
    println("Prandtl number:  Pr=$Pr")

    # Create a new model based on the old with the configuration we want (we
    # just need to change 𝜈 and κ, the boundary conditions, and clock). This is
    # a hack until the Model constructor can be made more flexible.
    𝜈h, 𝜈v, κh, κv = 𝜈, 𝜈, 𝜈/Pr, 𝜈/Pr  # Assuming isotropic 𝜈 and κ.
    configuration = _ModelConfiguration(𝜈h, 𝜈v, κh, κv)
    boundary_conditions = BoundaryConditions(:periodic, :periodic, :rigid_lid, :no_slip)

    time, time_step = 0, 0
    clock = Clock(time, time_step, Δt)

    model = Model(model.metadata, configuration, boundary_conditions,
                  model.constants, model.eos, model.grid,
                  model.velocities, model.tracers, model.pressures,
                  model.G, model.Gp, model.forcings,
                  model.stepper_tmp, model.operator_tmp, model.poisson_solver, clock,
                  model.output_writers, model.diagnostics)

    # Write temperature field to disk every 10 time steps.
    output_dir, output_prefix, output_freq = ".", "rayleigh_benard", 10
    field_writer = FieldWriter(output_dir, output_prefix, output_freq, [model.tracers.T], ["T"])
    push!(model.output_writers, field_writer)

    diag_freq, Nu_running_avg = 1, 0
    Nu_wT_diag = Nusselt_wT(diag_freq, Float64[], Float64[], Nu_running_avg)
    push!(model.diagnostics, Nu_wT_diag)

    Nu_Chi_diag = Nusselt_Chi(diag_freq, Float64[], Float64[], Nu_running_avg)
    push!(model.diagnostics, Nu_Chi_diag)

    # Small random perturbations are added to boundary conditions to ensure instability formation.
    top_T    = 283 .- (ΔT/2) .+ 0.001.*rand(Nx, Ny)
    bottom_T = 283 .+ (ΔT/2) .+ 0.001.*rand(Nx, Ny)

    for i in 1:Nt
        time_step!(model; Nt=1, Δt=model.clock.Δt)
        # Impose constant T boundary conditions at top and bottom every time step.
        @. model.tracers.T.data[:, :,   1] = top_T
        @. model.tracers.T.data[:, :, end] = bottom_T
    end

    make_temperature_movie(model, field_writer)
    plot_Nusselt_number_diagnostics(model, Nu_wT_diag, Nu_Chi_diag)
end

function make_temperature_movie(model::Model, fw::FieldWriter)
    n_frames = Int(model.clock.time_step / fw.output_frequency)

    xC, yC, zC = model.grid.xC, model.grid.yC, model.grid.zC

    print("Creating temperature movie... ($n_frames frames)\n")

    Plots.gr()
    default(dpi=300)
    movie = @animate for tidx in 0:n_frames
        print("\rframe = $tidx / $n_frames   ")
        temperature = read_output(model, fw, "T", tidx*fw.output_frequency*model.clock.Δt)
        Plots.heatmap(xC, zC, rotl90(temperature[:, Int(ceil(model.grid.Ny/2)), :]) .- 283, color=:balance,
                      clims=(-0.5, 0.5), aspect_ratio=:equal,
                      title="T @ t=$(tidx*fw.output_frequency*model.clock.Δt)")
    end

    mp4(movie, "rayleigh_benard_$(round(Int, time())).mp4", fps = 30)
end

function plot_Nusselt_number_diagnostics(model::Model, Nu_wT_diag::Nusselt_wT, Nu_Chi_diag::Nusselt_Chi)
    println("Plotting Nusselt number diagnostics...")

    Nx, Ny, Nz = model.grid.Nx, model.grid.Ny, model.grid.Nz
    t = 0:model.clock.Δt:model.clock.time

    PyPlot.plot(t, Nu_wT_diag.Nu, label="Nu_wT")
    PyPlot.plot(t, Nu_wT_diag.Nu_inst, label="Nu_wT_inst")
    PyPlot.plot(t, Nu_Chi_diag.Nu, label="Nu_Chi")
    PyPlot.plot(t, Nu_Chi_diag.Nu_inst, label="Nu_Chi_inst")

    PyPlot.title("Rayleigh–Bénard convection ($Nx×$Ny×$Nz) @ Ra=5000")
    PyPlot.xlabel("Time (s)")
    PyPlot.ylabel("Nusselt number Nu")
    PyPlot.legend()
    PyPlot.savefig("rayleigh_benard_nusselt_diag.png", dpi=300, format="png", transparent=false)
end


text:
"""
    Mask{N, T}(classes) <: Block

Block for an N-dimensional categorical mask. `obs` is valid for
`Mask{N, T}(classes)`
if it is an N-dimensional array with every element in `classes`.
"""
struct Mask{N,T} <: Block
    classes::AbstractVector{T}
end
Mask{N}(classes::AbstractVector{T}) where {N,T} = Mask{N,T}(classes)

function checkblock(block::Mask{N,T}, a::AbstractArray{T,N}) where {N,T}
    return all(map(x -> x ∈ block.classes, a))
end

mockblock(mask::Mask{N, T}) where {N, T} = rand(mask.classes, ntuple(_ -> 16, N))::AbstractArray{T, N}

# Visualization

function showblock!(io, ::ShowText, block::Mask{2}, obs)
    img = _maskimage(obs, block.classes)
    ImageInTerminal.imshow(io, img)
end

function _maskimage(mask, classes)
    classtoidx = Dict(class => i for (i, class) in enumerate(classes))
    colors = distinguishable_colors(length(classes), transform = deuteranopic)
    return map(x -> colors[classtoidx[x]], mask)
end

_maskimage(mask::AbstractArray{<:Gray{T}}, args...) where {T} =
    _maskimage(reinterpret(T, mask), args...)
_maskimage(mask::AbstractArray{<:Normed{T}}, args...) where {T} =
    _maskimage(reinterpret(T, mask), args...)


text:
module LazyInverses
using LinearAlgebra
using Base.Threads: @spawn

const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}}
const parallel_threshold = 1024 # used to toggle parallel implementation of inverse Cholesky dot product (in algebra.jl)

export inverse, pinverse, pseudoinverse
export Inverse, PseudoInverse, AbstractInverse

# IDEA: extend Zygote's logdet adjoint with a lazy inverse ...
include("inverse.jl")
include("pseudo_inverse.jl")

const AbstractInverse{T, M} = Union{Inverse{T, M}, PseudoInverse{T, M}}

include("algebra.jl")
include("cholesky.jl") # contains specializations for inverses of Cholesky factorizations

end # LazyInverses


text:
# This file was generated, do not modify it. # hide
coerce!(df0, :Churn => OrderedFactor)
levels(df0.Churn) # to check order

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module STL

export parse, parse_malloc

struct Vertex
    x::Float32
    y::Float32
    z::Float32
end

const Normal = Vertex

struct Triangle
    normal::Normal
    v1::Vertex
    v2::Vertex
    v3::Vertex
end

function parse(path)
    open(path; lock = false) do io
    	skip(io, 80) # skip header
        triangle_count = read(io, UInt32)
    	triangles = Vector{Triangle}(undef, triangle_count) # preallocate memory for triangles
    	dest = Base.unsafe_convert(Ptr{Triangle}, triangles) # destination pointer
    	unsafe_read(io, dest, sizeof(Triangle)) # copying first triangle
        for _ in 2:triangle_count
            skip(io, 2)
            dest += sizeof(Triangle) # moving to the next trianlge in dest
            unsafe_read(io, dest, sizeof(Triangle))
        end
        triangles
    end
end

function parse_malloc(path)
    open(path; lock = false) do io
        skip(io, 80)
        triangle_count = read(io, UInt32)
        triangles = convert(Ptr{Triangle}, Base.Libc.malloc(sizeof(Triangle)*triangle_count))
        dest = triangles
        unsafe_read(io, dest, sizeof(Triangle)) # copying first triangle
        for _ in 2:triangle_count
            skip(io, 2)
            dest += sizeof(Triangle) # moving to the next trianlge in dest
            unsafe_read(io, dest, sizeof(Triangle))
        end
        Base.unsafe_wrap(Array, triangles, triangle_count; own = true)
    end
end

end #module STL