text:
"""
`FDGZProblem` - Initalize the data structure
"""
function FDGZProblem(ϵ, n)
    return FDGZProblem(0, 1, ϵ, n,  DirichletBC(-1), DirichletBC(1))
end

function FDGZProblem(a, b, ϵ, n, bcl::TBCL, bcr::TBCR) where {TBCL<:DirichletBC, TBCR<:DirichletBC}
    Δx = (b-a)/(n+1);
    x = LinRange(a, b, n+2)[2:end-1];
    F = zeros(n);

    return FDGZProblem(a, b, Δx, n, x, ϵ, F, bcl, bcr)
end


text:
# tests for the triggering module

@testset "Classic STA/LTA" begin
    # use random SeisChannel with fs > 1.0
    S = randSeisData(s=1.0, c=0.0, nx=10000, fs_min=1.0)
    sta = rand() .* 10.0
    lta = sta * 2.0
    U = classic_sta_lta(S, sta, lta)

    # should type, length, output
    @test isa(U, SeisData)
    @test [eltype(U.x[ii]) for ii in 1:U.n] == [eltype(S.x[ii]) for ii in 1:S.n]
    @test [length(S[ii].x) for ii in 1:S.n] == [length(U[ii].x) for ii in 1:S.n]
    @test U.x != S.x

    # test when sta > lta
    @test_throws DomainError classic_sta_lta(S,lta,sta)

    # test in-place
    classic_sta_lta!(S, sta, lta)
    @test U == S

    # test SeisChannel shorter than STA or LTA
    S = randSeisData(s=1.0,c=0.0,nx=100,fs_min=1.)
    @test_throws DomainError classic_sta_lta(S,-0.5,lta)
    @test_throws DomainError classic_sta_lta(S,-0.5,-0.25)

end

@testset "Recursive STA/LTA" begin
    # use random SeisChannel with fs > 1.0
    S = randSeisData(s=1.0, c=0.0, nx=10000, fs_min=1.0)
    sta = rand() .* 10.0
    lta = sta * 2.0
    U = recursive_sta_lta(S, sta, lta)

    # should type, length, output
    @test isa(U, SeisData)
    @test [eltype(U.x[ii]) for ii in 1:U.n] == [eltype(S.x[ii]) for ii in 1:S.n]
    @test [length(S[ii].x) for ii in 1:S.n] == [length(U[ii].x) for ii in 1:S.n]
    @test U.x != S.x

    # test when sta > lta
    @test_throws DomainError recursive_sta_lta(S,lta,sta)

    # test in-place
    recursive_sta_lta!(S, sta, lta)
    @test U == S

    # test SeisChannel shorter than STA or LTA
    S = randSeisData(s=1.0,c=0.0,nx=100,fs_min=1.)
    @test_throws DomainError recursive_sta_lta(S,-0.5,lta)
    @test_throws DomainError recursive_sta_lta(S,-0.5,-0.25)
end

@testset "Delayed STA/LTA" begin
        # use random SeisChannel with fs > 1.0
        S = randSeisData(s=1.0, c=0.0, nx=10000, fs_min=1.0)
        sta = rand() .* 10.0
        lta = sta * 2.0
        U = delayed_sta_lta(S, sta, lta)

        # should type, length, output
        @test isa(U, SeisData)
        @test [eltype(U.x[ii]) for ii in 1:U.n] == [eltype(S.x[ii]) for ii in 1:S.n]
        @test [length(S[ii].x) for ii in 1:S.n] == [length(U[ii].x) for ii in 1:S.n]
        @test U.x != S.x

        # test when sta > lta
        @test_throws DomainError delayed_sta_lta(S,lta,sta)

        # test in-place
        delayed_sta_lta!(S, sta, lta)
        @test U == S

        # test SeisChannel shorter than STA or LTA
        S = randSeisData(s=1.0,c=0.0,nx=100,fs_min=1.)
        @test_throws DomainError delayed_sta_lta(S,-0.5,lta)
        @test_throws DomainError delayed_sta_lta(S,-0.5,-0.25)
end

@testset "Carl STA/LTA" begin
        # use random SeisChannel with fs > 1.0
        S = randSeisData(s=1.0, c=0.0, nx=10000, fs_min=1.0)
        sta = rand() .* 10.0
        lta = sta * 2.0
        ratio = 0.8
        quiet = 0.8
        U = carl_sta_trig(S, sta, lta, ratio, quiet)

        # should type, length, output
        @test isa(U, SeisData)
        @test [eltype(U.x[ii]) for ii in 1:U.n] == [eltype(S.x[ii]) for ii in 1:S.n]
        @test [length(S[ii].x) for ii in 1:S.n] == [length(U[ii].x) for ii in 1:S.n]
        @test U.x != S.x

        # test when sta > lta
        @test_throws DomainError carl_sta_trig(S, lta, sta, ratio, quiet)

        # test in-place
        carl_sta_trig!(S, sta, lta, ratio, quiet)
        @test U == S

        # test SeisChannel shorter than STA or LTA
        S = randSeisData(s=1.0,c=0.0,nx=100,fs_min=1.)
        @test_throws DomainError carl_sta_trig(S,-0.5,lta,ratio,quiet)
        @test_throws DomainError carl_sta_trig(S,-0.5,-0.25,ratio,quiet)
end

@testset "Z-detect STA/LTA" begin
        # use random SeisChannel with fs > 1.0
        S = randSeisData(s=1.0, c=0.0, nx=10000, fs_min=1.0)
        sta = rand() .* 10.0
        U = z_detect(S, sta)

        # should type, length, output
        @test isa(U, SeisData)
        @test [eltype(U.x[ii]) for ii in 1:U.n] == [eltype(S.x[ii]) for ii in 1:S.n]
        @test [length(S[ii].x) for ii in 1:S.n] == [length(U[ii].x) for ii in 1:S.n]
        @test U.x != S.x

        # test in-place
        z_detect!(S, sta)
        @test U == S

        # test SeisChannel shorter than STA or LTA
        S = randSeisData(s=1.0,c=0.0,nx=100,fs_min=1.)
        @test_throws DomainError z_detect(S,-0.5)
end

@testset "Trigger on/off" begin
    # use random SeisChannel with fs > 1.0
    C = randSeisChannel(s=true,c=false,nx=10000,fs_min=1.)
    # make a fake earthquake
    C.x[4000:6000] .*= exp.(-log.(0.001:0.0005:1.001)) .+ 1
    C.x[7000:9000] .*= exp.(-log.(0.001:0.0005:1.001)) .+ 1

    # choose random window
    sta = rand() .* 10.0
    lta = sta * 5.0
    classic_sta_lta!(C,sta,lta)
    cft = C.x

    # on/off thresholds
    thresh1 = 4.0
    thresh2 = 2.0

    @test isa(trigger_onset(cft,thresh1,thresh2),Matrix{Int64})

    # test with incredibly high threshold
    thresh1 = Inf
    @test isa(trigger_onset(cft,thresh1,thresh2),Matrix)
    @test isempty(trigger_onset(cft,thresh1,thresh2))

    # test max event delete
    thresh1 = 3.5
    @test isempty(trigger_onset(cft,thresh1,thresh2,max_len=2,max_len_delete=true))
end

ext: jl


text:
# Automatically generated test file
function OEIS1()

g = [
1.000000000000000,
2.000000000000000,
4.000000000000000,
5.000000000000000,
9.000000000000000,
14.000000000000000,
8.000000000000000,
10.000000000000000,
18.000000000000000,
27.000000000000000,
17.000000000000000,
6.000000000000000,
11.000000000000000,
15.000000000000000,
29.000000000000000,
44.000000000000000,
19.000000000000000,
36.000000000000000,
54.000000000000000,
35.000000000000000,
21.000000000000000,
42.000000000000000,
23.000000000000000,
46.000000000000000,
25.000000000000000,
50.000000000000000,
28.000000000000000,
55.000000000000000,
83.000000000000000,
112.000000000000000,
31.000000000000000,
62.000000000000000,
33.000000000000000,
66.000000000000000,
37.000000000000000,
72.000000000000000,
108.000000000000000,
71.000000000000000,
39.000000000000000,
78.000000000000000,
41.000000000000000,
82.000000000000000,
124.000000000000000,
45.000000000000000,
89.000000000000000,
134.000000000000000,
88.000000000000000,
48.000000000000000,
96.000000000000000,
51.000000000000000,
101.000000000000000,
152.000000000000000,
53.000000000000000,
106.000000000000000,
160.000000000000000]

ζ_vec = [27558.9, 26259.0, 25020.4, 23840.2, 22715.7, 21644.3, 20623.3, 19650.6, 18723.7, 17840.5]

I_sols = [
[1, 55],
[1, 54, 55],
[1, 52, 53, 55],
[1, 51, 52, 53, 55],
[1, 46, 50, 52, 53, 55]]

f_sols = [
30839.557287157251267,
27558.879537513508694,
25479.067332593433093,
22300.676380090502789,
20535.515068432461703]

return g, ζ_vec, I_sols, f_sols

end


text:
export Residual

mutable struct Residual
    blocks::Vector
    function Residual(sequence::Vector)
        blocknum = length(sequence)
        blocks = Vector(undef,blocknum)
        for i = 1:blocknum
            blocks[i] = sequence[i]
        end
        new(blocks)
    end
    function Residual(sequence...)
        blocknum = length(sequence)
        blocks = Vector(undef,blocknum)
        for i = 1:blocknum
            blocks[i] = sequence[i]
        end
        new(blocks)
    end
end

function paramsof(m::Residual)
    params = Vector{Variable}(undef,0)
    for i = 1:length(m.blocks)
        append!(params, paramsof(m.blocks[i]))
    end
    return params
end

function xparamsof(m::Residual)
    params = Vector{XVariable}(undef,0)
    for i = 1:length(m.blocks)
        append!(params, xparamsof(m.blocks[i]))
    end
    return params
end

function forward(r::Residual, input::Variable)
    x = forward(r.blocks[1], input)
    for i = 2:length(r.blocks)
        x = forward(r.blocks[i], x)
    end
    @assert size(x)==size(input) "the input's dims should be the same as output's dims"
    return x + input
end

function predict(r::Residual, input)
    x = predict(r.blocks[1], input)
    for i = 2:length(r.blocks)
        x = predict(r.blocks[i], x)
    end
    return x + input
end

function nparamsof(r::Residual)
    c = 0
    for i = 1:length(r.blocks)
        c += nparamsof(r.blocks[i])
    end
    return c
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
using CSV, JLD2, FileIO, DataFrames

root_dir = dirname(pwd())
unzipped_files_path = *(root_dir, "/Datasets/Dataset_Bovespa/Unzipped/")
hdf5_file_name = (unzipped_files_path * "COTAHIST" * ".jld2")

csv = open((unzipped_files_path*"COTAHIST_2019"*".csv"),"w")
matrix = load(hdf5_file_name, "COTAHIST/2019")

df = DataFrame(data = [], codneg = [], preult = [])

for acao in eachrow(matrix)
    #if (acao[4] == 10 || acao[4] == 20) #mercado lote padrao ou fracionario
    if (acao[4] == 10)
        push!(df, [acao[1], rstrip(acao[3]), acao[13]])
    end
end

df = unstack(df, :data, :codneg, :preult) #pivota a tabela para as acoes
df = coalesce.(df, 0)  #substitui missing por 0
CSV.write(csv, df)


text:
#accept prob for γᵢ = ΠΠ(e(ηᵒy)+1)/(e(ηy)+1) #* exp(((γold - mu)²-(γnew-mu)²)/2σ²)exp((γᵒ-γⁿ)²/2gsd²)
function q_gamma(sample_gamma_new::Array{Float64,1}, sample_gamma::Array{Float64,1}, xb::Array{Float64,1}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, mu::Vector{Float64}, sigmas::Vector{Float64}, L::IntegerVector, L_new::IntegerVector, llvec::Vector{Float64}, llvecnew::Vector{Float64},ll_nF::Vector{Float64}, nF::Int, N::Int)

    relocate!(llvec, sample_gamma, groupindex, N)
    Yeppp.add!(llvec, llvec, xb)
    relocate!(llvecnew, sample_gamma_new, groupindex, N)
    Yeppp.add!(llvecnew, llvecnew, xb)

    loglogistic!(llvec, Y, N)
    loglogistic!(llvecnew, Y, N)
    negate!(llvec, llvec, N)
    Yeppp.add!(llvec, llvec, llvecnew)

    for i in 1:nF
        @inbounds ll_nF[i] = -(sample_gamma_new[i] - mu[L_new[i]])^2/sigmas[L_new[i]]^2*0.5 - log(sigmas[L_new[i]]) + (sample_gamma[i] - mu[L[i]])^2/sigmas[L[i]]^2*0.5 + log(sigmas[L[i]])
    end

    for i in 1:N
        @inbounds ll_nF[groupindex[i]] += llvec[i]
    end
    Yeppp.exp!(ll_nF, ll_nF)
    nothing
end

#The first part of Q function for beta, divided by -N
#the goal and gradient function for estimation β
#goal = 1/N ∑∑∑mean [log(1+exp(-η_{im}yᵢ)) for m in 1:M ]
#gradient = 1/N ∑∑∑mean [yᵢx[i,:]/(1+exp(η_{im}yᵢ)) for m in 1:M ]
function Q1(beta_new::Array{Float64,1}, storage::Vector, X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, sample_gamma_mat::Matrix{Float64},groupindex::IntegerVector, llvec::Vector{Float64}, llvecnew::Vector{Float64}, xb::Vector)

    N,J = size(X)
    # xb = X*beta_new
    A_mul_B!(xb, X, beta_new)
    ll=0.0
    if length(storage)>0
        fill!(storage, 0.0)
    end
    M = size(sample_gamma_mat,2)
    for jcol in 1:M
        #llvec[:] = xb
        relocate!(llvec, sample_gamma_mat[:,jcol], groupindex, N)
        Yeppp.add!(llvec, llvec, xb)
        #add!(llvec, sample_gamma_mat[groupindex, jcol], llvec)

        # map!(Exp_xy(), llvec, llvec, Y)
        negateiftrue!(llvec, Y)
        Yeppp.exp!(llvec, llvec)
        if length(storage) > 0
            llvecnew[:] = llvec
            # map!(Xy1x(), llvecnew, llvecnew, Y)
            x1x!(llvecnew)
            negateiffalse!(llvecnew, Y)

            for irow in 1:N
                for j in 1:J
                    @inbounds storage[j] += llvecnew[irow] * X[irow,j]
                end
            end
        end
        # map!(Log1p(), llvec, llvec)
        log1p!(llvec)
        ll += sum(llvec)
    end

    if length(storage)>0
        for j in 1:J
            storage[j] = storage[j] / M
        end
    end
    -ll/M
end

function gibbsMH!(X::Matrix, Y::Vector{Bool}, groupindex::IntegerVector, wi::Vector, mu::Vector, sigmas::Vector, β::Vector, ncomponent::Int, nF::Int, N::Int, M::Int, M_discard::Int, proposingsigma::RealVector, wipool::Vector, mupool::Vector, sigmaspool::Vector, L::Vector, L_new::Vector, sample_gamma::Vector, sample_gamma_new::Vector, sample_gamma_mat::Matrix, xb::Vector, llvec::Vector, llvecnew::Vector, ll_nF::Vector, tmp_mu::Vector, wi_divide_sigmas::Vector, inv_2sigmas_sq::Vector, tmp_p::Vector)
    fill!(wipool, 0.0)
    fill!(mupool, 0.0)
    fill!(sigmaspool, 0.0)
    #Gibbs samping for M+M_discard times
    # xb  = X*β
    A_mul_B!(xb, X, β)
    fill!(wi_divide_sigmas, 0.0)
    fill!(inv_2sigmas_sq, 0.0)
    for i in 1:length(wi)
        if sigmas[i] < realmin(Float64)
            wi_divide_sigmas[i] = 0.0
            inv_2sigmas_sq[i] = wi[i]*realmax(Float64)
        else
            wi_divide_sigmas[i] = wi[i]/sigmas[i]
            inv_2sigmas_sq[i] = 0.5 / sigmas[i]^2
        end
    end
    for iter_gibbs in 1:(M+M_discard)
        #update Lᵢ
        for i in 1:nF
            for j in 1:ncomponent
                tmp_mu[j] = -(mu[j] - sample_gamma[i])^2 * inv_2sigmas_sq[j]
            end
            ratiosumexp!(tmp_mu, wi_divide_sigmas, tmp_p, ncomponent)
            L_new[i] = rand(Categorical(tmp_p))
            sample_gamma_new[i] = rand(Normal(sample_gamma[i], proposingsigma[i]))
        end

        #update γᵢ;
        #Calculate the accept probability, stored in ll_nF
        q_gamma(sample_gamma_new, sample_gamma, xb, Y,groupindex, mu, sigmas, L, L_new, llvec, llvecnew, ll_nF, nF, N)
        for i in 1:nF
            if rand() < ll_nF[i]
                sample_gamma[i] = sample_gamma_new[i]
            end
            L[i] = L_new[i]
        end

        #only keep samples after M_discard
        jcol = iter_gibbs - M_discard
        if jcol > 0
            sample_gamma_mat[:, jcol] = sample_gamma
            addcounts!(wipool, L, 1:ncomponent)
            sumby!(mupool, sample_gamma, L, 1:ncomponent)
            sumsqby!(sigmaspool, sample_gamma, L, 1:ncomponent)
        end
    end

end
#The main function
#X, Y, groupindex
#nF is the number of groups
#intial values of β, ω, μ and σ must be supplied

function latentgmm(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, ncomponent::Int, β_init::Vector{Float64}, wi_init::Vector{Float64}, mu_init::Vector{Float64}, sigmas_init::Vector{Float64}; Mmax::Int=10000, M_discard::Int=1000, maxiteration::Int=100, initial_iteration::Int=0, tol::Real=.005, proposingsigma::RealVector=ones(maximum(groupindex)), ngh::Int=1000, sn::Vector{Float64}=sigmas_init, an::Float64=1.0/maximum(groupindex), debuginfo::Bool=false,restartMCMCsampling::Bool=false)

    # initialize theta
    length(wi_init) == length(mu_init) == length(sigmas_init) == ncomponent || error("The length of initial values should be $ncomponent")
    N,J=size(X)
    length(β_init) == J || error("Initial values of fixed efffect coefficients should have same dimension as X")
    nF = maximum(groupindex)
    M = Mmax

    wi = copy(wi_init)
    mu = copy(mu_init)
    sigmas = copy(sigmas_init)
    β = copy(β_init)

    wi_old = ones(wi)./ncomponent
    mu_old = zeros(mu)
    sigmas_old = ones(sigmas)
    beta_old = randn(J)

    ghx, ghw = gausshermite(ngh)
    #Preallocate the storage space, reusable for each iteration
    # L_mat = zeros(Int64, (nF, M))
    L = rand(Categorical(wi), nF)
    L_new = rand(Categorical(wi), nF)
    sample_gamma = zeros(nF)
    sample_gamma_new = zeros(nF)
    sample_gamma_mat = zeros(nF, Mmax)
    sumlogmat = zeros(nF, ngh*ncomponent)
    llvec = zeros(N)
    llvecnew = zeros(N)
    ll_nF = zeros(nF)

    wipool = zeros(ncomponent)
    mupool = zeros(ncomponent)
    sigmaspool = zeros(ncomponent)
    tmp_p=ones(ncomponent) / ncomponent
    tmp_mu=zeros(ncomponent)
    wi_divide_sigmas = zeros(wi)
    inv_2sigmas_sq = ones(sigmas) .* 1e20
    xb = zeros(N)
    ################################################################################################################
    #iterattion begins here

    no_iter=1
    M = min(2000, Mmax)
    Q_maxiter = 2
    for iter_em in 1:maxiteration
        if iter_em == (initial_iteration + 1)
            M = Mmax
            Q_maxiter = 10
        end
        if restartMCMCsampling || (iter_em == 1)
            for i in 1:nF
                L[i] = rand(Categorical(wi))
                sample_gamma[i] = rand(Normal(mu[L[i]], sigmas[L[i]]))
            end
         end
        if any(isnan(sigmas)) || any(isnan(wi))
            warn("wi=$wi, sigmas = $sigmas")
            return(wi, mu, sigmas, β, -Inf, [0.0])
        end

         gibbsMH!(X, Y, groupindex, wi, mu, sigmas, β, ncomponent, nF, N, M, M_discard, proposingsigma, wipool, mupool, sigmaspool, L, L_new, sample_gamma, sample_gamma_new, sample_gamma_mat, xb, llvec, llvecnew, ll_nF, tmp_mu, wi_divide_sigmas, inv_2sigmas_sq, tmp_p)

        for i in 1:nF
            proposingsigma[i] = std(sample_gamma_mat[i,:])+.1
        end
        for j in 1:ncomponent
            if wipool[j] == 0
                wipool[j] = 1
            end
        end
        copy!(wi_old, wi)
        copy!(mu_old, mu)
        copy!(sigmas_old, sigmas)
        #update wi, mu and sigmas
        wipoolsum = sum(wipool)
        for ic in 1:ncomponent
            wi[ic] = wipool[ic] / wipoolsum
            mu[ic] = mupool[ic] / wipool[ic]
            sigmas[ic] = sqrt((sigmaspool[ic] - wipool[ic] * mu[ic] ^2 + 2 * M * an * sn[ic]^2) / (wipool[ic] + 2 * M * an))
        end
        if debuginfo
            println("At $(iter_em)th iteration:")
        end
        #no longer update beta if it already converged
        if !stopRule(β, beta_old, tol=tol/10)
            copy!(beta_old, β)
            opt = Opt(:LD_LBFGS, J)
            maxeval!(opt, Q_maxiter)
            max_objective!(opt, (beta_new, storage)->Q1(beta_new, storage, X,Y, sample_gamma_mat[:,1:M], groupindex, llvec, llvecnew, xb))
            optimize!(opt, β)
            if debuginfo
                println("beta=", β)
            end
        end

        if debuginfo
            # println(β, "\t", wi, "\t", mu, "\t", sigmas, "\t", marginallikelihood(β, X, Y, groupindex, nF, wi, mu, sigmas, ghx, ghw, llvec, ll_nF, xb, sumlogmat)+sum(pn(sigmas, sn, an=an)))
            println("wi=$wi")
            println("mu=$mu")
            println("sigma=$sigmas")
            println("ll=",marginallikelihood(β, X, Y, groupindex, nF, wi, mu, sigmas, ghx, ghw, llvec, ll_nF, xb, sumlogmat))
        end
        if (iter_em == maxiteration) && (maxiteration > 3)
            warn("latentgmm not converge!")
        end
        if stopRule(vcat(β, wi, mu, sigmas), vcat(beta_old, wi_old, mu_old, sigmas_old), tol=tol) && (iter_em > initial_iteration) && (iter_em > 3)
            if debuginfo
                println("latentgmm converged at $(iter_em)th iteration")
            end
            break
        end
    end

     # xb=X*β
     # llmc = Float64[conditionallikelihood(xb, sample_gamma_mat[:,i], Y, groupindex) for i in 1:M]
    #m = MixtureModel(map((u, v) -> Normal(u, v), mu, sigmas), wi)

    return(wi, mu, sigmas, β, marginallikelihood(β, X, Y, groupindex, nF, wi, mu, sigmas, ghx, ghw, llvec, ll_nF, xb, sumlogmat)+sum(pn(sigmas, sn, an=an)), sample_gamma_mat[:,1:M])
end

#For fixed wi
function latentgmm_ctau(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, ncomponent::Int, β_init::Vector{Float64}, wi_init::Vector{Float64}, mu_init::Vector{Float64}, sigmas_init::Vector{Float64}, whichtosplit::Int64, tau::Float64, ghx::Vector{Float64}, ghw::Vector{Float64}; mu_lb::Vector{Float64}=-Inf.*ones(wi_init), mu_ub::Vector{Float64}=Inf.*ones(wi_init), Mmax::Int=5000, M_discard::Int=1000, maxiteration::Int=100, initial_iteration::Int=0, tol::Real=.005, proposingsigma::RealVector=ones(maximum(groupindex)), sn::Vector{Float64}=sigmas_init, an::Float64=0.25, debuginfo::Bool=false, sample_gamma_mat::Matrix = zeros(maximum(groupindex), Mmax), sumlogmat::Matrix = zeros(maximum(groupindex), length(ghx)*ncomponent), llvec::Vector=zeros(length(Y)), llvecnew::Vector = zeros(length(Y)), ll_nF::Vector = zeros(maximum(groupindex)), xb::Vector=zeros(length(Y)), Q_maxiter::Int = 5, restartMCMCsampling::Bool=false)

    # initialize theta
    N,J=size(X)
    nF = maximum(groupindex)
    M = Mmax
    tau = min(tau, 1-tau)

    wi = copy(wi_init)
    mu = copy(mu_init)
    sigmas = copy(sigmas_init)
    β = copy(β_init)

    wi_tmp = wi[whichtosplit]+wi[whichtosplit+1]
    wi[whichtosplit] = wi_tmp*tau
    wi[whichtosplit+1] = wi_tmp*(1-tau)
    mu = min(max(mu, mu_lb), mu_ub)

    wi_old = ones(wi)./ncomponent
    mu_old = zeros(mu)
    sigmas_old = ones(sigmas)
    beta_old = randn(J)

    wi0=copy(wi) # wi0, mu0 is to store the best parameter
    mu0=copy(mu)
    sigmas0=copy(sigmas)
    β0 = copy(β)

    #ghx, ghw = gausshermite(ngh)
    #Preallocate the storage space, reusable for each iteration
    # L_mat = zeros(Int64, (nF, M))
    L = rand(Categorical(wi), nF)
    L_new = rand(Categorical(wi), nF)
    sample_gamma = zeros(nF)
    sample_gamma_new = zeros(nF)
    fill!(sample_gamma_mat, 0.0)
    fill!(sumlogmat, 0.0)
    #llvec = zeros(N)
    #llvecnew = zeros(N)
    #ll_nF = zeros(nF)

    wipool = zeros(ncomponent)
    mupool = zeros(ncomponent)
    sigmaspool = zeros(ncomponent)

    tmp_p=ones(ncomponent) / ncomponent
    tmp_mu = zeros(ncomponent)
    wi_divide_sigmas = zeros(wi)
    inv_2sigmas_sq = ones(sigmas) .* 1e20
    ml0=-Inf

    ################################################################################################################
    #iterattion begins here

    no_iter=1
    #M = min(2000, Mmax)
    #Q_maxiter = 2
    lessthanmax = 0
    for iter_em in 1:maxiteration
        if restartMCMCsampling || (iter_em == 1)
            for i in 1:nF
                L[i] = rand(Categorical(wi))
                sample_gamma[i] = rand(Normal(mu[L[i]], sigmas[L[i]]))
            end
         end
         if any(isnan(sigmas)) || any(isnan(wi))
             warn("wi=$wi, sigmas = $sigmas")
             return(wi, mu, sigmas, β, -Inf)
         end
         gibbsMH!(X, Y, groupindex, wi, mu, sigmas, β, ncomponent, nF, N, M, M_discard, proposingsigma, wipool, mupool, sigmaspool, L, L_new, sample_gamma, sample_gamma_new, sample_gamma_mat, xb, llvec, llvecnew, ll_nF, tmp_mu, wi_divide_sigmas, inv_2sigmas_sq, tmp_p)
         for i in 1:nF
             proposingsigma[i] = std(sample_gamma_mat[i,:]) + .1
         end
         copy!(wi_old, wi)
         copy!(mu_old, mu)
         copy!(sigmas_old, sigmas)
         #update wi, mu and sigmas
         for j in 1:ncomponent
             if wipool[j] <= 1
                 warn("wi contains 0!")
                 return(wi, mu, sigmas, β, -Inf)
             end
         end
         wipoolsum = sum(wipool)
         for ic in 1:ncomponent
             wi[ic] = wipool[ic] / wipoolsum
             mu[ic] = mupool[ic] / wipool[ic]
             sigmas[ic] = sqrt((sigmaspool[ic] - wipool[ic] * mu[ic] ^2 + 2 * an * sn[ic].^2) / (wipool[ic] + 2 * an))
         end

        wi_tmp = wi[whichtosplit]+wi[whichtosplit+1]
        wi[whichtosplit] = wi_tmp*tau
        wi[whichtosplit+1] = wi_tmp*(1-tau)
        Yeppp.max!(mu, mu, mu_lb)
        Yeppp.min!(mu, mu, mu_ub)
        if debuginfo
            println("At $(iter_em)th iteration:")
        end
        #no longer update beta if it already converged
         if !stopRule(β, beta_old, tol=tol/10) #(mod(iter_em, 5) == 1 )
             copy!(beta_old, β)
             opt = Opt(:LD_LBFGS, J)
             maxeval!(opt, Q_maxiter)
             max_objective!(opt, (beta_new, storage)->Q1(beta_new, storage, X,Y, sample_gamma_mat[:, 1:M], groupindex, llvec, llvecnew, xb))
             optimize!(opt, β)
             if debuginfo
                 println("beta=", β)
             end
         end

        if (iter_em == maxiteration) && (maxiteration > 20)
            warn("latentgmm_ctau not yet converge!")
        end
        ml1 = marginallikelihood(β, X, Y, groupindex, nF, wi, mu, sigmas, ghx, ghw, llvec, ll_nF, xb, sumlogmat) + sum(pn(sigmas, sn, an=an))
        if debuginfo
            # println(wi, " ", mu, " ", sigmas, " ", sum(pn(sigmas, sn, an=an)), " ", ml1)
            println("wi=$wi")
            println("mu=$mu")
            println("sigma=$sigmas")
            println("ll=",marginallikelihood(β, X, Y, groupindex, nF, wi, mu, sigmas, ghx, ghw, llvec, ll_nF, xb, sumlogmat))
        end
        if ml1 > ml0
            ml0 = ml1
            copy!(wi0, wi)
            copy!(mu0, mu)
            copy!(sigmas0, sigmas)
            copy!(β0, β)
            lessthanmax = 0
        else
            lessthanmax += 1
        end
        if lessthanmax > 3
            if debuginfo
                println("latentgmm_ctau stop at $(iter_em)th iteration")
            end
            break
        end

    end
    #For fixed wi, no need to output gamma_mat
    return(wi0, mu0, sigmas0, β0, ml0)
end

#Starting from 25 initial values, find the best for fixed wi, used as start of the next 2 more iterations
function loglikelihoodratio_ctau(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, ncomponent1::Int,  betas0::Vector{Float64}, wi_C1::Vector{Float64},  whichtosplit::Int64, tau::Float64, mu_lb::Vector{Float64}, mu_ub::Vector{Float64}, sigmas_lb::Vector{Float64}, sigmas_ub::Vector{Float64}, gamma0::Vector{Float64}; ntrials::Int=25, ngh::Int=1000, sn::Vector{Float64}=sigmas_ub ./ 2, an=.25, debuginfo::Bool=false, sample_gamma_mat::Matrix = zeros(maximum(groupindex), Mmax), sumlogmat::Matrix = zeros(maximum(groupindex), length(ghx)*ncomponent), llvec::Vector=zeros(length(Y)), llvecnew::Vector = zeros(length(Y)), xb::Vector=zeros(length(Y)), Mctau::Int=1000, restartMCMCsampling::Bool=false)

    nF = maximum(groupindex)
    tau = min(tau, 1-tau)
    ghx, ghw = gausshermite(ngh)

    wi = repmat(wi_C1, 1, 4*ntrials)
    mu = zeros(ncomponent1, 4*ntrials)
    sigmas = ones(ncomponent1, 4*ntrials)
    betas = repmat(betas0, 1, 4*ntrials)
    ml = -Inf .* ones(4*ntrials)
    for i in 1:4*ntrials
        mu[:, i] = rand(ncomponent1) .* (mu_ub .- mu_lb) .+ mu_lb
        sigmas[:, i] = rand(ncomponent1) .* (sigmas_ub .- sigmas_lb) .+ sigmas_lb

        wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] = latentgmm_ctau(X, Y, groupindex, ncomponent1, betas0, wi[:, i], mu[:, i], sigmas[:, i], whichtosplit, tau, ghx, ghw, mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=5, Mmax=Mctau, M_discard=500, sn=sn, an=an, sample_gamma_mat = sample_gamma_mat, sumlogmat=sumlogmat, llvec=llvec, llvecnew=llvecnew, xb=xb, Q_maxiter=2, restartMCMCsampling=restartMCMCsampling)
    end

    mlperm = sortperm(ml)
    for j in 1:ntrials
        i = mlperm[4*ntrials+1 - j] # start from largest ml
        wi[:, i], mu[:, i], sigmas[:, i], betas[:, i], ml[i] = latentgmm_ctau(X, Y, groupindex, ncomponent1, betas[:, i], wi[:, i], mu[:, i], sigmas[:, i], whichtosplit, tau, ghx, ghw, mu_lb=mu_lb,mu_ub=mu_ub, maxiteration=200, initial_iteration=0, Mmax=Mctau, M_discard=500, sn=sn, an=an, debuginfo=debuginfo, sample_gamma_mat = sample_gamma_mat, sumlogmat=sumlogmat, llvec=llvec, llvecnew=llvecnew, xb=xb, Q_maxiter=10, restartMCMCsampling=restartMCMCsampling)
    end

    mlmax, imax = findmax(ml[mlperm[(3*ntrials+1):4*ntrials]])
    imax = mlperm[3*ntrials+imax]

    re=latentgmm(X, Y, groupindex, ncomponent1, betas[:, imax], wi[:, imax], mu[:, imax], sigmas[:, imax], Mmax=5000, maxiteration=3, initial_iteration=0, an=an, sn=sn, debuginfo=debuginfo, restartMCMCsampling=restartMCMCsampling)

    return(re[5])
end

function loglikelihoodratio(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, ncomponent1::Int; vtau::Vector{Float64}=[.5,.3,.1;], ntrials::Int=25, ngh::Int=1000, debuginfo::Bool=false, Mctau::Int=2000, restartMCMCsampling::Bool=false, ctauparallel=true)
    C0 = ncomponent1 - 1
    C1 = ncomponent1
    nF = maximum(groupindex)
    an1 = 1/nF
    gamma_init, betas_init, sigmas_tmp = maxposterior(X, Y, groupindex)
    wi_init, mu_init, sigmas_init, ml_tmp = gmm(gamma_init, C0, ones(C0)/C0, quantile(gamma_init, linspace(0, 1, C0+2)[2:end-1]), ones(C0), an=an1)

    wi_init, mu_init, sigmas_init, betas_init, ml_C0, gamma_mat = latentgmm(X, Y, groupindex, C0, betas_init, wi_init, mu_init, sigmas_init, Mmax=5000, initial_iteration=10, maxiteration=100, an=an1, sn=std(gamma_init).*ones(C0), restartMCMCsampling=restartMCMCsampling)

    if C0 > 1
        trand=LatentGaussianMixtureModel.asymptoticdistribution(X, Y, groupindex, wi_init, mu_init, sigmas_init, betas_init)
    end

    gamma0 = vec(mean(gamma_mat, 2))
    mingamma = minimum(gamma0)
    maxgamma = maximum(gamma0)

    or = sortperm(mu_init)
    wi0 = wi_init[or]
    mu0 = mu_init[or]
    sigmas0 = sigmas_init[or]
    betas0 = betas_init
    an = decidepenalty(wi0, mu0, sigmas0, nF)

    N,J=size(X)
    sample_gamma_mat = zeros(nF, Mctau)
    sumlogmat = zeros(nF, ngh*ncomponent1)
    llvec = zeros(N)
    llvecnew = zeros(N)
    xb = zeros(N)
    if ctauparallel
        lr=@parallel (max) for irun in 1:(C0*length(vtau))

            whichtosplit = mod1(irun, C0)
            i = cld(irun, C0)
            ind = [1:whichtosplit, whichtosplit:C0;]
            if C1==2
                mu_lb = mingamma .* ones(2)
                mu_ub = maxgamma .* ones(2)
            elseif C1>2
                mu_lb = [mingamma, (mu0[1:(C0-1)] .+ mu0[2:C0])./2;]
                mu_ub = [(mu0[1:(C0-1)] .+ mu0[2:C0])./2, maxgamma;]
                mu_lb = mu_lb[ind]
                mu_ub = mu_ub[ind]
            end
            sigmas_lb = 0.25 .* sigmas0[ind]
            sigmas_ub = 2 .* sigmas0[ind]

            wi_C1 = wi0[ind]
            wi_C1[whichtosplit] = wi_C1[whichtosplit]*vtau[i]
            wi_C1[whichtosplit+1] = wi_C1[whichtosplit+1]*(1-vtau[i])

            loglikelihoodratio_ctau(X, Y, groupindex, ncomponent1, betas0, wi_C1, whichtosplit, vtau[i], mu_lb, mu_ub,sigmas_lb, sigmas_ub, gamma0, ntrials=ntrials, ngh=ngh, sn=sigmas0[ind], an=an, debuginfo=debuginfo, sample_gamma_mat = sample_gamma_mat, sumlogmat=sumlogmat, llvec=llvec, llvecnew=llvecnew, Mctau=Mctau, xb=xb, restartMCMCsampling=restartMCMCsampling)

        end

        Tvalue = 2*(lr - ml_C0)
        if C0 == 1
            pvalue = 1 - cdf(Chisq(2), Tvalue)
        else
            pvalue = mean(trand .> 2*(lr - ml_C0))
        end
        return [Tvalue, pvalue;]
    else
        lr = zeros(length(vtau), C0)
        for whichtosplit in 1:C0, i in 1:length(vtau)

             #whichtosplit = mod1(irun, C0)
             #i = cld(irun, C0)
             ind = [1:whichtosplit, whichtosplit:C0;]
             if C1==2
                 mu_lb = mingamma .* ones(2)
                 mu_ub = maxgamma .* ones(2)
             elseif C1>2
                 mu_lb = [mingamma, (mu0[1:(C0-1)] .+ mu0[2:C0])./2;]
                 mu_ub = [(mu0[1:(C0-1)] .+ mu0[2:C0])./2, maxgamma;]
                 mu_lb = mu_lb[ind]
                 mu_ub = mu_ub[ind]
             end
             sigmas_lb = 0.25 .* sigmas0[ind]
             sigmas_ub = 2 .* sigmas0[ind]

             wi_C1 = wi0[ind]
             wi_C1[whichtosplit] = wi_C1[whichtosplit]*vtau[i]
             wi_C1[whichtosplit+1] = wi_C1[whichtosplit+1]*(1-vtau[i])

             lr[i, whichtosplit]=loglikelihoodratio_ctau(X, Y, groupindex, ncomponent1, betas0, wi_C1, whichtosplit, vtau[i], mu_lb, mu_ub,sigmas_lb, sigmas_ub, gamma0, ntrials=ntrials, ngh=ngh, sn=sigmas0[ind], an=an, debuginfo=debuginfo, sample_gamma_mat = sample_gamma_mat, sumlogmat=sumlogmat, llvec=llvec, llvecnew=llvecnew, Mctau=Mctau, xb=xb, restartMCMCsampling=restartMCMCsampling)
         end

         Tvalue = 2*(maximum(lr) - ml_C0)
         if C0 == 1
             pvalue = 1 - cdf(Chisq(2), Tvalue)
         else
             pvalue = mean(trand .> 2*(maximum(lr) - ml_C0))
         end
         return [Tvalue, pvalue;]
    end
end


text:
function assemble_bulk_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    cellquads,
    mesh,
)

    ncells = number_of_cells(mesh)
    dim = dimension(mesh)
    jac = jacobian(mesh)
    detjac = determinant_jacobian(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)
        if cellsign == +1 || cellsign == 0
            pquad = cellquads[+1, cellid]
            rhs = bulk_transformation_rhs(
                transfstress,
                basis,
                pquad,
                dim,
                jac,
                detjac,
                vectosymmconverter,
            )
            nodeids = nodal_connectivity(mesh, +1, cellid)
            assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
        end
    end
end

function assemble_interelement_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    facequads,
    mesh,
)

    dim = dimension(mesh)
    ncells = number_of_cells(mesh)
    nfaces = number_of_faces_per_cell(facequads)
    faceids = 1:nfaces
    nbrfaceids = [opposite_face(faceid) for faceid in faceids]
    normals = reference_face_normals()
    product_cellsign = +1
    facedetjac = face_determinant_jacobian(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)

        if cellsign == product_cellsign || cellsign == 0
            nodeids1 = nodal_connectivity(mesh, product_cellsign, cellid)

            for faceid in faceids
                nbrcellid = cell_connectivity(mesh, faceid, cellid)
                if cellid < nbrcellid &&
                   solution_cell_id(mesh, product_cellsign, cellid) !=
                   solution_cell_id(mesh, product_cellsign, nbrcellid)

                    nbrcellsign = cell_sign(mesh, nbrcellid)

                    if nbrcellsign == product_cellsign || nbrcellsign == 0
                        quad1 = facequads[product_cellsign, faceid, cellid]
                        quad2 = facequads[
                            product_cellsign,
                            nbrfaceids[faceid],
                            nbrcellid,
                        ]

                        nodeids2 = nodal_connectivity(
                            mesh,
                            product_cellsign,
                            nbrcellid,
                        )

                        assemble_face_interelement_transformation_linear_form!(
                            systemrhs,
                            -transfstress,
                            basis,
                            quad1,
                            normals[faceid],
                            nodeids1,
                            dim,
                            facedetjac[faceid],
                            vectosymmconverter,
                        )
                        assemble_face_interelement_transformation_linear_form!(
                            systemrhs,
                            transfstress,
                            basis,
                            quad2,
                            normals[faceid],
                            nodeids2,
                            dim,
                            facedetjac[nbrfaceids[faceid]],
                            vectosymmconverter,
                        )
                    end
                end
            end
        end
    end
end

function interelement_transformation_linear_form(
    transfstress,
    basis,
    quad,
    normal,
    dim,
    facedetjac,
    vectosymmconverter,
)

    numqp = length(quad)
    normals = repeat(normal, inner = (1, numqp))
    scalearea = repeat([facedetjac], numqp)
    return coherent_interface_transformation_linear_form(
        transfstress,
        basis,
        quad,
        normals,
        dim,
        scalearea,
        vectosymmconverter,
    )
end

function component_transformation_linear_form(
    transfstress,
    basis,
    quad,
    component,
    normal,
    dim,
    facedetjac,
    vectosymmconverter,
)

    numqp = length(quad)
    normals = repeat(normal, inner = (1, numqp))
    components = repeat(component, inner = (1, numqp))
    scalearea = repeat([facedetjac], numqp)
    return component_interface_transformation_linear_form(
        transfstress,
        basis,
        quad,
        components,
        normals,
        dim,
        scalearea,
        vectosymmconverter,
    )
end

function assemble_face_interelement_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    quad,
    normal,
    nodeids,
    dim,
    facedetjac,
    vectosymmconverter,
)

    rhs = interelement_transformation_linear_form(
        transfstress,
        basis,
        quad,
        normal,
        dim,
        facedetjac,
        vectosymmconverter,
    )
    assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
end

function assemble_face_component_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    quad,
    component,
    normal,
    nodeids,
    dim,
    facedetjac,
    vectosymmconverter,
)

    rhs = component_transformation_linear_form(
        transfstress,
        basis,
        quad,
        component,
        normal,
        dim,
        facedetjac,
        vectosymmconverter,
    )
    assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
end

function assemble_penalty_displacement_transformation_linear_form!(
    sysrhs,
    transfstress,
    basis,
    facequads,
    mesh,
    onboundary,
)

    facemidpoints = reference_face_midpoints()
    normals = reference_face_normals()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    dim = dimension(mesh)
    ncells = number_of_cells(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)

        if cellsign == +1 || cellsign == 0
            cellmap = cell_map(mesh, +1, cellid)
            for faceid = 1:nfaces
                if cell_connectivity(mesh, faceid, cellid) == 0 &&
                   onboundary(cellmap(facemidpoints[faceid]))

                    quad = facequads[+1, faceid, cellid]
                    nodeids = nodal_connectivity(mesh, +1, cellid)

                    assemble_face_interelement_transformation_linear_form!(
                        sysrhs,
                        -transfstress,
                        basis,
                        quad,
                        normals[faceid],
                        nodeids,
                        dim,
                        facedetjac[faceid],
                        vectosymmconverter,
                    )
                end
            end
        end
    end
end

function assemble_penalty_displacement_component_transformation_linear_form!(
    sysrhs,
    transfstress,
    basis,
    facequads,
    mesh,
    onboundary,
    component,
)

    facemidpoints = reference_face_midpoints()
    normals = reference_face_normals()
    nfaces = number_of_faces_per_cell(facequads)
    facedetjac = face_determinant_jacobian(mesh)
    dim = dimension(mesh)
    ncells = number_of_cells(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)

        if cellsign == +1 || cellsign == 0
            cellmap = cell_map(mesh, +1, cellid)
            for faceid = 1:nfaces
                if cell_connectivity(mesh, faceid, cellid) == 0 &&
                   onboundary(cellmap(facemidpoints[faceid]))

                    quad = facequads[+1, faceid, cellid]
                    nodeids = nodal_connectivity(mesh, +1, cellid)

                    assemble_face_component_transformation_linear_form!(
                        sysrhs,
                        -transfstress,
                        basis,
                        quad,
                        component,
                        normals[faceid],
                        nodeids,
                        dim,
                        facedetjac[faceid],
                        vectosymmconverter,
                    )
                end
            end
        end
    end
end

function assemble_coherent_interface_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    interfacequads,
    mesh,
)

    ncells = number_of_cells(mesh)
    dim = dimension(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh,cellid)

        if cellsign == 0
            normals = interface_normals(interfacequads, cellid)
            scalearea = interface_scale_areas(interfacequads, cellid)

            for s in [+1, -1]
                quad = interfacequads[s, cellid]
                rhs =
                    s *
                    0.5 *
                    coherent_interface_transformation_linear_form(
                        transfstress,
                        basis,
                        quad,
                        normals,
                        dim,
                        scalearea,
                        vectosymmconverter,
                    )
                nodeids = nodal_connectivity(mesh, s, cellid)
                assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
            end
        end
    end
end

function assemble_incoherent_interface_transformation_linear_form!(
    systemrhs,
    transfstress,
    basis,
    interfacequads,
    mesh,
)

    ncells = number_of_cells(mesh)
    dim = dimension(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh,cellid)

        if cellsign == 0
            normals = interface_normals(interfacequads, cellid)
            components = normals
            scalearea = interface_scale_areas(interfacequads, cellid)

            for s in [+1, -1]
                quad = interfacequads[s, cellid]
                rhs =
                    s *
                    0.5 *
                    component_interface_transformation_linear_form(
                        transfstress,
                        basis,
                        quad,
                        components,
                        normals,
                        dim,
                        scalearea,
                        vectosymmconverter,
                    )
                nodeids = nodal_connectivity(mesh, s, cellid)
                assemble_cell_rhs!(systemrhs, nodeids, dim, rhs)
            end
        end
    end
end


text:
module OpenPetris
import Base: ∈
using Petri
import SemanticModels.ModelTools.PetriModels: PetriModel
import SemanticModels.ModelTools.OpenModels: OpenModel
using ModelingToolkit
import ModelingToolkit: Constant
using Catlab.Doctrines
import Catlab.Doctrines: ⊗, compose, otimes

export OpenPetri, eye, ⊕

MAX_STATES = 20
X = @variables(X[1:MAX_STATES])[1]
STATELOOKUP = Dict(s.op.name=>i for (i,s) in enumerate(X))

const OpenPetri{V} = OpenModel{V, Petri.Model}

⊕(v::Vector, w::Vector) = vcat(v,w)
# ⊕(v::Vector{Int}, w::Vector{Int}) = vcat(v,w.+length(v))

function otimes(f::OpenPetri{T}, g::OpenPetri{T}) where {T<: Vector}
    f.model.S
    g.model.S
    # TODO: need to renumber the states of g
    M = Petri.Model(f.model.S ⊕ g.model.S, f.model.Δ ⊕ g.model.Δ)
    return OpenModel(f.dom ⊕ g.dom, M, f.codom ⊕ g.codom)
end

function otimes(f::OpenModel{T, Md}, g::OpenModel{T,Md}) where {T<: Vector{Int}, Md<:Petri.Model}
    f.model.S
    g.model.S
    nf = length(f.model.S)
    ng = length(g.model.S)
    newstates = Dict(X[s]=>X[s+nf] for (i, s) in enumerate(g.model.S))
    replace(t::Tuple{Operation, Operation}) = (replace(t[1]), replace(t[2]))
    replace(c::Constant) = c
    replace(op::Operation) = begin
        if op.op == (+)
            return sum(map(replace, op.args))
        end
        if op.op == (*)
            return prod(map(replace, op.args))
        end
        if length( op.args ) == 0
            return newstates[op]
        end
        return op
    end
    newtransitions = f.model.Δ
    if length(g.model.Δ) > 0
        newtransitions = newtransitions ⊕ map(g.model.Δ) do t
            replace(t)
        end
    end

    newstatespace = collect(1:(nf+ng))
    M = Petri.Model(newstatespace, newtransitions)
    return OpenModel(f.dom ⊕ (g.dom .+ nf), M, f.codom ⊕ (g.codom .+ nf))
end

function equnion(a::Vector, b::Vector)
    x = copy(a)
    for item in b
        if !any(item2 -> isequal(item2, item), x)
            push!(x, item)
        end
    end
    return x
end

∈(x::Operation, S::Vector{Operation}) = any(isequal.(x,S))

function compose(f::OpenModel{T, Md}, g::OpenModel{T, Md}) where {T<: Vector, Md<:Petri.Model}
    Y = f.codom
    Y′ = g.dom
    @assert length(Y) == length(Y′)
    Z = g.codom
    M = f.model
    N = g.model

    states = vcat(M.S, ( 1:length(filter(s->!(s ∈ Y′), N.S)) ) .+ length(M.S))
    newstates = Dict(X[Y′[i]]=>X[Y[i]] for i in 1:length(Y))
    i = 0
    newstates′ = map(N.S) do s
        if s ∈ Y′
            return nothing
        end
        i+=1
        X[s] => X[i+length(M.S)]
    end |> l-> filter(x-> x != nothing, l) |> Dict
    newstates = union(newstates, newstates′) |> Dict

    replace(t::Tuple{Operation, Operation}) = (replace(t[1]), replace(t[2]))
    replace(c::Constant) = c
    replace(op::Operation) = begin
        if op.op == (+)
            return sum(map(replace, op.args))
        end
        if op.op == (*)
            return prod(map(replace, op.args))
        end
        if length( op.args ) == 0
            # op ∈ keys(newstates), but for Operations
            if any(isequal.(keys(newstates), op))
                return newstates[op]
            end
            return op
        end
        return op
    end
    newtransitions = f.model.Δ
    if length(g.model.Δ) > 0
        newtransitions = newtransitions ⊕ map(g.model.Δ) do t
            replace(t)
        end
    end
    Δ = newtransitions
    Λ = vcat(M.Λ, N.Λ)
    Φ = vcat(M.Φ, N.Φ)
    Mp_yN = Petri.Model(states, Δ, Λ, Φ)
    Z′ = map(Z) do z
        findfirst(x->isequal(x, newstates[X[z]]), X)
    end
    return OpenModel(f.dom, Mp_yN, Z′)
end

eye(n::Int) = foldr(otimes, [OpenModel([1], NullPetri(1), [1]) for i in 1:n])

end


text:
module PLite

include("common_defs.jl")
include("common_checks.jl")
include("solvers.jl")

end


text:
# # [Quantum Circuit Born Machine](@id qcbm)

# Yao is designed with variational quantum circuits in mind, and this tutorial
# will introduce how to use Yao for this kind of task by implementing a quantum
# circuit born machine described in [Jin-Guo Liu, Lei Wang (2018)](https://arxiv.org/abs/1804.04168)

# let's use the packages first

using Yao, LinearAlgebra, Plots

# ## Training Target

# In this tutorial, we will ask the variational circuit to learn the most basic
# distribution: a guassian distribution. It is defined as follows:

# ```math
# f(x \left| \mu, \sigma^2\right) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
# ```

# We implement it as `gaussian_pdf`:

function gaussian_pdf(x, μ::Real, σ::Real)
    pl = @. 1 / sqrt(2pi * σ^2) * exp(-(x - μ)^2 / (2 * σ^2))
    pl / sum(pl)
end
pg = gaussian_pdf(1:1<<6, 1<<5-0.5, 1<<4);

# We can plot the distribution, it looks like

plot(pg)

# ## Create the Circuit

# A quantum circuit born machine looks like the following:

# ![differentiable ciruit](assets/differentiable.png)

# It is composited by two different layers: **rotation layer** and **entangler layer**.

# ## Rotation Layer

# Arbitrary rotation is built with **Rotation Gate** on **Z, X, Z** axis
# with parameters.

# ```math
# Rz(\theta) \cdot Rx(\theta) \cdot Rz(\theta)
# ```

# Since our input will be a ``|0\dots 0\rangle`` state.
# The first layer of arbitrary rotation can just
# use ``Rx(\theta) \cdot Rz(\theta)`` and the last
# layer of arbitrary rotation could just
# use ``Rz(\theta)\cdot Rx(\theta)``

# In **幺**, every Hilbert operator is a **block** type, this
# ncludes all **quantum gates** and **quantum oracles**.
# In general, operators appears in a quantum circuit
# can be divided into **Composite Blocks** and **Primitive Blocks**.

# We follow the low abstraction principle and
# thus each block represents a certain approach
# of calculation. The simplest **Composite Block**
# is a **Chain Block**, which chains other blocks
# (oracles) with the same number of qubits together.
# It is just a simple mathematical composition of
# operators with same size. e.g.

# ```math
# \text{chain(X, Y, Z)} \iff X \cdot Y \cdot Z
# ```

# We can construct an arbitrary rotation block by chain ``Rz``, ``Rx``, ``Rz`` together.

chain(Rz(0.0), Rx(0.0), Rz(0.0))

# `Rx`, `Rz` will construct new rotation gate,
# which are just shorthands for `rot(X, 0.0)`, etc.

# Then let's chain them up

layer(nbit::Int, x::Symbol) = layer(nbit, Val(x))
layer(nbit::Int, ::Val{:first}) = chain(nbit, put(i=>chain(Rx(0), Rz(0))) for i = 1:nbit);

# We do not need to feed the first `n` parameter into `put` here.
# All factory methods can be **lazy** evaluate **the first arguements**, which is the number of qubits.
# It will return a lambda function that requires a single interger input.
# The instance of desired block will only be constructed until all the information is filled.
# When you filled all the information in somewhere of the declaration, 幺 will be able to infer the others.
# We will now define the rest of rotation layers

layer(nbit::Int, ::Val{:last}) = chain(nbit, put(i=>chain(Rz(0), Rx(0))) for i = 1:nbit)
layer(nbit::Int, ::Val{:mid}) = chain(nbit, put(i=>chain(Rz(0), Rx(0), Rz(0))) for i = 1:nbit);

# ## Entangler

# Another component of quantum circuit born machine are
# several **CNOT** operators applied on different qubits.

entangler(pairs) = chain(control(ctrl, target=>X) for (ctrl, target) in pairs);

# We can then define such a born machine

function build_circuit(n, nlayers, pairs)
    circuit = chain(n)
    push!(circuit, layer(n, :first))
    for i in 2:nlayers
        push!(circuit, cache(entangler(pairs)))
        push!(circuit, layer(n, :mid))
    end
    push!(circuit, cache(entangler(pairs)))
    push!(circuit, layer(n, :last))
    return circuit
end

# We use the method `cache` here to tag the entangler block that it
# should be cached after its first run, because it is actually a
# constant oracle. Let's see what will be constructed

build_circuit(4, 1, [1=>2, 2=>3, 3=>4])

# ## MMD Loss & Gradients

# The MMD loss is describe below:

# ```math
# \begin{aligned}
# \mathcal{L} &= \left| \sum_{x} p \theta(x) \phi(x) - \sum_{x} \pi(x) \phi(x) \right|^2\\
#             &= \langle K(x, y) \rangle_{x \sim p_{\theta}, y\sim p_{\theta}} - 2 \langle K(x, y) \rangle_{x\sim p_{\theta}, y\sim \pi} + \langle K(x, y) \rangle_{x\sim\pi, y\sim\pi}
# \end{aligned}
# ```

# We will use a squared exponential kernel here.

struct RBFKernel
    σ::Float64
    m::Matrix{Float64}
end

function RBFKernel(σ::Float64, space)
    dx2 = (space .- space').^2
    return RBFKernel(σ, exp.(-1/2σ * dx2))
end

kexpect(κ::RBFKernel, x, y) = x' * κ.m * y

# There are two different way to define the loss:

# In simulation we can use the probability distribution of the state directly

get_prob(qcbm) = probs(zero_state(nqubits(qcbm)) |> qcbm)

function loss(κ, c, target)
    p = get_prob(c) - target
    return kexpect(κ, p, p)
end

# Or if you want to simulate the whole process with measurement (which is entirely
# physical), you should define the loss with measurement results, for convenience
# we directly use the simulated results as our loss

# ### Gradients
# the gradient of MMD loss is

# ```math
# \begin{aligned}
# \frac{\partial \mathcal{L}}{\partial \theta^i_l} &= \langle K(x, y) \rangle_{x\sim p_{\theta^+}, y\sim p_{\theta}} - \langle K(x, y) \rangle_{x\sim p_{\theta}^-, y\sim p_{\theta}}\\
# &- \langle K(x, y) \rangle _{x\sim p_{\theta^+}, y\sim\pi} + \langle K(x, y) \rangle_{x\sim p_{\theta^-}, y\sim\pi}
# \end{aligned}
# ```

# which can be implemented as

function gradient(qcbm, κ, ptrain)
    n = nqubits(qcbm)
    prob = get_prob(qcbm)
    grad = zeros(Float64, nparameters(qcbm))

    count = 1
    for k in 1:2:length(qcbm), each_line in qcbm[k], gate in content(each_line)
        dispatch!(+, gate, π/2)
        prob_pos = probs(zero_state(n) |> qcbm)

        dispatch!(-, gate, π)
        prob_neg = probs(zero_state(n) |> qcbm)

        dispatch!(+, gate, π/2) # set back

        grad_pos = kexpect(κ, prob, prob_pos) - kexpect(κ, prob, prob_neg)
        grad_neg = kexpect(κ, ptrain, prob_pos) - kexpect(κ, ptrain, prob_neg)
        grad[count] = grad_pos - grad_neg
        count += 1
    end
    return grad
end

# Now let's setup the training

using Flux.Optimise
qcbm = build_circuit(6, 10, [1=>2, 3=>4, 5=>6, 2=>3, 4=>5, 6=>1])
dispatch!(qcbm, :random) # initialize the parameters

κ = RBFKernel(0.25, 0:2^6-1)
pg = gaussian_pdf(1:1<<6, 1<<5-0.5, 1<<4);
opt = ADAM()

function train(qcbm, κ, opt, target)
    history = Float64[]
    for _ in 1:100
        push!(history, loss(κ, qcbm, target))
        ps = parameters(qcbm)
        Optimise.update!(opt, ps, gradient(qcbm, κ, target))
        popdispatch!(qcbm, ps)
    end
    return history
end

history = train(qcbm, κ, opt, pg)
trained_pg = probs(zero_state(nqubits(qcbm)) |> qcbm)

# The history of training looks like below

title!("training history")
xlabel!("steps"); ylabel!("loss")
plot(history)

# And let's check what we got

fig2 = plot(1:1<<6, trained_pg; label="trained")
plot!(fig2, 1:1<<6, pg; label="target")
title!("distribution")
xlabel!("x"); ylabel!("p")

# So within 50 steps, we got a pretty close estimation of our
# target distribution!


text:
module TestKMPP

using ParallelKMeans
using Test
using Random

# resulting indices obtained from `Clustering` implementation of kmpp
@testset "singlethread kmpp" begin
    Random.seed!(1980)
    design_matrix = rand(10, 1000)

    Random.seed!(2020)
    @test ParallelKMeans.smart_init(design_matrix, 10, 1).indices == [33, 931, 853, 940, 926, 528, 644, 552, 460, 433]
end

@testset "multithread kmpp" begin
    Random.seed!(1980)
    design_matrix = rand(10, 1000)

    Random.seed!(2020)
    @test ParallelKMeans.smart_init(design_matrix, 10, 2).indices == [33, 931, 853, 940, 926, 528, 644, 552, 460, 433]
end

end # module


text:
create_udp_entry(packageid::Integer) = create_package(rdb_pkg_id_to_type(packageid))
function create_udp_entry(::Type{RDB_SYNC_t},
    mask::UInt32=0x00000000,
    cmdMask::UInt32=RDB_SYNC_CMD_RENDER_SINGLE_FRAME, #RDB_SYNC_CMD_RENDER_CONTINUOUS,
    systemTime::Cdouble=time()
    )

    RDB_SYNC_t(mask, cmdMask, systemTime)
end
function create_udp_entry(::Type{RDB_TRIGGER_t},
    deltaT::Cfloat,
    frameNo::UInt32,
    spare::Int32=convert(Int32, 0)
    )

    RDB_TRIGGER_t(deltaT, frameNo, spare)
end

function create_udp_entry_header{T<:RDB_PACKAGE_ELEMENT}(elements::Vector{T}, flags::UInt16=0x0000)
    headerSize = sizeof(RDB_MSG_ENTRY_HDR_t)
    #dataSize = sizeof(elements)
    #elementSize = div(dataSize, length(elements))
    dataSize = sizeof(elements[1])*length(elements)
    elementSize = sizeof(elements[1])
    #pkgId = rdb_type_to_pkg_id(T)
    pkgId = rdb_type_to_pkg_id(elements[1])
    RDB_MSG_ENTRY_HDR_t(headerSize, dataSize, elementSize, pkgId, flags)
end

function create_udp_entry(::Type{RDB_DRIVER_CTRL_t},
    playerId::UInt32;
    steeringWheel::Cfloat=convert(Cfloat,0),
    steeringSpeed::Cfloat=convert(Cfloat,0),
    throttlePedal::Cfloat=convert(Cfloat,0),
    brakePedal::Cfloat=convert(Cfloat,0),
    clutchPedal::Cfloat=convert(Cfloat,0),
    accelTgt::Cfloat=convert(Cfloat,0),
    steeringTgt::Cfloat=convert(Cfloat,0),
    curvatureTgt::Cdouble=convert(Cdouble,0),
    steeringTorque::Cfloat=convert(Cfloat,0),
    engineTorqueTgt::Cfloat=convert(Cfloat,0),
    speedTgt::Cfloat=convert(Cfloat,0),
    gear::UInt8=convert(UInt8,0),
    sourceId::UInt8=convert(UInt8,0),
    spare0::Tuple{UInt8, UInt8}=(convert(UInt8,0),convert(UInt8,0)),
    validityFlags::UInt32=convert(UInt32,0),
    flags::UInt32=convert(UInt32,0),
    mockupInput0::UInt32=convert(UInt32,0),
    mockupInput1::UInt32=convert(UInt32,0),
    mockupInput2::UInt32=convert(UInt32,0),
    spare::UInt32=convert(UInt32,0)
    )

    RDB_DRIVER_CTRL_t(    playerId,
    steeringWheel,
    steeringSpeed,
    throttlePedal,
    brakePedal,
    clutchPedal,
    accelTgt,
    steeringTgt,
    curvatureTgt,
    steeringTorque,
    engineTorqueTgt,
    speedTgt,
    gear,
    sourceId,
    spare0,
    validityFlags,
    flags,
    mockupInput0,
    mockupInput1,
    mockupInput2,
    spare)

end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
using SpecialFunctions, ForwardDiff, FiniteDifferences
include("shared.jl")

const BIG_FD = central_fdm(10,1)

dwbesselk(v,x) = ForwardDiff.derivative(_v->BesselK._besselk(_v,x,100,1e-12,false), v)

dbesselk(v, x) = BIG_FD(_v->besselk(_v,x), v)

const _h = 1e-6
fastfdbesselk(v, x) = (besselk(v+_h, x) - besselk(v, x))/_h

const BASELINE = [dbesselk(z[1], z[2])  for z in Iterators.product(VGRID, XGRID)]
const OURSOL   = [dwbesselk(z[1], z[2]) for z in Iterators.product(VGRID, XGRID)]
const FASTFD   = [fastfdbesselk(z[1], z[2])  for z in Iterators.product(VGRID, XGRID)]
const TOLS     = atolfun.(zip(BASELINE, OURSOL))
const FDTOLS   = atolfun.(zip(BASELINE, FASTFD))
const DIFTOLS  = log10.(TOLS) .- log10.(FDTOLS)

let res = findmax(DIFTOLS)
  println("Worst case derivative difference: $(res[1])")
  println("Value of dfun: $(BASELINE[res[2]])")
end

gnuplot_save_matrix!("../plotdata/atols_deriv_fd.csv",    FDTOLS, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_deriv_ad.csv",      TOLS, VGRID, XGRID)
gnuplot_save_matrix!("../plotdata/atols_deriv_fdad.csv", DIFTOLS, VGRID, XGRID)


text:
module TimeDataFrames

using CSV
using DataFrames
using ExtendedDates

import Base: eachcol, eachrow, show

export TimeDataFrame, periods, firstperiod, lastperiod, dataframe, innerjoin, outerjoin, lag, lead, ncol, nrow, index

mutable struct TimeDataFrame
    data::DataFrame
    periods::AbstractVector{ExtendedDates.SimpleDate}
    continuous::Bool
    function TimeDataFrame(data::DataFrame,
                           periods::AbstractVector{T},
                           continuous::Bool) where T <: ExtendedDates.SimpleDate
        new(data, periods, continuous)
    end
end;

TimeDataFrame() = TimeDataFrame(DataFrame(), Vector{T}(), true) where T <: ExtendedDates.SimpleDate

function TimeDataFrame(dataframe::AbstractDataFrame, firstperiod::T;
                       copycols::Bool=true) where T <: ExtendedDates.SimpleDate
    periods = range(firstperiod, length=DataFrames.nrow(dataframe), step = typeof(firstperiod - firstperiod)(1))
    TimeDataFrame(DataFrame(dataframe; copycols), periods, true)
end

function TimeDataFrame(filename::String, firstperiod::T) where T <: ExtendedDates.SimpleDate
    data = DataFrame(CSV.File(filename))
    continuous = true
    periods = range(firstperiod, length=size(data, 1), step = typeof(firstperiod - firstperiod)(1))
    TimeDataFrame(data, periods, true)
end

TimeDataFrame(firstperiod::T, pairs::Pair{Symbol,<:Any}...;
              makeunique::Bool=false, copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(pairs; makeunique, copycols),
                                firstperiod)

TimeDataFrame(firstperiod::T, pairs::Pair{<:AbstractString,<:Any}...;
              makeunique::Bool=false, copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(pairs; makeunique, copycols),
                                firstperiod)

# these two are needed as a workaround Tables.jl dispatch
TimeDataFrame(pairs::AbstractVector{<:Pair}, firstperiod::T; makeunique::Bool=false,
              copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(pairs..., makeunique=makeunique, copycols=copycols),
                                firstperiod)

TimeDataFrame(pairs::NTuple{N, Pair}, firstperiod::T; makeunique::Bool=false,
              copycols::Bool=true) where {N,T <: ExtendedDates.SimpleDate} =
                  TimeDataFrame(DataFrame(pairs..., makeunique=makeunique, copycols=copycols),
                                firstperiod)

TimeDataFrame(d::AbstractDict, firstperiod::T; copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
    TimeDataFrame(DataFrame(d; copycols), firstperiod)

#=
TimeDataFrame(frequency::ExtendedDates.Frequency, firstperiod::ExtendedDates.SimpleDate; kwargs...)
    if isempty(kwargs)
        TimeDataFrame()
    else
        cnames = Symbol[]
        columns = Any[]
        copycols = true
        for (kw, val) in kwargs
            if kw == :copycols
                if val isa Bool
                    copycols = val
                else
                    throw(ArgumentError("the `copycols` keyword argument must be Boolean"))
                end
            else
                push!(cnames, kw)
                push!(columns, val)
            end
        end
        TimeDataFrame(DataFrame(columns, Index(cnames), copycols=copycols), Undated, 1)
    end
end
=#

TimeDataFrame(columns::AbstractVector, cnames::AbstractVector{Symbol},
              firstperiod::T; makeunique::Bool=false, copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(columns, cnames; makeunique, copycols), firstperiod)

TimeDataFrame(columns::AbstractVector, cnames::AbstractVector{<:AbstractString},
              firstperiod::T; makeunique::Bool=false, copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(columns, Symbol.(cnames), makeunique=makeunique, copycols=copycols),
                                firstperiod)

TimeDataFrame(columns::AbstractVector{<:AbstractVector}, firstperiod::T,
              cnames::AbstractVector{Symbol}=gennames(length(columns));
              makeunique::Bool=false, copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(collect(AbstractVector, columns),
                                          Index(convert(Vector{Symbol}, cnames), makeunique=makeunique),
                                          copycols=copycols), firstperiod)

TimeDataFrame(columns::AbstractVector{<:AbstractVector},
              cnames::AbstractVector{<:AbstractString}, firstperiod::T;
              makeunique::Bool=false, copycols::Bool=true) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(columns, Symbol.(cnames); makeunique=makeunique, copycols=copycols),
                                firstperiod)

TimeDataFrame(columns::NTuple{N, AbstractVector}, cnames::NTuple{N, Symbol}, firstperiod::T;
              makeunique::Bool=false, copycols::Bool=true) where {N, T <: ExtendedDates.SimpleDate} =
                  TimeDataFrame(DataFrame(collect(AbstractVector, columns), collect(Symbol, cnames),
                                          makeunique=makeunique, copycols=copycols), firstperiod)

TimeDataFrame(columns::NTuple{N, AbstractVector}, cnames::NTuple{N, AbstractString}, firstperiod::T;
          makeunique::Bool=false, copycols::Bool=true) where {N, T <: ExtendedDates.SimpleDate} =
              TimeDataFrame(DataFrame(columns, Symbol.(cnames); makeunique=makeunique, copycols=copycols),
                            firstperiod)

TimeDataFrame(columns::NTuple{N, AbstractVector}, firstperiod::T; copycols::Bool=true) where {N, T <: ExtendedDates.SimpleDate} =
    TimeDataFrame(DataFrame(collect(AbstractVector, columns), gennames(length(columns)),
                            copycols=copycols), firstperiod)

TimeDataFrame(columns::AbstractMatrix, firstperiod::T,
          cnames::AbstractVector{Symbol} = gennames(size(columns, 2)); makeunique::Bool=false) where T <: ExtendedDates.SimpleDate =
              TimeDataFrame(DataFrame(AbstractVector[columns[:, i] for i in 1:size(columns, 2)], cnames,
                                      makeunique=makeunique, copycols=false), firstperiod)

TimeDataFrame(columns::AbstractMatrix, cnames::AbstractVector{<:AbstractString}, firstperiod::T;
              makeunique::Bool=false) where T <: ExtendedDates.SimpleDate =
                  TimeDataFrame(DataFrame(columns, Symbol.(cnames); makeunique=makeunique), firstperiod)

TimeDataFrame(column_eltypes::AbstractVector{T}, cnames::AbstractVector{Symbol},
              firstperiod::P, nrows::Integer=0; makeunique::Bool=false) where {T<:Type, P <: ExtendedDates.SimpleDate} =
                  TimeDataFrame(DataFrame(column_eltypes, cnames, nrows; makeunique), firstperiod)

TimeDataFrame(column_eltypes::AbstractVector{T}, cnames::AbstractVector{S}, firstperiod::P,
              nrows::Integer=0; makeunique::Bool=false) where {T <: Type, P <: ExtendedDates.SimpleDate, S <: AbstractString} =
                  TimeDataFrame(DataFrame(column_eltypes, Symbol.(cnames), nrows; makeunique=makeunique),
                                firstperiod)

# Redefining getproperty breaks tdf.data !
function Base.getproperty(tdf::TimeDataFrame, symbol::Symbol)
    data = getfield(tdf, :data)
    x = getproperty(data, symbol)
    periods = getfield(tdf, :periods)
    continuous = getfield(tdf, :continuous)
    df = DataFrame()
    df[!, symbol] = x
    TimeDataFrame(df,  periods, continuous)
end

Base.names(tdf::TimeDataFrame) = DataFrames.names(DataFrames.index(getfield(tdf, :data)))
_names(tdf::TimeDataFrame) = DataFrames._names(DataFrames.index(getfield(tdf, :data)))

function Base.setproperty!(tdf::TimeDataFrame, symbol::Symbol, x::AbstractVector)
    data = getfield(tdf, :data)
    periods = getfield(tdf, :periods)
    data[!, symbol] = x
    continuous = getfield(tdf, :continuous)
    TimeDataFrame(data, periods, continuous)
end

function Base.setproperty!(tdf1::TimeDataFrame, col_ind::Symbol, tdf2::TimeDataFrame)
    n1 = nrow(tdf1)
    n2 = nrow(tdf2)
    continuous2 = getfield(tdf2, :continuous)
    data1 = getfield(tdf1, :data)
    data2 = getfield(tdf2, :data)
    periods2 = getfield(tdf2, :periods)
    if ncol(tdf1) == 0
        insertcols!(data1, 1, col_ind => missings(n2))
    else
        continuous1 = getfield(tdf1, :continuous)
        periods1 = getfield(tdf1, :periods)
        if typeof(periods1) != typeof(periods2)
            error("The frequency must be the same in both TimeDataFrame")
        end
        minperiod = (periods1[1] < periods2[1]) ? periods1[1] : periods2[1]
        maxperiod = (periods1[n1] > periods2[n2]) ? periods1[n1] : periods2[n2]
        # Missing expanding number of periods in data1 if necessary
    end
#    if !(col_ind in DataFrames.names(data1))
#    end
    setindex!(data1, data2[!, 1], ! , col_ind)
#    rename!(data1, [col_ind])
    setfield!(tdf1, :data, data1)
    setfield!(tdf1, :periods, periods2)
    setfield!(tdf1, :continuous, continuous2)
end

Base.getindex(df::TimeDataFrame, id1::Integer, id2::Integer) = getindex(getfield(df, :data), id1, id2)
Base.getindex(df::TimeDataFrame, idx::CartesianIndex{2}) = df[idx[1], idx[2]]
Base.getindex(df::TimeDataFrame, ::typeof(!), col_ind::Symbol) = getindex(getfield(df, :data), !, col_ind)
Base.getindex(df::TimeDataFrame, ::typeof(!), col_ind::Union{Signed, Unsigned}) = getindex(getfield(df, :data), !, col_ind)
Base.view(df::TimeDataFrame, idx::CartesianIndex{2}) = view(df, idx[1], idx[2])
Base.setindex!(tdf::TimeDataFrame, val, idx::CartesianIndex{2}) =
    (tdf[idx[1], idx[2]] = val)
Base.setindex!(tdf::TimeDataFrame, x::AbstractArray{Float64, 1}, ::typeof(!), col_ind::Symbol) =
               setindex!(getfield(tdf, :data), x, !, col_ind)
Base.setindex!(tdf::TimeDataFrame, x::AbstractArray{Union{Missing, Float64}, 1}, ::typeof(!), col_ind::Symbol) =
               setindex!(getfield(tdf, :data), x, !, col_ind)
Base.setindex!(tdf::TimeDataFrame, x::AbstractArray{Float64, 1}, ::typeof(!), col_ind::Union{Signed, Unsigned}) =
    setindex!(getfield(tdf, :data), x, !, col_ind)

Base.broadcastable(tdf::TimeDataFrame) = tdf
struct TimeDataFrameStyle <: Base.Broadcast.BroadcastStyle end

Base.Broadcast.BroadcastStyle(::Type{<:TimeDataFrame}) =
    TimeDataFrameStyle()

Base.Broadcast.BroadcastStyle(::TimeDataFrameStyle, ::Base.Broadcast.BroadcastStyle) = TimeDataFrameStyle()
Base.Broadcast.BroadcastStyle(::Base.Broadcast.BroadcastStyle, ::TimeDataFrameStyle) = TimeDataFrameStyle()
Base.Broadcast.BroadcastStyle(::TimeDataFrameStyle, ::TimeDataFrameStyle) = TimeDataFrameStyle()

function copyto_widen!(res::AbstractVector{T}, bc::Base.Broadcast.Broadcasted,
                       pos, col) where T
    for i in pos:length(axes(bc)[1])
        val = bc[CartesianIndex(i, col)]
        S = typeof(val)
        if S <: T || promote_type(S, T) <: T
            res[i] = val
        else
            newres = similar(Vector{promote_type(S, T)}, length(res))
            copyto!(newres, 1, res, 1, i-1)
            newres[i] = val
            return copyto_widen!(newres, bc, i + 1, col)
        end
    end
    return res
end

function getcolbc(bcf::Base.Broadcast.Broadcasted{Style}, colind) where {Style}
    # we assume that bcf is already flattened and unaliased
    newargs = map(bcf.args) do x
        Base.Broadcast.extrude(x isa TimeDataFrame ? x[!, colind] : x)
    end
    Base.Broadcast.Broadcasted{Style}(bcf.f, newargs, bcf.axes)
end

function Base.copy(bc::Base.Broadcast.Broadcasted{TimeDataFrameStyle})
    ndim = length(axes(bc))
    if ndim != 2
        throw(DimensionMismatch("cannot broadcast a time data frame into $ndim dimensions"))
    end

    bcf = Base.Broadcast.flatten(bc)
    first_tdf = true
    local periods, continuous
    colnames = []
    for tdf in bcf.args
        if tdf isa TimeDataFrame
            if first_tdf
                periods = getfield(tdf, :periods)
                continuous = getfield(tdf, :continuous)
                first_tdf = false
            elseif getfield(tdf, :periods) != periods
                error("Time data frames don't have the same periods")
            elseif getfield(tdf, :continuous) != continuous
                error("TimeDataFrames don't have the same continuous status")
            end
        end
    end
    colnames = unique!([_names(df) for df in bcf.args if df isa TimeDataFrame])
    if any(length(c) != 1 for c in colnames)
        wrongnames = setdiff(union(colnames...), intersect(colnames...))
        if isempty(wrongnames)
            throw(ArgumentError("Column names in broadcasted time data frames " *
                                "must have the same order"))
        else
            msg = join(wrongnames, ", ", " and ")
            throw(ArgumentError("Column names in broadcasted time data frames must match. " *
                                "Non matching column names are $msg"))
        end
    end
    nrows = length(axes(bcf)[1])
    tdf = TimeDataFrame(DataFrame(), periods, continuous)
    for i in axes(bcf)[2]
        if nrows == 0
            col = Any[]
        else
            bcf′ = getcolbc(bcf, i)
            v1 = bcf′[CartesianIndex(1, i)]
            startcol = similar(Vector{typeof(v1)}, nrows)
            startcol[1] = v1
            col = copyto_widen!(startcol, bcf′, 2, i)
        end
        tdf[!, colnames[1][i]] = col
    end
    return tdf
end

function Base.copyto!(tdf::TimeDataFrame, bc::Base.Broadcast.Broadcasted)
    TimeDataFrame(copyto!(dataframe(tdf), bc),
                  periods(tdf),
                  iscontinuous(tdf))
end

Base.ndims(::TimeDataFrame) = 2
Base.ndims(::Type{<:TimeDataFrame}) = 2
index(df::TimeDataFrame) = getfield(getfield(df, :data), :colindex)
_columns(df::TimeDataFrame) = getfield(getfield(df, :data), :columns)

import Base.size
Base.size(df::TimeDataFrame) = (nrow(df), ncol(df))
function Base.size(df::TimeDataFrame, i::Integer)
    if i == 1
        nrow(df)
    elseif i == 2
        ncol(df)
    else
        throw(ArgumentError("DataFrames only have two dimensions"))
    end
end

function innerjoin(d1::TimeDataFrame, d2::TimeDataFrame)
    continuous1 = getfield(d1, :continuous)
    continuous2 = getfield(d2, :continuous)
    data1 = getfield(d1, :data)
    data2 = getfield(d2, :data)
    periods1 = getfield(d1, :periods)
    periods2 = getfield(d2, :periods)
    if typeof(periods1) != typeof(periods2)
        error("innerjoin: both TimeDataFrames must have the same frequency")
    end
    return TimeDataFrame(sort!(DataFrames.innerjoin(data1, data2, on=:Column1),1), intersect(periods1, periods2), true)
end

function outerjoin(d1::TimeDataFrame, d2::TimeDataFrame)
    continuous1 = getfield(d1, :continuous)
    continuous2 = getfield(d2, :continuous)
    data1 = getfield(d1, :data)
    data2 = getfield(d2, :data)
    periods1 = getfield(d1, :periods)
    periods2 = getfield(d2, :periods)
    if typeof(periods1) != typeof(periods2)
        error("innerjoin: both TimeDataFrames must have the same frequency")
    end
    return TimeDataFrame(sort!(DataFrames.outerjoin(data1, data2, on=:Column1),1), union(periods1, periods2), true)
end

function Base.show(io::IO, tdf::TimeDataFrame)
    df = getfield(tdf, :data)
    dfcopy = copy(df)
    periods = getfield(tdf, :periods)
    insertcols!(dfcopy, 1, :Periods => periods)
    show(io, dfcopy, show_row_number = false, eltypes = false, summary = false)
end

function Base.show(tdf::TimeDataFrame)
    df = getfield(tdf, :data)
    dfcopy = copy(df)
    periods = getfield(tdf, :periods)
    insertcols!(dfcopy, 1, :Periods => periods)
    show(dfcopy, show_row_number = false, eltypes = false, summary = false)
end

function Base.isequal(tdf1::TimeDataFrame, tdf2::TimeDataFrame)
    isequal(getfield(tdf1, :data), getfield(tdf2, :data)) || return false
    isequal(getfield(tdf1, :periods), getfield(tdf2, :periods)) || return false
    getfield(tdf1, :continuous) == getfield(tdf2, :continuous) || return false
    return true
end

"""
    copy(tdf::TimeDataFrame; copycols::Bool=true)

Copy time data frame `tdf`.
If `copycols=true` (the default), return a new  `TimeDataFrame` holding
copies of column vectors in `tdf`.
If `copycols=false`, return a new `TimeDataFrame` sharing column vectors with `tdf`.
"""
function Base.copy(tdf::TimeDataFrame; copycols::Bool=true)
    return TimeDataFrame(copy(dataframe(tdf), copycols=copycols), copy(periods(tdf)), iscontinuous(tdf))
end
##############################################################################
##
## Equality
##
##############################################################################

function Base.:(==)(tdf1::TimeDataFrame, tdf2::TimeDataFrame)
    return isequal(tdf1, tdf2)
end

"""
    isapprox(tdf1::TimeDataFrame, tdf2::TimeDataFrame;
             rtol::Real=atol>0 ? 0 : √eps, atol::Real=0,
             nans::Bool=false, norm::Function=norm)

Inexact equality comparison. `tdf1` and `tdf2` must have the same size, column names and periods.
Return  `true` if `isapprox` with given keyword arguments
applied to all pairs of columns stored in `df1` and `df2` returns `true`.
"""
function Base.isapprox(tdf1::TimeDataFrame, tdf2::TimeDataFrame;
                       atol::Real=0, rtol::Real=atol>0 ? 0 : √eps(),
                       nans::Bool=false, norm::Function=norm)
    if size(tdf1) != size(tdf2)
        throw(DimensionMismatch("dimensions must match: a has dims " *
                                "$(size(tdf1)), b has dims $(size(tdf2))"))
    end
    if !isequal(index(tdf1), index(tdf2))
        throw(ArgumentError("column names of passed time data frames do not match"))
    end
    if !isequal(periods(tdf1), periods(tdf2))
        throw(ArgumentError("periods  of passed time data frames do not match"))
    end
    return isapprox(dataframe(tdf1), dataframe(tdf2), atol=atol, rtol=rtol, nans=nans, norm=norm)
end

Base.Matrix(tdf::TimeDataFrame) = Matrix(dataframe(tdf))
Base.Matrix{T}(tdf::TimeDataFrame) where T = Matrix{T}(dataframe(tdf))
Base.Array(tdf::TimeDataFrame) = Array(dataframe(tdf))
Base.Array{T}(tdf::TimeDataFrame) where T = Array{T}(dataframe(tdf))

include("accessors.jl")
export iscontinuous, dataframe, firstperiod, lastperiod, periods
include("dataframe_functions.jl")
include("timeseries_functions.jl")
export lag, lead, align!
include("plot.jl")

end # module


text:
using AbstractDifferentiation
using Distributions
using ForwardDiff
using GalacticOptim
using LinearAlgebra
using Optim
using Pathfinder
using Random
using ReverseDiff
using SciMLBase
using Test
using Transducers

@testset "single path pathfinder" begin
    @testset "IsoNormal" begin
        # here pathfinder finds the exact solution after 1 iteration
        logp(x) = -sum(abs2, x) / 2
        ∇logp(x) = -x
        ndraws = 100
        rngs = if VERSION ≥ v"1.7"
            [MersenneTwister(), Random.default_rng()]
        else
            [MersenneTwister()]
        end
        seed = 42
        @testset for dim in [1, 5, 10, 100], rng in rngs
            executor = rng isa MersenneTwister ? SequentialEx() : ThreadedEx()

            init = randn(dim)
            Random.seed!(rng, seed)
            result = @inferred pathfinder(logp, ∇logp; init, ndraws, rng, executor)
            @test result isa PathfinderResult
            @test result.input === (logp, ∇logp)
            @test result.optim_prob isa SciMLBase.OptimizationProblem
            @test result.logp(init) ≈ logp(init)
            @test result.rng === rng
            @test result.optimizer ===
                Pathfinder.default_optimizer(Pathfinder.DEFAULT_HISTORY_LENGTH)
            fit_distribution = result.fit_distribution
            @test fit_distribution isa MvNormal
            @test fit_distribution.μ ≈ zeros(dim)
            @test fit_distribution.Σ isa Pathfinder.WoodburyPDMat
            @test fit_distribution.Σ ≈ I
            @test size(fit_distribution.Σ.B) == (dim, 2) # history contains only 1 iteration
            @test result.draws isa AbstractMatrix
            @test size(result.draws) == (dim, ndraws)
            @test result.fit_distribution_transformed === result.fit_distribution
            @test result.draws_transformed === result.draws
            @test result.num_tries ≥ 1
            @test result.optim_solution isa SciMLBase.OptimizationSolution
            @test result.optim_trace isa Pathfinder.OptimizationTrace
            @test result.fit_distributions isa Vector{typeof(fit_distribution)}
            @test length(result.fit_distributions) == length(result.optim_trace)
            @test result.fit_distributions[result.fit_iteration + 1] == fit_distribution
            @test result.fit_iteration ==
                argmax(getproperty.(result.elbo_estimates, :value))

            Random.seed!(rng, seed)
            result2 = pathfinder(logp, ∇logp; init, ndraws, rng, executor)
            @test result2.fit_iteration == result.fit_iteration
            @test result2.draws == result.draws
            @test getproperty.(result2.elbo_estimates, :value) ==
                getproperty.(result.elbo_estimates, :value)

            ndraws = 2
            result3 = pathfinder(logp, ∇logp; init, ndraws, executor)
            @test size(result3.draws) == (dim, ndraws)
        end
    end
    @testset "MvNormal" begin
        #! format: off
        Σ = [
            2.71   0.5    0.19   0.07   1.04
            0.5    1.11  -0.08  -0.17  -0.08
            0.19  -0.08   0.26   0.07  -0.7
            0.07  -0.17   0.07   0.11  -0.21
            1.04  -0.08  -0.7   -0.21   8.65
        ]
        #! format: on
        P = inv(Symmetric(Σ))
        logp(x) = -dot(x, P, x) / 2
        ∇logp(x) = -(P * x)
        dim = 5
        ad_backend = AD.ReverseDiffBackend()
        ndraws_elbo = 100
        rngs = if VERSION ≥ v"1.7"
            [MersenneTwister(), Random.default_rng()]
        else
            [MersenneTwister()]
        end
        seed = 38
        optimizer = Optim.LBFGS(; m=6)
        @testset for rng in rngs
            executor = rng isa MersenneTwister ? SequentialEx() : ThreadedEx()

            Random.seed!(rng, seed)
            result = @inferred pathfinder(
                logp; rng, dim, optimizer, ndraws_elbo, ad_backend, executor
            )
            @test result.input === logp
            @test result.fit_distribution.Σ ≈ Σ rtol = 1e-1
            @test result.optimizer == optimizer
            Random.seed!(rng, seed)
            result2 = pathfinder(
                logp; rng, dim, optimizer, ndraws_elbo, ad_backend, executor
            )
            @test result2.fit_distribution == result.fit_distribution
            @test result2.draws == result.draws
            @test getproperty.(result2.elbo_estimates, :value) ==
                getproperty.(result.elbo_estimates, :value)
        end
        @testset "kwargs forwarded to solve" begin
            Random.seed!(42)
            i = 0
            callback = (args...,) -> (i += 1; false)
            pathfinder(logp; dim, callback)
            @test i ≠ 6

            Random.seed!(42)
            i = 0
            pathfinder(logp; dim, callback, maxiters=5)
            @test i == 6
        end
    end
    @testset "retries" begin
        @testset "logp returning NaN" begin
            dim = 5
            nfail = 20
            logp(x) = i ≤ nfail ? NaN : -sum(abs2, x) / 2
            callback = (args...,) -> (i += 1; true)
            i = 1
            result = pathfinder(logp; dim, callback)
            @test result.fit_distribution.μ ≈ zeros(dim) atol = 1e-6
            @test result.fit_distribution.Σ ≈ diagm(ones(dim)) atol = 1e-6
            @test result.num_tries == nfail + 1
            @test result.optim_prob.u0 == result.optim_trace.points[1]
            i = 1
            init = randn(dim)
            result2 = pathfinder(logp; init, callback, ntries=nfail)
            @test !isapprox(result2.fit_distribution.μ, zeros(dim); atol=1e-6)
            @test result2.fit_iteration == 0
            @test isempty(result2.elbo_estimates)
            @test result2.num_tries == nfail
            @test result2.optim_prob.u0 == result2.optim_trace.points[1]
        end
    end
    @testset "UniformSampler" begin
        @test_throws DomainError Pathfinder.UniformSampler(-1.0)
        @test_throws DomainError Pathfinder.UniformSampler(0.0)
        @testset for scale in [1, 2], seed in [42, 38]
            sampler_fun = Pathfinder.UniformSampler(scale)
            rng = Random.seed!(Random.default_rng(), seed)
            x = zeros(100)
            sampler_fun(rng, x)
            @test all(-scale .≤ x .≤ scale)
            Random.seed!(rng, seed)
            x2 = zeros(100)
            x2 .= rand.(rng) .* 2scale .- scale
            @test x2 == x
        end
    end
    @testset "errors if no gradient provided" begin
        logp(x) = -sum(abs2, x) / 2
        init = randn(5)
        prob = SciMLBase.OptimizationProblem(logp, init, nothing)
        @test_throws ArgumentError pathfinder(prob)
        fun = SciMLBase.OptimizationFunction(logp, GalacticOptim.AutoForwardDiff())
        prob = SciMLBase.OptimizationProblem(fun, init, nothing)
        @test_throws ArgumentError pathfinder(prob)
    end
    @testset "errors if neither dim nor init valid" begin
        logp(x) = -sum(abs2, x) / 2
        @test_throws ArgumentError pathfinder(logp)
        @test_throws ArgumentError pathfinder(logp; dim=0)
        pathfinder(logp; dim=3)
        pathfinder(logp; init=randn(5))
    end
end


text:
import Base: fill!, similar, eye, zeros, ones, fill

function fill(X::Type{<: GPUArray}, val, dims::Integer...)
    fill(X, val, dims)
end
function fill(X::Type{<: GPUArray}, val::T, dims::NTuple{N, Integer}) where {T, N}
    res = similar(X, T, dims)
    fill!(res, val)
end

function fill!(A::GPUArray{T, N}, val) where {T, N}
    valconv = T(val)
    gpu_call(const_kernel2, A, (A, valconv, UInt32(length(A))))
    A
end

zeros(T::Type{<: GPUArray}, dims::NTuple{N, Integer}) where N = fill(T, zero(eltype(T)), dims)
ones(T::Type{<: GPUArray}, dims::NTuple{N, Integer}) where N = fill(T, one(eltype(T)), dims)

function eyekernel(state, res::AbstractArray{T}, stride) where T
    i = linear_index(state)
    i > stride && return
    ilin = (stride * (i - UInt32(1))) + i
    @inbounds res[ilin] = one(T)
    return
end

eye(T::Type{<: GPUArray}, i1::Integer) = eye(T, (i1, i1))
eye(T::Type{<: GPUArray}, i1::Integer, i2::Integer) = eye(T, (i1, i2))
function eye(T::Type{<: GPUArray}, dims::NTuple{2, Integer})
    res = zeros(T, dims)
    gpu_call(eyekernel, res, (res, UInt32(size(res, 1))), minimum(dims))
    res
end

(T::Type{<: GPUArray})(dims::Integer...) = T(dims)
(T::Type{<: GPUArray{X} where X})(dims::NTuple{N, Integer}) where N = similar(T, eltype(T), dims)

similar(x::X, ::Type{T}, size::Base.Dims{N}) where {X <: GPUArray, T, N} = similar(X, T, size)

convert(AT::Type{<: GPUArray{T, N}}, A::GPUArray{T, N}) where {T, N} = A

function indexstyle(x::T) where T
    style = try
        Base.IndexStyle(x)
    catch
        nothing
    end
    style
end

function collect_kernel(state, A, iter, ::IndexCartesian)
    idx = @cartesianidx(A, state)
    @inbounds A[idx...] = iter[idx...]
    return
end

function collect_kernel(state, A, iter, ::IndexLinear)
    idx = linear_index(state)
    @inbounds A[idx] = iter[idx]
    return
end

eltype_or(::Type{<: GPUArray}, or) = or
eltype_or(::Type{<: GPUArray{T}}, or) where T = T
eltype_or(::Type{<: GPUArray{T, N}}, or) where {T, N} = T

function convert(AT::Type{<: GPUArray}, iter)
    isize = Base.iteratorsize(iter)
    style = indexstyle(iter)
    ettrait = Base.iteratoreltype(iter)
    if isbits(iter) && isize == Base.HasShape() && style != nothing && ettrait == Base.HasEltype()
        # We can collect on the GPU
        A = similar(AT, eltype_or(AT, eltype(iter)), size(iter))
        gpu_call(collect_kernel, A, (A, iter, style))
        A
    else
        convert(AT, collect(iter))
    end
end

function convert(AT::Type{<: GPUArray{T, N}}, A::DenseArray{T, N}) where {T, N}
    copy!(AT(Base.size(A)), A)
end

function convert(AT::Type{<: GPUArray{T1}}, A::DenseArray{T2, N}) where {T1, T2, N}
    copy!(similar(AT, T1, size(A)), convert(Array{T1, N}, A))
end
function convert(AT::Type{<: GPUArray}, A::DenseArray{T2, N}) where {T2, N}
    copy!(similar(AT, T2, size(A)), A)
end

function convert(AT::Type{Array{T, N}}, A::GPUArray{CT, CN}) where {T, N, CT, CN}
    convert(AT, copy!(Array{CT, CN}(Int.(Base.size(A))), A))
end


text:
module MLJIteration

using MLJBase
using IterationControl
import IterationControl: debug, skip, composite
import IterationControl: update!, done, takedown, train!

export IteratedModel

const CONTROLS = vcat(IterationControl.CONTROLS,
                      [:WithIterationsDo,
                       :WithEvaluationDo,
                       :WithFittedParamsDo,
                       :WithReportDo,
                       :WithMachineDo,
                       :WithModelDo,
                       :CycleLearningRate])

const TRAINING_CONTROLS = [:Step, ]

# export all control types:
for control in CONTROLS
    eval(:(export $control))
end

const CONTROLS_DEFAULT = [Step(10),
                          Patience(5),
                          GL(),
                          TimeLimit(0.03), # about 2 mins
                          InvalidValue()]

include("utilities.jl")
include("controls.jl")

const Control = Union{[@eval($c) for c in CONTROLS]...}
const TrainingControl = Union{[@eval($c) for c in TRAINING_CONTROLS]...}

include("constructors.jl")
include("traits.jl")
include("ic_model.jl")
include("core.jl")

end # module


text:
struct WrappingInt32
    val::UInt32
end

@kwdef mutable struct TCPHeader
    sport::UInt16 = 0
    dport::UInt16 = 0
    seqno::WrappingInt32 = WrappingInt32(0)
    ackno::WrappingInt32 = WrappingInt32(0)
    doff::UInt8 = 5
    urg::Bool = false
    ack::Bool = false
    psh::Bool = false
    rst::Bool = false
    syn::Bool = false
    fin::Bool = false
    win::UInt16 = 0
    cksum::UInt16 = 0
    uptr::UInt16 = 0
end

struct Buffer
    storage::Union{Base.RefValue{String}, Base.RefValue{SubString{String}}}
end

struct TCPSegment
    header::TCPHeader
    payload::Buffer
    TCPSegment(data::String) = new(TCPHeader(), Buffer(Ref(data)))
    TCPSegment(header::TCPHeader, data) = new(header, Buffer(Ref(data)))
end

@inline Base.length(bf::Buffer) = length(bf.storage[])
@inline length_in_sequence_space(seg::TCPSegment) = length(seg.payload.storage[]) + seg.header.syn + seg.header.fin

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


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

"""
    struct MCMCChainPoolInit <: MCMCInitAlgorithm

MCMC chain pool initialization strategy.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct MCMCChainPoolInit <: MCMCInitAlgorithm
    init_tries_per_chain::ClosedInterval{Int64} = ClosedInterval(8, 128)
    nsteps_init::Int64 = 1000
end

export MCMCChainPoolInit

function _construct_chain(
    rngpart::RNGPartition,
    id::Integer,
    algorithm::MCMCAlgorithm,
    density::AbstractDensity,
    initval_alg::InitvalAlgorithm
)
    rng = AbstractRNG(rngpart, id)
    v_init = unshaped(bat_initval(rng, density, initval_alg).result, varshape(density))
    MCMCIterator(rng, algorithm, density, id, v_init)
end

_gen_chains(
    rngpart::RNGPartition,
    ids::AbstractRange{<:Integer},
    algorithm::MCMCAlgorithm,
    density::AbstractDensity,
    initval_alg::InitvalAlgorithm
) = [_construct_chain(rngpart, id, algorithm, density, initval_alg) for id in ids]

function mcmc_init!(
    rng::AbstractRNG,
    algorithm::MCMCAlgorithm,
    density::AbstractDensity,
    nchains::Integer,
    init_alg::MCMCChainPoolInit,
    tuning_alg::MCMCTuningAlgorithm,
    nonzero_weights::Bool,
    callback::Function
)
    @info "Trying to generate $nchains viable MCMC chain(s)."

    initval_alg = InitFromTarget()

    min_nviable::Int = minimum(init_alg.init_tries_per_chain) * nchains
    max_ncandidates::Int = maximum(init_alg.init_tries_per_chain) * nchains

    rngpart = RNGPartition(rng, Base.OneTo(max_ncandidates))

    ncandidates::Int = 0

    dummy_initval = unshaped(bat_initval(rng, density, InitFromTarget()).result, varshape(density))
    dummy_chain = MCMCIterator(deepcopy(rng), algorithm, density, 1, dummy_initval)
    dummy_tuner = tuning_alg(dummy_chain)

    chains = similar([dummy_chain], 0)
    tuners = similar([dummy_tuner], 0)
    outputs = similar([DensitySampleVector(dummy_chain)], 0)
    cycle::Int = 1

    while length(tuners) < min_nviable && ncandidates < max_ncandidates
        n = min(min_nviable, max_ncandidates - ncandidates)
        @debug "Generating $n $(cycle > 1 ? "additional " : "")MCMC chain(s)."

        new_chains = _gen_chains(rngpart, ncandidates .+ (one(Int64):n), algorithm, density, initval_alg)

        filter!(isvalidchain, new_chains)

        new_tuners = tuning_alg.(new_chains)
        new_outputs = DensitySampleVector.(new_chains)
        next_cycle!.(new_chains)
        tuning_init!.(new_tuners, new_chains, init_alg.nsteps_init)
        ncandidates += n

        @debug "Testing $(length(new_tuners)) MCMC chain(s)."

        mcmc_iterate!(
            new_outputs, new_chains, new_tuners;
            max_nsteps = clamp(div(init_alg.nsteps_init, 5), 10, 50),
            callback = callback,
            nonzero_weights = nonzero_weights
        )

        viable_idxs = findall(isviablechain.(new_chains))
        viable_tuners = new_tuners[viable_idxs]
        viable_chains = new_chains[viable_idxs]
        viable_outputs = new_outputs[viable_idxs]

        @debug "Found $(length(viable_idxs)) viable MCMC chain(s)."

        if !isempty(viable_tuners)
            mcmc_iterate!(
                viable_outputs, viable_chains, viable_tuners;
                max_nsteps = init_alg.nsteps_init,
                callback = callback,
                nonzero_weights = nonzero_weights
            )

            nsamples_thresh = floor(Int, 0.8 * median([nsamples(chain) for chain in viable_chains]))
            good_idxs = findall(chain -> nsamples(chain) >= nsamples_thresh, viable_chains)
            @debug "Found $(length(viable_tuners)) MCMC chain(s) with at least $(nsamples_thresh) unique accepted samples."

            append!(chains, view(viable_chains, good_idxs))
            append!(tuners, view(viable_tuners, good_idxs))
            append!(outputs, view(viable_outputs, good_idxs))
        end

        cycle += 1
    end

    length(tuners) < min_nviable && error("Failed to generate $min_nviable viable MCMC chains")

    m = nchains
    tidxs = LinearIndices(tuners)
    n = length(tidxs)

    modes = hcat(broadcast(samples -> Array(bat_findmode(samples, MaxDensitySampleSearch()).result), outputs)...)

    final_chains = similar(chains, 0)
    final_tuners = similar(tuners, 0)
    final_outputs = similar(outputs, 0)

    if 2 <= m < size(modes, 2)
        clusters = kmeans(modes, m, init = KmCentralityAlg())
        clusters.converged || error("k-means clustering of MCMC chains did not converge")

        mincosts = fill(Inf, m)
        chain_sel_idxs = fill(0, m)

        for i in tidxs
            j = clusters.assignments[i]
            if clusters.costs[i] < mincosts[j]
                mincosts[j] = clusters.costs[i]
                chain_sel_idxs[j] = i
            end
        end

        @assert all(j -> j in tidxs, chain_sel_idxs)

        for i in sort(chain_sel_idxs)
            push!(final_chains, chains[i])
            push!(final_tuners, tuners[i])
            push!(final_outputs, outputs[i])
        end
    else
        @assert length(chains) == nchains
        resize!(final_chains, nchains)
        copyto!(final_chains, chains)

        @assert length(tuners) == nchains
        resize!(final_tuners, nchains)
        copyto!(final_tuners, outputs)

        @assert length(outputs) == nchains
        resize!(final_outputs, nchains)
        copyto!(final_outputs, outputs)
    end

    @info "Selected $(length(final_tuners)) MCMC chain(s)."
    tuning_postinit!.(final_tuners, final_chains, final_outputs)

    (chains = final_chains, tuners = final_tuners, outputs = final_outputs)
end


text:
######
##FEM Poisson dx Convergence Tests
######
using FiniteElementDiffEq, DiffEqDevTools,DiffEqProblemLibrary

cs = cs_fempoisson_wave

sim = test_convergence(cs,FEMDiffEqPoisson())

#Plot Result
#plot(sim,xguide="Delta x")

#Returns true if convergence is like dx^2 in L2
sim.𝒪est[:L2]-2 <.1


text:
export runSimulation!, simulate!

@doc doc"""
Function for running a simulation from the current model definition and parsing the output data into new SimulationData-structure

    ModelInstance = runsimulation(ModelInstance)
"""->
function runSimulation!(M::Model, name; plot = false, pars = false, vars = false, xlim = false, ylim =false, colors = false, linewidth = 2)
    # Save current version of mpdel to odefile
    toOdeFile(M)
    #Update the vars list to match the order in the ode file
    M.vars = getVariables(M)
    #Run the simulation
    odefile = M.name * ".ode"
    xppcall = os[OS_NAME].xppcall
    options = os[OS_NAME].options
    run(`$xppcall $odefile $options`)
    #Open and parse output file into new SimulationData-structure
    f = open("output.dat")
    M = parseOutputFile(f, M, name)
    close(f)
    if plot
        plotModel(M, name; pars = pars, vars = vars, xlim = xlim, ylim = ylim, colors =colors, linewidth = linewidth)
    end
    run(`rm output.dat`)
    n = M.name
    run(`rm $n.ode`)
    run(`rm log.txt`)
end

"""
Simulate the system using ODE.jl package
"""
function simulate!(M::Model, name, trange; solver = ode23s, plot = false, pars = false, vars = false, xlim = false, ylim =false, colors = false, linewidth = 2)
    #Update parameters
    M.p = Float64[M.pars[p] for p in M.p_names]
    #Update initials
    M.y0 = Float64[M.init[y] for y in M.y_names]
    #Convert F into necessary form (only dependent on t and y)
    # for (i,n) in enumerate(M.p_names)
    #   println(n, "\t", p[i])
    # end
    F(t,y) = M.F(t, y, M.p)
    #Uses ode23s by default for other available solver consult ODE module: https://github.com/JuliaLang/ODE.jl
    t,y = solver(F, M.y0, trange)
    if plot != false
        plotModel(M, name; pars = pars, vars = vars, xlim = xlim, ylim = ylim, colors =colors, linewidth = linewidth)
    end
    #Retrieve results
    M = parseODEOutput(t, y, M, name)
    return(M)
end

"""
Get vectorised output and store it in SimulationData structure
"""
function parseODEOutput(t, ysim, M::Model, name)
    if name == false
        #Get the new key for the dict
        k = length(M.sims) + 1
    else
        #Overwrite the last simulation
        k = name
    end
    #Instantiate new SimulationData-structure
    M.sims[k] = SimulationData(M, name)
    M.sims[k].D["t"] = t
    for (i,y) in enumerate(M.y_names)
      M.sims[k].D[y] = [y[i] for y in ysim]
    end

    for (i,a) in enumerate(M.a_names)
      M.sims[k].D[a] = zeros(length(t))
    end
    for (ix,ti) in enumerate(t)
      algs = M.A(ti, ysim[ix], M.p)
      for (i,a) in enumerate(M.a_names)
        M.sims[k].D[a][ix] = deepcopy(algs[i])
      end
    end
    return(M)
end

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl