text:
#: ========================================================================
#global const __DEFAULT__R__ = (__DEFAULT_HASH_CODE__, TPattern[])  #+ modify
#global const SIMILARITY_LEVEL = 0.9
#: ========================================================================

@inline processelemif(p,c,L) = for el ∈ L if c(el) p(el); end end

global const exp_m_i = [exp(-i) for i=1:100]
global const coeff_norm = [1.0/((1-exp(-i))/(ℯ-1)) for i=1:100]

#@inline exp_weighted(f,X::Int) = coeff_norm[X]*mapreduce(i->exp_m_i[i]*f(i), +, 1:X)
#@inline exp_weighted(f,X::Int) = coeff_norm[X]*sum([exp_m_i[i]*f(i) for i=1:X])
function exp_weighted(f::Function,X::Int)
    s = 0.0
    for i=1:X
        @inbounds s += exp_m_i[i]*f(i)
    end
    return coeff_norm[X]*s
end

@inline first_elem_iden_weighted(f,X::Int) = (f(1) ? exp_weighted(f,X) : 0.0)

@inline last_elem_sim_weighted(f,X::Int) = (f(X)>0.8 ? exp_weighted(f,X) : 0.0)

#TODO optimize
function patt_similarity2(p1::TPattern, p2::TPattern)::Float64
    len1 = length(p1)
    len2 = length(p2)
    if len1==len2==0  return 1.0  end
    lmin = min(len1,len2)
    return (lmin==0 ? 0.0 : (first_elem_iden_weighted(i->(p1[i]==p2[i]), lmin)))
end

#TODO optimize
function similarity2(a,b)::Float64
    if a[1]==b[1]
        return 1.0
    else
        N = length(a[2])
        if N!=length(b[2])
            return 0.0
        else
            return exp_weighted(i->patt_similarity2(a[2][i],b[2][i]), N)
        end
    end
end

##
#p1 = (:x,TPattern[[13, 0, 0, 0, 0, 0, 0], [14, 15, 0], [16, 17, 18, 19], [20, 21, 0], [-1]])
#p2 = (:y,TPattern[[22, 0, 0, 23], [14, 15, 0], [16, 17, 18, 19], [20, 21, 0], [-1]])
#similarity2(p1,p2)
#patt_similarity2([13, 0, 0, 0, 0, 0, 0],[22, 0, 0, 23])
#patt_similarity2([10, 11, 0, 12],[10, 11, 0, 0, 18, 19, 20, 17, 21, 22])^2
##

function good_label_crit(q::Int, A::Vector)::Bool
    if length(A[q][2][1])>1  return true  end
    return !(@inbounds A[q][2][1]==__M0_PATT__::TPattern || A[q][2][1]==__PATT__all_number_line__::TPattern)
end

# search substr in fullstr FROM THE RIGHTEST POSITION
function nonoverlapping(
    substr::Vector{TR},
    fullstr::Vector{TR},
    fullstr_good_pos::Vector{Bool}
    )::Vector{IntRange}  where  TR

    Nsub    = length(substr)
    sim_min = Nsub * SIMILARITY_LEVEL
    sim     = -1.0
    q       = -1
    posL    = length(fullstr)-Nsub+1
    R       = IntRange[]
    j       = -1
    while posL >= 1
        q = -1
        # test similarity FROM THE RIGHTEST POSITION
        for p ∈ posL:-1:1
            #if good_label_crit(p,fullstr)
            if fullstr_good_pos[p]
                @inbounds sim = similarity2(substr[1], fullstr[p]) + similarity2(substr[Nsub], fullstr[p+Nsub-1])
                #=
                if Nsub==2 && sim>1.9
                    println("Nsub==2 && sim>1.9")
                    @show sim
                    @show substr
                    @show fullstr[p:p+Nsub-1]
                    @show similarity2(fullstr[p], fullstr[p+Nsub-1])
                    println("")
                end
                =#
                if sim<1.5 || similarity2(fullstr[p], fullstr[p+Nsub-1])>0.4
                    # if first and last pattern doesn't match at high fidelity
                    # or the first and last in the fullstr to be matched are too similar
                    continue
                end
                # going FROM THE LEFTEST POSITION
                # j = 2 to Nsub-1
                j = 2
                while j<=Nsub-1
                    @inbounds sim += similarity2(substr[j], fullstr[p+j-1])
                    if (sim+(Nsub-j)+1e-12<sim_min)
                        # if, even the remaining part all matches exactly (sim+=1 for each)
                        # the sim value still below minimum
                        break
                    end
                    j += 1 # going FROM THE LEFTEST POSITION
                end
                if j==Nsub
                    # the above while loop reaches the end
                    # overall similarity larger than fixed level
                    q = p
                    push!(R, q:q+Nsub-1)
                    posL = q-Nsub
                    break
                end
            end
        end
        if q < 0
            break
        end
    end
    return R
end

function MostFreqSimilarSubsq(
    str::Vector{TR};
    # meaning of Lmin and Lmax :
    # min and max lengths of the pattern for
    # a group of blocks
    # experiments show that Lmin=3 is efficient for QE files
    Lmin=2,
    Lmax=20
    )::Vector{IntRange} where TR

    #: output
    # println("MostFreqSimilarSubsq:")

    N = length(str)
    if N<3 || allunique(str)
        return IntRange[]
    end

    #+------------- lift the "degenereacy" --------------
    # most apperance >> earliest appearance
    sortf1(x) = 100*length(x) - (first(last(x))/N)
    # most apperance >> longest range
    sortf2(x) = 1000*length(x) + length(first(x))
    function crit_i_l(i::Int, l::Int)
        if length(str[i][2][1])>1   return true   end
        @inbounds cnd = (str[i][2][1]==__M0_PATT__::TPattern || str[i][2][1]==__PATT__all_number_line__::TPattern)
        if cnd
            return false
        else
            return (@inbounds similarity2(str[i],str[i+l-1])<0.5)
        end
    end
    #+---------------------------------------------------

    good_p_label= Bool[good_label_crit(i,str) for i=1:N]
    RES         = Vector{IntRange}[]
    B           = nothing
    lenLASTBmax = -10

    all_blk_l   = Vector{IntRange}[]
    max_len_B   = -1
    curr_len_B  = 0
    lB          = 0
    for l=Lmin:min(Lmax,N÷2)
        U = Dict{Vector{TR},Int}()  # record the unique substrings by num of reps
        updU(i::Int) = (@inbounds increaseindex!(U, str[i:i+l-1]))  # dict for num of repetitions
        crit_i(i::Int) = crit_i_l(i, l)
        @inbounds processelemif(updU, crit_i, 1:N-l+1) ;

        all_blk_l  = Vector{IntRange}[]
        max_len_B  = 2
        curr_len_B = 0
        lB         = 0
        for (s,m) ∈ U
            if m>1
                B  = nonoverlapping(s, str, good_p_label)
                lB = length(B)
                if lB>=max_len_B
                    max_len_B = lB
                    push!(all_blk_l, B)
                end
            end
        end
        if max_len_B < lenLASTBmax  # increasing the sub-pattern length l
            break  # won't get longer nonoverlapping blocks B
        end
        if length(all_blk_l)>0
            (_,_i) = findmax(sortf1.(all_blk_l))
            lenLASTBmax = max(lenLASTBmax, length(all_blk_l[_i]))
            push!(RES, all_blk_l[_i])  # only record the best for each l
        end

        #: output
        # println("------------ l = $l ------------")
        # A = sort(all_blk_l,by=sortf1)
        # [(println((length(A[k]),last(A[k]), str[last(A[k])])); 0)  for k=max(1,length(A)-10):length(A)]

    end

    #: output
    # println("-------------------------------")
    # println("")
    # println("")

    if length(RES)==0  return IntRange[]  end
    # only return the best
    (_,_i) = findmax(sortf2.(RES))
    return RES[_i]
end

#=

function MostFreqSimilarSubsq_test_ver(
    str::Vector{TR};
    # meaning of Lmin and Lmax :
    # min and max lengths of the pattern for
    # a group of blocks
    # experiments show that Lmin=3 is efficient for QE files
    Lmin=3,
    Lmax=20
    )::Vector{IntRange} where TR

    N = length(str)
    if N<3 || allunique(str)
        return IntRange[]
    end

    #+------------- lift the "degenereacy" --------------
    # most apperance >> earliest appearance
    sortf1(x) = length(x) - (first(last(x))/N)
    # most apperance >> longest range
    sortf2(x) = 10*length(x) + length(first(x))
    #+---------------------------------------------------

    RES = Vector{IntRange}[]
    B   = nothing
    lenLASTBmax = -10
    all_blk_l = Vector{IntRange}[]
    max_len_B  = -1
    curr_len_B = 0
    lB  = 0
    sortf1_val = 0.0
    for l=Lmin:min(Lmax,N÷2)
        #println("------ l = $l ----------")
        U = Dict{Vector{TR},Int}()  # record the unique substrings by num of reps
        updU(i)   = (@inbounds increaseindex!(U, str[i:i+l-1]))  # dict for num of repetitions
        crit_i(i) = (@inbounds  good_label_crit(str[i]) && similarity2(str[i],str[i+l-1])<0.5)
        processelemif(updU, crit_i, 1:N-l+1) ;

        all_blk_l = Vector{IntRange}[]
        max_len_B = 2
        curr_len_B = 0
        lB = 0
        for (s,m) ∈ U
            if m>1
                B  = nonoverlapping(s, str)
                lB = length(B)
                #if lB>1 && lB>=max_len_B
                if lB>=max_len_B
                    max_len_B = lB
                    push!(all_blk_l, B)
                end
            end
        end
        if max_len_B < lenLASTBmax
            # increasing the sub-pattern length l
            # won't get longer nonoverlapping blocks B
            break
        end
        if length(all_blk_l)>0
            #A = sort(all_blk_l, by=sortf1)
            #LASTB = last(A)
            #LASTB = last(sort(all_blk_l, by=sortf1))
            #@assert max_len_B==length(LASTB)
            (_,LASTB_i) = findmax(sortf1.(all_blk_l))
            curr_len_B = length(all_blk_l[LASTB_i])
            @assert max_len_B==curr_len_B
            #[(println((length(A[k]),last(A[k]), str[last(A[k])])); 0)  for k=max(1,length(A)-10):length(A)]
            #lenLASTBmax = max(lenLASTBmax, length(LASTB))
            lenLASTBmax = max(lenLASTBmax, curr_len_B)
            #push!(RES, LASTB)  # only record the best for each l
            push!(RES, all_blk_l[LASTB_i])  # only record the best for each l
        end
    end
    return length(RES)==0 ? IntRange[] : last(sort(RES,by=sortf2))  # only return the best
end

# use heap , #! slow
function MostFreqSimilarSubsq_heap_ver(
    str::Vector{TR};
    # meaning of Lmin and Lmax :
    # min and max lengths of the pattern for
    # a group of blocks
    # experiments show that Lmin=3 is efficient for QE files
    Lmin=3,
    Lmax=20
    )::Vector{IntRange} where TR

    N = length(str)

    if N<3 || allunique(str)
        return IntRange[]
    end

    #+------------- lift the "degenereacy" --------------
    # most apperance >> longest range >> earliest appearance
    sortf3(x::Vector{IntRange}) = -(2Lmax*length(x) + length(first(x)) - (first(last(x))/N))
    #+---------------------------------------------------

    B        = nothing
    l_B_max  = -10
    lB       = 0
    blk_H    = BinaryHeap(Base.By(sortf3), Vector{IntRange}[])
    for l=Lmin:min(Lmax,N÷2)
        l_B_max_at_l = -1
        U = Dict{Vector{TR},Int}()  # record the unique substrings by num of reps
        updU(i::Int)   = (@inbounds increaseindex!(U, str[i:i+l-1]))  # dict for num of repetitions
        crit_i(i::Int) = (@inbounds  (good_label_crit(str[i]) && similarity2(str[i],str[i+l-1])<0.5))
        processelemif(updU, crit_i, 1:N-l+1) ;

        lB = 0
        for (s,m) ∈ U
            if m>1
                B  = nonoverlapping(s, str)
                lB = length(B)
                if lB>1 && lB>=l_B_max_at_l
                    l_B_max_at_l = lB
                    push!(blk_H, B)
                end
            end
        end
        if l_B_max_at_l < l_B_max
            # increasing the sub-pattern length l
            # won't get longer nonoverlapping blocks B
            break
        end
    end
    return length(blk_H)==0 ? IntRange[] : first(blk_H)  # only return the best
end

=#

#=

# not used
avg(l) = sum(l)/length(l)

# not used
function patt_similarity(p1::TPattern, p2::TPattern, f=identity)::Float64
    len1 = length(p1)
    len2 = length(p2)
    lmin = min(len1,len2)
    return (len1==len2==0 ? 1.0 : f(avg(p1[1:lmin].==p2[1:lmin])))
end

# not used
similarity(a,b)::Float64 = (a[1]==b[1] ? 1.0 : (length(a[2])!=length(b[2]) ? 0.0 : avg([patt_similarity(x,y) for (x,y) ∈ zip(a[2],b[2])])))

# not used
function sim_each_char(s1::Vector{TR}, i1::Int, s2::Vector{TR}, i2::Int)::Bool where TR
    n1 = length(s1)
    j = 0
    while j < n1
        if similarity(s1[i1+j], s2[i2+j]) < SIMILARITY_LEVEL   return false   end
        j += 1
    end
    return j==n1
end

=#

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
const zhuyin_trie = Trie{String}([
    ("a", "ㄚ"),
    ("ai", "ㄞ"),
    ("an", "ㄢ"),
    ("ang", "ㄤ"),
    ("ao", "ㄠ"),
    ("ba", "ㄅㄚ"),
    ("bai", "ㄅㄞ"),
    ("ban", "ㄅㄢ"),
    ("bang", "ㄅㄤ"),
    ("bao", "ㄅㄠ"),
    ("bei", "ㄅㄟ"),
    ("ben", "ㄅㄣ"),
    ("beng", "ㄅㄥ"),
    ("bi", "ㄅㄧ"),
    ("bian", "ㄅㄧㄢ"),
    ("biao", "ㄅㄧㄠ"),
    ("bie", "ㄅㄧㄝ"),
    ("bin", "ㄅㄧㄣ"),
    ("bing", "ㄅㄧㄥ"),
    ("bo", "ㄅㄛ"),
    ("bu", "ㄅㄨ"),
    ("ca", "ㄘㄚ"),
    ("cai", "ㄘㄞ"),
    ("can", "ㄘㄢ"),
    ("cang", "ㄘㄤ"),
    ("cao", "ㄘㄠ"),
    ("ce", "ㄘㄜ"),
    ("cen", "ㄘㄣ"),
    ("ceng", "ㄘㄥ"),
    ("cha", "ㄔㄚ"),
    ("chai", "ㄔㄞ"),
    ("chan", "ㄔㄢ"),
    ("chang", "ㄔㄤ"),
    ("chao", "ㄔㄠ"),
    ("che", "ㄔㄜ"),
    ("chen", "ㄔㄣ"),
    ("cheng", "ㄔㄥ"),
    ("chi", "ㄔ"),
    ("chong", "ㄔㄨㄥ"),
    ("chou", "ㄔㄡ"),
    ("chu", "ㄔㄨ"),
    ("chua", "ㄔㄨㄚ"),
    ("chuai", "ㄔㄨㄞ"),
    ("chuan", "ㄔㄨㄢ"),
    ("chuang", "ㄔㄨㄤ"),
    ("chui", "ㄔㄨㄟ"),
    ("chun", "ㄔㄨㄣ"),
    ("chuo", "ㄔㄨㄛ"),
    ("ci", "ㄘ"),
    ("cong", "ㄘㄨㄥ"),
    ("cou", "ㄘㄡ"),
    ("cu", "ㄘㄨ"),
    ("cuan", "ㄘㄨㄢ"),
    ("cui", "ㄘㄨㄟ"),
    ("cun", "ㄘㄨㄣ"),
    ("cuo", "ㄘㄨㄛ"),
    ("da", "ㄉㄚ"),
    ("dai", "ㄉㄞ"),
    ("dan", "ㄉㄢ"),
    ("dang", "ㄉㄤ"),
    ("dao", "ㄉㄠ"),
    ("de", "ㄉㄜ"),
    ("deng", "ㄉㄥ"),
    ("di", "ㄉㄧ"),
    ("dian", "ㄉㄧㄢ"),
    ("diang", "ㄉㄧㄤ"),
    ("diao", "ㄉㄧㄠ"),
    ("die", "ㄉㄧㄝ"),
    ("ding", "ㄉㄧㄥ"),
    ("diu", "ㄉㄧㄡ"),
    ("dong", "ㄉㄨㄥ"),
    ("dou", "ㄉㄡ"),
    ("du", "ㄉㄨ"),
    ("duan", "ㄉㄨㄢ"),
    ("dui", "ㄉㄨㄟ"),
    ("dun", "ㄉㄨㄣ"),
    ("duo", "ㄉㄨㄛ"),
    ("e", "ㄜ"),
    ("ei", "ㄟ"),
    ("en", "ㄣ"),
    ("er", "ㄦ"),
    ("fa", "ㄈㄚ"),
    ("fan", "ㄈㄢ"),
    ("fang", "ㄈㄤ"),
    ("fei", "ㄈㄟ"),
    ("fen", "ㄈㄣ"),
    ("feng", "ㄈㄥ"),
    ("fo", "ㄈㄛ"),
    ("fou", "ㄈㄡ"),
    ("fu", "ㄈㄨ"),
    ("ga", "ㄍㄚ"),
    ("gai", "ㄍㄞ"),
    ("gan", "ㄍㄢ"),
    ("gang", "ㄍㄤ"),
    ("gao", "ㄍㄠ"),
    ("ge", "ㄍㄜ"),
    ("gei", "ㄍㄟ"),
    ("gen", "ㄍㄣ"),
    ("geng", "ㄍㄥ"),
    ("gong", "ㄍㄨㄥ"),
    ("gou", "ㄍㄡ"),
    ("gu", "ㄍㄨ"),
    ("gua", "ㄍㄨㄚ"),
    ("guai", "ㄍㄨㄞ"),
    ("guan", "ㄍㄨㄢ"),
    ("guang", "ㄍㄨㄤ"),
    ("gui", "ㄍㄨㄟ"),
    ("gun", "ㄍㄨㄣ"),
    ("guo", "ㄍㄨㄛ"),
    ("ha", "ㄏㄚ"),
    ("hai", "ㄏㄞ"),
    ("han", "ㄏㄢ"),
    ("hang", "ㄏㄤ"),
    ("hao", "ㄏㄠ"),
    ("he", "ㄏㄜ"),
    ("hei", "ㄏㄟ"),
    ("hen", "ㄏㄣ"),
    ("heng", "ㄏㄥ"),
    ("hong", "ㄏㄨㄥ"),
    ("hou", "ㄏㄡ"),
    ("hu", "ㄏㄨ"),
    ("hua", "ㄏㄨㄚ"),
    ("huai", "ㄏㄨㄞ"),
    ("huan", "ㄏㄨㄢ"),
    ("huang", "ㄏㄨㄤ"),
    ("hui", "ㄏㄨㄟ"),
    ("hun", "ㄏㄨㄣ"),
    ("huo", "ㄏㄨㄛ"),
    ("ji", "ㄐㄧ"),
    ("jia", "ㄐㄧㄚ"),
    ("jian", "ㄐㄧㄢ"),
    ("jiang", "ㄐㄧㄤ"),
    ("jiao", "ㄐㄧㄠ"),
    ("jie", "ㄐㄧㄝ"),
    ("jin", "ㄐㄧㄣ"),
    ("jing", "ㄐㄧㄥ"),
    ("jiong", "ㄐㄩㄥ"),
    ("jiu", "ㄐㄧㄡ"),
    ("ju", "ㄐㄩ"),
    ("juan", "ㄐㄩㄢ"),
    ("jue", "ㄐㄩㄝ"),
    ("jun", "ㄐㄩㄣ"),
    ("ka", "ㄎㄚ"),
    ("kai", "ㄎㄞ"),
    ("kan", "ㄎㄢ"),
    ("kang", "ㄎㄤ"),
    ("kao", "ㄎㄠ"),
    ("ke", "ㄎㄜ"),
    ("ken", "ㄎㄣ"),
    ("keng", "ㄎㄥ"),
    ("kong", "ㄎㄨㄥ"),
    ("kou", "ㄎㄡ"),
    ("ku", "ㄎㄨ"),
    ("kua", "ㄎㄨㄚ"),
    ("kuai", "ㄎㄨㄞ"),
    ("kuan", "ㄎㄨㄢ"),
    ("kuang", "ㄎㄨㄤ"),
    ("kui", "ㄎㄨㄟ"),
    ("kun", "ㄎㄨㄣ"),
    ("kuo", "ㄎㄨㄛ"),
    ("la", "ㄌㄚ"),
    ("lai", "ㄌㄞ"),
    ("lan", "ㄌㄢ"),
    ("lang", "ㄌㄤ"),
    ("lao", "ㄌㄠ"),
    ("le", "ㄌㄜ"),
    ("lei", "ㄌㄟ"),
    ("leng", "ㄌㄥ"),
    ("li", "ㄌㄧ"),
    ("lian", "ㄌㄧㄢ"),
    ("liang", "ㄌㄧㄤ"),
    ("liao", "ㄌㄧㄠ"),
    ("lie", "ㄌㄧㄝ"),
    ("lin", "ㄌㄧㄣ"),
    ("ling", "ㄌㄧㄥ"),
    ("liu", "ㄌㄧㄡ"),
    ("long", "ㄌㄨㄥ"),
    ("lou", "ㄌㄡ"),
    ("lu", "ㄌㄨ"),
    ("lü", "ㄌㄩ"),
    ("lu:", "ㄌㄩ"),  # manual addition for alternate typing
    ("luan", "ㄌㄨㄢ"),
    ("lüe", "ㄌㄩㄝ"),
    ("lu:e", "ㄌㄩㄝ"),  # manual addition for alternate typing
    ("lun", "ㄌㄨㄣ"),
    ("luo", "ㄌㄨㄛ"),
    ("ma", "ㄇㄚ"),
    ("mai", "ㄇㄞ"),
    ("man", "ㄇㄢ"),
    ("mang", "ㄇㄤ"),
    ("mao", "ㄇㄠ"),
    ("mei", "ㄇㄟ"),
    ("men", "ㄇㄣ"),
    ("meng", "ㄇㄥ"),
    ("mi", "ㄇㄧ"),
    ("mian", "ㄇㄧㄢ"),
    ("miao", "ㄇㄧㄠ"),
    ("mie", "ㄇㄧㄝ"),
    ("min", "ㄇㄧㄣ"),
    ("ming", "ㄇㄧㄥ"),
    ("miu", "ㄇㄧㄡ"),
    ("mo", "ㄇㄨㄛ"),
    ("mou", "ㄇㄡ"),
    ("mu", "ㄇㄨ"),
    ("na", "ㄋㄚ"),
    ("nai", "ㄋㄞ"),
    ("nan", "ㄋㄢ"),
    ("nang", "ㄋㄤ"),
    ("nao", "ㄋㄠ"),
    ("ne", "ㄋㄛ"),
    ("nei", "ㄋㄟ"),
    ("nen", "ㄋㄣ"),
    ("neng", "ㄋㄥ"),
    ("ni", "ㄋㄧ"),
    ("nia", "ㄋㄧㄚ"),
    ("nian", "ㄋㄧㄢ"),
    ("niang", "ㄋㄧㄤ"),
    ("niao", "ㄋㄧㄠ"),
    ("nie", "ㄋㄧㄝ"),
    ("nin", "ㄋㄧㄣ"),
    ("ning", "ㄋㄧㄥ"),
    ("niu", "ㄋㄧㄡ"),
    ("nong", "ㄋㄨㄥ"),
    ("nou", "ㄋㄡ"),
    ("nu", "ㄋㄨ"),
    ("nü", "ㄋㄩ"),
    ("nuan", "ㄋㄨㄢ"),
    ("nüe", "ㄋㄩㄝ"),
    ("nuo", "ㄋㄨㄛ"),
    ("nun", "ㄋㄨㄣ"),
    ("ou", "ㄡ"),
    ("pa", "ㄆㄚ"),
    ("pai", "ㄆㄞ"),
    ("pan", "ㄆㄢ"),
    ("pang", "ㄆㄤ"),
    ("pao", "ㄆㄠ"),
    ("pei", "ㄆㄟ"),
    ("pen", "ㄆㄣ"),
    ("peng", "ㄆㄥ"),
    ("pi", "ㄆㄧ"),
    ("pian", "ㄆㄧㄢ"),
    ("piao", "ㄆㄧㄠ"),
    ("pie", "ㄆㄧㄝ"),
    ("pin", "ㄆㄧㄣ"),
    ("ping", "ㄆㄧㄥ"),
    ("po", "ㄆㄨㄛ"),
    ("pou", "ㄆㄡ"),
    ("pu", "ㄆㄨ"),
    ("qi", "ㄑㄧ"),
    ("qia", "ㄑㄧㄚ"),
    ("qian", "ㄑㄧㄢ"),
    ("qiang", "ㄑㄧㄤ"),
    ("qiao", "ㄑㄧㄠ"),
    ("qie", "ㄑㄧㄝ"),
    ("qin", "ㄑㄧㄣ"),
    ("qing", "ㄑㄧㄥ"),
    ("qiong", "ㄑㄩㄥ"),
    ("qiu", "ㄑㄧㄡ"),
    ("qu", "ㄑㄩ"),
    ("quan", "ㄑㄩㄢ"),
    ("que", "ㄑㄩㄝ"),
    ("qun", "ㄑㄩㄣ"),
    ("ran", "ㄖㄢ"),
    ("rang", "ㄖㄤ"),
    ("rao", "ㄖㄠ"),
    ("re", "ㄖㄜ"),
    ("ren", "ㄖㄣ"),
    ("reng", "ㄖㄥ"),
    ("ri", "ㄖ"),
    ("rong", "ㄖㄨㄥ"),
    ("rou", "ㄖㄡ"),
    ("ru", "ㄖㄨ"),
    ("ruan", "ㄖㄨㄢ"),
    ("rui", "ㄖㄨㄟ"),
    ("run", "ㄖㄨㄣ"),
    ("ruo", "ㄖㄨㄛ"),
    ("sa", "ㄙㄚ"),
    ("sai", "ㄙㄞ"),
    ("san", "ㄙㄢ"),
    ("sang", "ㄙㄤ"),
    ("sao", "ㄙㄠ"),
    ("se", "ㄙㄜ"),
    ("sei", "ㄙㄟ"),
    ("sen", "ㄙㄣ"),
    ("seng", "ㄙㄥ"),
    ("sha", "ㄕㄚ"),
    ("shai", "ㄕㄞ"),
    ("shan", "ㄕㄢ"),
    ("shang", "ㄕㄤ"),
    ("shao", "ㄕㄠ"),
    ("she", "ㄕㄜ"),
    ("shei", "ㄕㄟ"),
    ("shen", "ㄕㄣ"),
    ("sheng", "ㄕㄥ"),
    ("shi", "ㄕ"),
    ("shong", "ㄕㄡㄥ"),
    ("shou", "ㄕㄡ"),
    ("shu", "ㄕㄨ"),
    ("shua", "ㄕㄨㄚ"),
    ("shuai", "ㄕㄨㄞ"),
    ("shuan", "ㄕㄨㄢ"),
    ("shuang", "ㄕㄨㄤ"),
    ("shui", "ㄕㄨㄟ"),
    ("shun", "ㄕㄨㄣ"),
    ("shuo", "ㄕㄨㄛ"),
    ("si", "ㄙ"),
    ("song", "ㄙㄨㄥ"),
    ("sou", "ㄙㄡ"),
    ("su", "ㄙㄨ"),
    ("suan", "ㄙㄨㄢ"),
    ("sui", "ㄙㄨㄟ"),
    ("sun", "ㄙㄨㄣ"),
    ("suo", "ㄙㄨㄛ"),
    ("ta", "ㄊㄚ"),
    ("tai", "ㄊㄞ"),
    ("tan", "ㄊㄢ"),
    ("tang", "ㄊㄤ"),
    ("tao", "ㄊㄠ"),
    ("te", "ㄊㄜ"),
    ("teng", "ㄊㄥ"),
    ("ti", "ㄊㄧ"),
    ("tian", "ㄊㄧㄢ"),
    ("tiao", "ㄊㄧㄠ"),
    ("tie", "ㄊㄧㄝ"),
    ("ting", "ㄊㄧㄥ"),
    ("tong", "ㄊㄨㄥ"),
    ("tou", "ㄊㄡ"),
    ("tu", "ㄊㄨ"),
    ("tuan", "ㄊㄨㄢ"),
    ("tui", "ㄊㄨㄟ"),
    ("tun", "ㄊㄨㄣ"),
    ("tuo", "ㄊㄨㄛ"),
    ("wa", "ㄨㄚ"),
    ("wai", "ㄨㄞ"),
    ("wan", "ㄨㄢ"),
    ("wang", "ㄨㄤ"),
    ("wei", "ㄨㄟ"),
    ("wen", "ㄨㄣ"),
    ("weng", "ㄨㄥ"),
    ("wo", "ㄨㄛ"),
    ("wu", "ㄨ"),
    ("xi", "ㄒㄧ"),
    ("xia", "ㄒㄧㄚ"),
    ("xian", "ㄒㄧㄢ"),
    ("xiang", "ㄒㄧㄤ"),
    ("xiao", "ㄒㄧㄠ"),
    ("xie", "ㄒㄧㄝ"),
    ("xin", "ㄒㄧㄣ"),
    ("xing", "ㄒㄧㄥ"),
    ("xiong", "ㄒㄩㄥ"),
    ("xiu", "ㄒㄧㄡ"),
    ("xu", "ㄒㄩ"),
    ("xuan", "ㄒㄩㄢ"),
    ("xue", "ㄒㄩㄝ"),
    ("xun", "ㄒㄩㄣ"),
    ("ya", "ㄧㄚ"),
    ("yai", "ㄧㄞ"),
    ("yan", "ㄧㄢ"),
    ("yang", "ㄧㄤ"),
    ("yao", "ㄧㄠ"),
    ("ye", "ㄧㄝ"),
    ("yi", "ㄧ"),
    ("yin", "ㄧㄣ"),
    ("ying", "ㄧㄥ"),
    ("yo", "ㄧㄛ"),
    ("yong", "ㄩㄥ"),
    ("you", "ㄧㄡ"),
    ("yu", "ㄩ"),
    ("yuan", "ㄩㄢ"),
    ("yue", "ㄩㄝ"),
    ("yun", "ㄩㄣ"),
    ("za", "ㄗㄚ"),
    ("zai", "ㄗㄞ"),
    ("zan", "ㄗㄢ"),
    ("zang", "ㄗㄤ"),
    ("zao", "ㄗㄠ"),
    ("ze", "ㄗㄜ"),
    ("zei", "ㄗㄟ"),
    ("zen", "ㄗㄣ"),
    ("zeng", "ㄗㄥ"),
    ("zha", "ㄓㄚ"),
    ("zhai", "ㄓㄞ"),
    ("zhan", "ㄓㄢ"),
    ("zhang", "ㄓㄤ"),
    ("zhao", "ㄓㄠ"),
    ("zhe", "ㄓㄜ"),
    ("zhei", "ㄓㄟ"),
    ("zhen", "ㄓㄣ"),
    ("zheng", "ㄓㄥ"),
    ("zhi", "ㄓ"),
    ("zhong", "ㄓㄨㄥ"),
    ("zhou", "ㄓㄡ"),
    ("zhu", "ㄓㄨ"),
    ("zhua", "ㄓㄨㄚ"),
    ("zhuai", "ㄓㄨㄞ"),
    ("zhuan", "ㄓㄨㄢ"),
    ("zhuang", "ㄓㄨㄤ"),
    ("zhui", "ㄓㄨㄟ"),
    ("zhun", "ㄓㄨㄣ"),
    ("zhuo", "ㄓㄨㄛ"),
    ("zi", "ㄗ"),
    ("zong", "ㄗㄨㄥ"),
    ("zou", "ㄗㄡ"),
    ("zu", "ㄗㄨ"),
    ("zuan", "ㄗㄨㄢ"),
    ("zui", "ㄗㄨㄟ"),
    ("zun", "ㄗㄨㄣ"),
    ("zuo", "ㄗㄨㄛ")
])

const pinyin_trie = Trie([
    ("ㄚ", "a"),
    ("ㄞ", "ai"),
    ("ㄢ", "an"),
    ("ㄤ", "ang"),
    ("ㄠ", "ao"),
    ("ㄅㄚ", "ba"),
    ("ㄅㄞ", "bai"),
    ("ㄅㄢ", "ban"),
    ("ㄅㄤ", "bang"),
    ("ㄅㄠ", "bao"),
    ("ㄅㄟ", "bei"),
    ("ㄅㄣ", "ben"),
    ("ㄅㄥ", "beng"),
    ("ㄅㄧ", "bi"),
    ("ㄅㄧㄢ", "bian"),
    ("ㄅㄧㄠ", "biao"),
    ("ㄅㄧㄝ", "bie"),
    ("ㄅㄧㄣ", "bin"),
    ("ㄅㄧㄥ", "bing"),
    ("ㄅㄛ", "bo"),
    ("ㄅㄨ", "bu"),
    ("ㄘㄚ", "ca"),
    ("ㄘㄞ", "cai"),
    ("ㄘㄢ", "can"),
    ("ㄘㄤ", "cang"),
    ("ㄘㄠ", "cao"),
    ("ㄘㄜ", "ce"),
    ("ㄘㄣ", "cen"),
    ("ㄘㄥ", "ceng"),
    ("ㄔㄚ", "cha"),
    ("ㄔㄞ", "chai"),
    ("ㄔㄢ", "chan"),
    ("ㄔㄤ", "chang"),
    ("ㄔㄠ", "chao"),
    ("ㄔㄜ", "che"),
    ("ㄔㄣ", "chen"),
    ("ㄔㄥ", "cheng"),
    ("ㄔ", "chi"),
    ("ㄔㄨㄥ", "chong"),
    ("ㄔㄡ", "chou"),
    ("ㄔㄨ", "chu"),
    ("ㄔㄨㄚ", "chua"),
    ("ㄔㄨㄞ", "chuai"),
    ("ㄔㄨㄢ", "chuan"),
    ("ㄔㄨㄤ", "chuang"),
    ("ㄔㄨㄟ", "chui"),
    ("ㄔㄨㄣ", "chun"),
    ("ㄔㄨㄛ", "chuo"),
    ("ㄘ", "ci"),
    ("ㄘㄨㄥ", "cong"),
    ("ㄘㄡ", "cou"),
    ("ㄘㄨ", "cu"),
    ("ㄘㄨㄢ", "cuan"),
    ("ㄘㄨㄟ", "cui"),
    ("ㄘㄨㄣ", "cun"),
    ("ㄘㄨㄛ", "cuo"),
    ("ㄉㄚ", "da"),
    ("ㄉㄞ", "dai"),
    ("ㄉㄢ", "dan"),
    ("ㄉㄤ", "dang"),
    ("ㄉㄠ", "dao"),
    ("ㄉㄜ", "de"),
    ("ㄉㄥ", "deng"),
    ("ㄉㄧ", "di"),
    ("ㄉㄧㄢ", "dian"),
    ("ㄉㄧㄤ", "diang"),
    ("ㄉㄧㄠ", "diao"),
    ("ㄉㄧㄝ", "die"),
    ("ㄉㄧㄥ", "ding"),
    ("ㄉㄧㄡ", "diu"),
    ("ㄉㄨㄥ", "dong"),
    ("ㄉㄡ", "dou"),
    ("ㄉㄨ", "du"),
    ("ㄉㄨㄢ", "duan"),
    ("ㄉㄨㄟ", "dui"),
    ("ㄉㄨㄣ", "dun"),
    ("ㄉㄨㄛ", "duo"),
    ("ㄜ", "e"),
    ("ㄟ", "ei"),
    ("ㄣ", "en"),
    ("ㄦ", "er"),
    ("ㄈㄚ", "fa"),
    ("ㄈㄢ", "fan"),
    ("ㄈㄤ", "fang"),
    ("ㄈㄟ", "fei"),
    ("ㄈㄣ", "fen"),
    ("ㄈㄥ", "feng"),
    ("ㄈㄛ", "fo"),
    ("ㄈㄡ", "fou"),
    ("ㄈㄨ", "fu"),
    ("ㄍㄚ", "ga"),
    ("ㄍㄞ", "gai"),
    ("ㄍㄢ", "gan"),
    ("ㄍㄤ", "gang"),
    ("ㄍㄠ", "gao"),
    ("ㄍㄜ", "ge"),
    ("ㄍㄟ", "gei"),
    ("ㄍㄣ", "gen"),
    ("ㄍㄥ", "geng"),
    ("ㄍㄨㄥ", "gong"),
    ("ㄍㄡ", "gou"),
    ("ㄍㄨ", "gu"),
    ("ㄍㄨㄚ", "gua"),
    ("ㄍㄨㄞ", "guai"),
    ("ㄍㄨㄢ", "guan"),
    ("ㄍㄨㄤ", "guang"),
    ("ㄍㄨㄟ", "gui"),
    ("ㄍㄨㄣ", "gun"),
    ("ㄍㄨㄛ", "guo"),
    ("ㄏㄚ", "ha"),
    ("ㄏㄞ", "hai"),
    ("ㄏㄢ", "han"),
    ("ㄏㄤ", "hang"),
    ("ㄏㄠ", "hao"),
    ("ㄏㄜ", "he"),
    ("ㄏㄟ", "hei"),
    ("ㄏㄣ", "hen"),
    ("ㄏㄥ", "heng"),
    ("ㄏㄨㄥ", "hong"),
    ("ㄏㄡ", "hou"),
    ("ㄏㄨ", "hu"),
    ("ㄏㄨㄚ", "hua"),
    ("ㄏㄨㄞ", "huai"),
    ("ㄏㄨㄢ", "huan"),
    ("ㄏㄨㄤ", "huang"),
    ("ㄏㄨㄟ", "hui"),
    ("ㄏㄨㄣ", "hun"),
    ("ㄏㄨㄛ", "huo"),
    ("ㄐㄧ", "ji"),
    ("ㄐㄧㄚ", "jia"),
    ("ㄐㄧㄢ", "jian"),
    ("ㄐㄧㄤ", "jiang"),
    ("ㄐㄧㄠ", "jiao"),
    ("ㄐㄧㄝ", "jie"),
    ("ㄐㄧㄣ", "jin"),
    ("ㄐㄧㄥ", "jing"),
    ("ㄐㄩㄥ", "jiong"),
    ("ㄐㄧㄡ", "jiu"),
    ("ㄐㄩ", "ju"),
    ("ㄐㄩㄢ", "juan"),
    ("ㄐㄩㄝ", "jue"),
    ("ㄐㄩㄣ", "jun"),
    ("ㄎㄚ", "ka"),
    ("ㄎㄞ", "kai"),
    ("ㄎㄢ", "kan"),
    ("ㄎㄤ", "kang"),
    ("ㄎㄠ", "kao"),
    ("ㄎㄜ", "ke"),
    ("ㄎㄣ", "ken"),
    ("ㄎㄥ", "keng"),
    ("ㄎㄨㄥ", "kong"),
    ("ㄎㄡ", "kou"),
    ("ㄎㄨ", "ku"),
    ("ㄎㄨㄚ", "kua"),
    ("ㄎㄨㄞ", "kuai"),
    ("ㄎㄨㄢ", "kuan"),
    ("ㄎㄨㄤ", "kuang"),
    ("ㄎㄨㄟ", "kui"),
    ("ㄎㄨㄣ", "kun"),
    ("ㄎㄨㄛ", "kuo"),
    ("ㄌㄚ", "la"),
    ("ㄌㄞ", "lai"),
    ("ㄌㄢ", "lan"),
    ("ㄌㄤ", "lang"),
    ("ㄌㄠ", "lao"),
    ("ㄌㄜ", "le"),
    ("ㄌㄟ", "lei"),
    ("ㄌㄥ", "leng"),
    ("ㄌㄧ", "li"),
    ("ㄌㄧㄢ", "lian"),
    ("ㄌㄧㄤ", "liang"),
    ("ㄌㄧㄠ", "liao"),
    ("ㄌㄧㄝ", "lie"),
    ("ㄌㄧㄣ", "lin"),
    ("ㄌㄧㄥ", "ling"),
    ("ㄌㄧㄡ", "liu"),
    ("ㄌㄨㄥ", "long"),
    ("ㄌㄡ", "lou"),
    ("ㄌㄨ", "lu"),
    ("ㄌㄩ", "lü"),
    ("ㄌㄨㄢ", "luan"),
    ("ㄌㄩㄝ", "lüe"),
    ("ㄌㄨㄣ", "lun"),
    ("ㄌㄨㄛ", "luo"),
    ("ㄇㄚ", "ma"),
    ("ㄇㄞ", "mai"),
    ("ㄇㄢ", "man"),
    ("ㄇㄤ", "mang"),
    ("ㄇㄠ", "mao"),
    ("ㄇㄟ", "mei"),
    ("ㄇㄣ", "men"),
    ("ㄇㄥ", "meng"),
    ("ㄇㄧ", "mi"),
    ("ㄇㄧㄢ", "mian"),
    ("ㄇㄧㄠ", "miao"),
    ("ㄇㄧㄝ", "mie"),
    ("ㄇㄧㄣ", "min"),
    ("ㄇㄧㄥ", "ming"),
    ("ㄇㄧㄡ", "miu"),
    ("ㄇㄨㄛ", "mo"),
    ("ㄇㄡ", "mou"),
    ("ㄇㄨ", "mu"),
    ("ㄋㄚ", "na"),
    ("ㄋㄞ", "nai"),
    ("ㄋㄢ", "nan"),
    ("ㄋㄤ", "nang"),
    ("ㄋㄠ", "nao"),
    ("ㄋㄛ", "ne"),
    ("ㄋㄟ", "nei"),
    ("ㄋㄣ", "nen"),
    ("ㄋㄥ", "neng"),
    ("ㄋㄧ", "ni"),
    ("ㄋㄧㄚ", "nia"),
    ("ㄋㄧㄢ", "nian"),
    ("ㄋㄧㄤ", "niang"),
    ("ㄋㄧㄠ", "niao"),
    ("ㄋㄧㄝ", "nie"),
    ("ㄋㄧㄣ", "nin"),
    ("ㄋㄧㄥ", "ning"),
    ("ㄋㄧㄡ", "niu"),
    ("ㄋㄨㄥ", "nong"),
    ("ㄋㄡ", "nou"),
    ("ㄋㄨ", "nu"),
    ("ㄋㄩ", "nü"),
    ("ㄋㄨㄢ", "nuan"),
    ("ㄋㄩㄝ", "nüe"),
    ("ㄋㄨㄛ", "nuo"),
    ("ㄋㄨㄣ", "nun"),
    ("ㄡ", "ou"),
    ("ㄆㄚ", "pa"),
    ("ㄆㄞ", "pai"),
    ("ㄆㄢ", "pan"),
    ("ㄆㄤ", "pang"),
    ("ㄆㄠ", "pao"),
    ("ㄆㄟ", "pei"),
    ("ㄆㄣ", "pen"),
    ("ㄆㄥ", "peng"),
    ("ㄆㄧ", "pi"),
    ("ㄆㄧㄢ", "pian"),
    ("ㄆㄧㄠ", "piao"),
    ("ㄆㄧㄝ", "pie"),
    ("ㄆㄧㄣ", "pin"),
    ("ㄆㄧㄥ", "ping"),
    ("ㄆㄨㄛ", "po"),
    ("ㄆㄡ", "pou"),
    ("ㄆㄨ", "pu"),
    ("ㄑㄧ", "qi"),
    ("ㄑㄧㄚ", "qia"),
    ("ㄑㄧㄢ", "qian"),
    ("ㄑㄧㄤ", "qiang"),
    ("ㄑㄧㄠ", "qiao"),
    ("ㄑㄧㄝ", "qie"),
    ("ㄑㄧㄣ", "qin"),
    ("ㄑㄧㄥ", "qing"),
    ("ㄑㄩㄥ", "qiong"),
    ("ㄑㄧㄡ", "qiu"),
    ("ㄑㄩ", "qu"),
    ("ㄑㄩㄢ", "quan"),
    ("ㄑㄩㄝ", "que"),
    ("ㄑㄩㄣ", "qun"),
    ("ㄖㄢ", "ran"),
    ("ㄖㄤ", "rang"),
    ("ㄖㄠ", "rao"),
    ("ㄖㄜ", "re"),
    ("ㄖㄣ", "ren"),
    ("ㄖㄥ", "reng"),
    ("ㄖ", "ri"),
    ("ㄖㄨㄥ", "rong"),
    ("ㄖㄡ", "rou"),
    ("ㄖㄨ", "ru"),
    ("ㄖㄨㄢ", "ruan"),
    ("ㄖㄨㄟ", "rui"),
    ("ㄖㄨㄣ", "run"),
    ("ㄖㄨㄛ", "ruo"),
    ("ㄙㄚ", "sa"),
    ("ㄙㄞ", "sai"),
    ("ㄙㄢ", "san"),
    ("ㄙㄤ", "sang"),
    ("ㄙㄠ", "sao"),
    ("ㄙㄜ", "se"),
    ("ㄙㄟ", "sei"),
    ("ㄙㄣ", "sen"),
    ("ㄙㄥ", "seng"),
    ("ㄕㄚ", "sha"),
    ("ㄕㄞ", "shai"),
    ("ㄕㄢ", "shan"),
    ("ㄕㄤ", "shang"),
    ("ㄕㄠ", "shao"),
    ("ㄕㄜ", "she"),
    ("ㄕㄟ", "shei"),
    ("ㄕㄣ", "shen"),
    ("ㄕㄥ", "sheng"),
    ("ㄕ", "shi"),
    ("ㄕㄡㄥ", "shong"),
    ("ㄕㄡ", "shou"),
    ("ㄕㄨ", "shu"),
    ("ㄕㄨㄚ", "shua"),
    ("ㄕㄨㄞ", "shuai"),
    ("ㄕㄨㄢ", "shuan"),
    ("ㄕㄨㄤ", "shuang"),
    ("ㄕㄨㄟ", "shui"),
    ("ㄕㄨㄣ", "shun"),
    ("ㄕㄨㄛ", "shuo"),
    ("ㄙ", "si"),
    ("ㄙㄨㄥ", "song"),
    ("ㄙㄡ", "sou"),
    ("ㄙㄨ", "su"),
    ("ㄙㄨㄢ", "suan"),
    ("ㄙㄨㄟ", "sui"),
    ("ㄙㄨㄣ", "sun"),
    ("ㄙㄨㄛ", "suo"),
    ("ㄊㄚ", "ta"),
    ("ㄊㄞ", "tai"),
    ("ㄊㄢ", "tan"),
    ("ㄊㄤ", "tang"),
    ("ㄊㄠ", "tao"),
    ("ㄊㄜ", "te"),
    ("ㄊㄥ", "teng"),
    ("ㄊㄧ", "ti"),
    ("ㄊㄧㄢ", "tian"),
    ("ㄊㄧㄠ", "tiao"),
    ("ㄊㄧㄝ", "tie"),
    ("ㄊㄧㄥ", "ting"),
    ("ㄊㄨㄥ", "tong"),
    ("ㄊㄡ", "tou"),
    ("ㄊㄨ", "tu"),
    ("ㄊㄨㄢ", "tuan"),
    ("ㄊㄨㄟ", "tui"),
    ("ㄊㄨㄣ", "tun"),
    ("ㄊㄨㄛ", "tuo"),
    ("ㄨㄚ", "wa"),
    ("ㄨㄞ", "wai"),
    ("ㄨㄢ", "wan"),
    ("ㄨㄤ", "wang"),
    ("ㄨㄟ", "wei"),
    ("ㄨㄣ", "wen"),
    ("ㄨㄥ", "weng"),
    ("ㄨㄛ", "wo"),
    ("ㄨ", "wu"),
    ("ㄒㄧ", "xi"),
    ("ㄒㄧㄚ", "xia"),
    ("ㄒㄧㄢ", "xian"),
    ("ㄒㄧㄤ", "xiang"),
    ("ㄒㄧㄠ", "xiao"),
    ("ㄒㄧㄝ", "xie"),
    ("ㄒㄧㄣ", "xin"),
    ("ㄒㄧㄥ", "xing"),
    ("ㄒㄩㄥ", "xiong"),
    ("ㄒㄧㄡ", "xiu"),
    ("ㄒㄩ", "xu"),
    ("ㄒㄩㄢ", "xuan"),
    ("ㄒㄩㄝ", "xue"),
    ("ㄒㄩㄣ", "xun"),
    ("ㄧㄚ", "ya"),
    ("ㄧㄞ", "yai"),
    ("ㄧㄢ", "yan"),
    ("ㄧㄤ", "yang"),
    ("ㄧㄠ", "yao"),
    ("ㄧㄝ", "ye"),
    ("ㄧ", "yi"),
    ("ㄧㄣ", "yin"),
    ("ㄧㄥ", "ying"),
    ("ㄧㄛ", "yo"),
    ("ㄩㄥ", "yong"),
    ("ㄧㄡ", "you"),
    ("ㄩ", "yu"),
    ("ㄩㄢ", "yuan"),
    ("ㄩㄝ", "yue"),
    ("ㄩㄣ", "yun"),
    ("ㄗㄚ", "za"),
    ("ㄗㄞ", "zai"),
    ("ㄗㄢ", "zan"),
    ("ㄗㄤ", "zang"),
    ("ㄗㄠ", "zao"),
    ("ㄗㄜ", "ze"),
    ("ㄗㄟ", "zei"),
    ("ㄗㄣ", "zen"),
    ("ㄗㄥ", "zeng"),
    ("ㄓㄚ", "zha"),
    ("ㄓㄞ", "zhai"),
    ("ㄓㄢ", "zhan"),
    ("ㄓㄤ", "zhang"),
    ("ㄓㄠ", "zhao"),
    ("ㄓㄜ", "zhe"),
    ("ㄓㄟ", "zhei"),
    ("ㄓㄣ", "zhen"),
    ("ㄓㄥ", "zheng"),
    ("ㄓ", "zhi"),
    ("ㄓㄨㄥ", "zhong"),
    ("ㄓㄡ", "zhou"),
    ("ㄓㄨ", "zhu"),
    ("ㄓㄨㄚ", "zhua"),
    ("ㄓㄨㄞ", "zhuai"),
    ("ㄓㄨㄢ", "zhuan"),
    ("ㄓㄨㄤ", "zhuang"),
    ("ㄓㄨㄟ", "zhui"),
    ("ㄓㄨㄣ", "zhun"),
    ("ㄓㄨㄛ", "zhuo"),
    ("ㄗ", "zi"),
    ("ㄗㄨㄥ", "zong"),
    ("ㄗㄡ", "zou"),
    ("ㄗㄨ", "zu"),
    ("ㄗㄨㄢ", "zuan"),
    ("ㄗㄨㄟ", "zui"),
    ("ㄗㄨㄣ", "zun"),
    ("ㄗㄨㄛ", "zuo")
])

ext: jl


text:
using Dates

# Test data is 3 days after a starting reference on Jan 1:
# The actual phi phase is [0, 1, 2, 3]

geolist = [Date(2019, 1, 1), Date(2019, 1, 3), Date(2019, 1, 6), Date(2019, 1, 10)]
timediffs = InsarTimeseries.day_diffs(geolist)
@test timediffs == [2, 3, 4]

intlist = [(geolist[i], geolist[j]) for i = 1:length(geolist)-1 for j = i+1:length(geolist)]

# julia> intlist
#  (2019-01-01, 2019-01-03)
#  (2019-01-01, 2019-01-06)
#  (2019-01-01, 2019-01-10)
#  (2019-01-03, 2019-01-06)
#  (2019-01-03, 2019-01-10)
#  (2019-01-06, 2019-01-10)

A = InsarTimeseries.build_A_matrix(geolist, intlist)
@test A == [
    1 0 0
    0 1 0
    0 0 1
    -1 1 0
    -1 0 1
    0 -1 1
]

B = InsarTimeseries.build_B_matrix(geolist, intlist)
@test B ≈ [
    2.0 0.0 0.0
    2.0 3.0 0.0
    2.0 3.0 4.0
    0.0 3.0 0.0
    0.0 3.0 4.0
    0.0 0.0 4.0
]

truth_phi = [0.0, 1.0, 2.0, 3.0]
truth_v = [1 / 2, 1 / 3, 1 / 4]
dphis = [1.0, 2.0, 3.0, 1.0, 2.0, 1.0]
# Since we add no noise here, we should perfectly recover the phi phase
@test isapprox(A \ dphis, truth_phi[2:end], rtol = 1e-7)

v_hat = B \ dphis
@test isapprox(v_hat, truth_v, rtol = 1e-7)

phi_hat = InsarTimeseries.integrate_velocities(v_hat, timediffs)
@test isapprox(phi_hat, truth_phi, rtol = 1e-7)

dpr = dphis .+ rand(eltype(dphis), size(dphis))
@test isapprox(
    A \ dpr,
    InsarTimeseries.integrate_velocities(B \ dpr, timediffs)[2:end],
    rtol = 1e-7,
)
# TODO: make the const. vel into a function to test it


text:
## Utility functions.

"Get results of variables in a vector."
function value_vec(vec_x::Union{Array{VariableRef,1}, VariableRef})
    if isa(vec_x, VariableRef)
        vec_value = value(vec_x)
    else
        vec_value = [value(vec_x[i]) for i = 1:length(vec_x)]
    end

    return vec_value
end

"Get results of dual variables in a vector."
function dual_vec(vec_cons)
    return [dual(vec_cons[i]) for i = 1:length(vec_cons)]
end

@enum ColRow begin
   row = 1
   col = 2
end

"Check if the array is a column vector."
function check_col_vec(vec::Array{Int64,2}, str_name::String)
    if size(vec)[Int(col)] != 1
        throw("$str_name is not a column vector")
    end
end

# """

# Check if the size of the matrix if right.

# # Arguments:
# -
# """
# function check_mat_match(
#     vec::Array{Float64,2},
#     array::Array{Int64,2},
#     whi::ColRow,
#     str_len::String,
#     str_array::String,
# )
#     if size(array)[Int(whi)] != len
#         throw("The $(str(whi)) of $str_array doesn't match the $str_len.")
#     end
# end

"Check if the length of the list is right."
function check_li_match(
        len::Float64, list::Array{Int64,1}, str_len::String, str_list::String,
        )
    if length(list_1) != length(list_2)
        throw("The length of list $str_list doesn't match the $str_len.")
    end
end


text:
export Line, PeriodicLine

## Periodic line

# angle is (false==0) and π (true==1)
# or ranges from (-1,1]
struct PeriodicLine{angle,T} <: PeriodicDomain{T}
    center::T
    L::Float64
    PeriodicLine{angle,T}(c,L) where {angle,T} = new{angle,T}(c,L)
    PeriodicLine{angle,T}(c) where {angle,T} = new{angle,T}(c,1.)
    PeriodicLine{angle,T}(d::PeriodicLine) where {angle,T} = new{angle,T}(d.center,d.L)
    PeriodicLine{angle,T}() where {angle,T} = new{angle,T}(0.,1.)
end

(::Type{PeriodicLine{a}})(c,L) where {a} = PeriodicLine{a,typeof(c)}(c,L)

PeriodicLine(c,a) = PeriodicLine{a/π,eltype(c)}(c,1.)
PeriodicLine() = PeriodicLine{false,Float64}(0.,1.)
PeriodicLine(b::Bool) = PeriodicLine{b,Float64}()

isambiguous(d::PeriodicLine) = isnan(d.center) && isnan(d.angle)
convert(::Type{Domain{T}},d::PeriodicLine{a}) where {a,T<:Number} = PeriodicLine{a,T}(d.center, d.L)
convert(::Type{PeriodicLine{T,TT}},::AnyDomain) where {T<:Number,TT} = PeriodicLine{T,TT}(NaN,NaN)
convert(::Type{IT},::AnyDomain) where {IT<:PeriodicLine} = PeriodicLine(NaN,NaN)

angle(d::PeriodicLine{a}) where {a}=a*π

reverseorientation(d::PeriodicLine{true})=PeriodicLine{false}(d.center,d.L)
reverseorientation(d::PeriodicLine{false})=PeriodicLine{true}(d.center,d.L)
reverseorientation(d::PeriodicLine{a}) where {a}=PeriodicLine{a-1}(d.center,d.L)

tocanonical(d::PeriodicLine{false},x) = real(2atan((x-d.center)/d.L))
fromcanonical(d::PeriodicLine{false},θ) = d.L*tan(θ/2) + d.center

tocanonical(d::PeriodicLine{a},x) where {a} = tocanonical(PeriodicLine{false,Float64}(0.,d.L),exp(-π*im*a)*(x-d.center))
fromcanonical(d::PeriodicLine{a},v::AbstractArray) where {a} =
    [fromcanonical(d,vk) for vk in v]
fromcanonical(d::PeriodicLine{a},x) where {a} =
    exp(π*im*a)*fromcanonical(PeriodicLine{false,Float64}(0.,d.L),x)+d.center

function invfromcanonicalD(d::PeriodicLine{false})
    @assert d.center==0  && d.L==1.0
    a=Fun(PeriodicSegment(),[1.,0,1])
end

mappoint(a::PeriodicLine{false},b::Circle,x) = b.radius*((a.L*im-(x-a.center))./(a.L*im+(x-a.center)))+b.center
function mappoint(b::Circle,a::PeriodicLine{false},x)
    y=(x-b.center)./b.radius
    a.center+a.L*im*(1-y)./(y+1)
end

# algebra
*(c::Number,d::PeriodicLine)=PeriodicLine(isapprox(d.center,0) ? d.center : c*d.center,angle(d)+angle(c))
*(d::PeriodicLine,c::Number)=c*d
for OP in (:+,:-)
    @eval begin
        $OP(c::Number,d::PeriodicLine{a,T}) where {a,T}=PeriodicLine{a,promote_type(eltype(c),T)}($OP(c,d.center),d.L)
        $OP(d::PeriodicLine{a,T},c::Number) where {a,T}=PeriodicLine{a,promote_type(eltype(c),T)}($OP(d.center,c),d.L)
    end
end

arclength(d::PeriodicLine) = Inf
leftendpoint(d::PeriodicLine)= -Inf
rightendpoint(d::PeriodicLine)= Inf
complexlength(d::PeriodicLine)=Inf

## vectorized

function convert(::Type{PeriodicLine},d::ClosedInterval)
    a,b=d.left,d.right
    @assert abs(a) == abs(b) == Inf

    if isa(a,Real) && isa(b,Real)
        if a==Inf
            @assert b==-Inf
            PeriodicLine(true)
        else
            @assert a==-Inf&&b==Inf
            PeriodicLine(false)
        end
    elseif abs(real(a)) < Inf
        @assert real(a)==real(b)
        @assert sign(imag(a))==-sign(imag(b))

        PeriodicLine(real(b),angle(b))
    elseif isnan(real(a)) && isnan(real(b))  # hack for -im*Inf
        PeriodicLine([imag(a)*im,imag(b)*im])
    elseif isnan(real(a))  # hack for -im*Inf
        PeriodicLine([real(b)+imag(a)*im,b])
    elseif isnan(real(b))  # hack for -im*Inf
        PeriodicLine([a,real(a)+imag(b)*im])
    elseif abs(imag(a)) < Inf
        @assert imag(a)==imag(b)
        @assert sign(real(a))==-sign(real(b))

        PeriodicLine(imag(b),angle(b))
    else
        @assert angle(b) == -angle(a)

        PeriodicLine(0.,angle(b))
    end
end


text:
using Pkg, RData

cd("C:\\Users\\UOS\\Desktop\\portfolio\\dataset")

df_kospi = load("kospi_1018.RData")

typeof(df_kospi)

keys(df_kospi) ### Dict type의 data에 key를 확인

kospi_mat = df_kospi["df_kospi_1018"][2:end,]
group_inx = df_kospi["group_index"]

### Data type 변환 DF -> Matrix
kospi_mat = convert(Matrix, kospi_mat);

### log return으로  변환
log_return_mat= function(x)

    nrow = size(x)[1]
    ncol = size(x)[2]

    lag1 = x[2:end,:]

    log_return_mat = log.(lag1) .- log.(kospi_mat[1:(end-1),:])

    return log_return_mat
end

rmat = 100 * log_return_mat(kospi_mat)[1:240,:]

using Statistics

μ̂ = mean(rmat, dims = 1)

nrow = size(rmat)[1]
ncol = size(rmat)[2]

p = ncol; n = nrow

mean(μ̂)

using Random, Distributions, LinearAlgebra

μ₀ = 0.05

Amat = vcat(Diagonal(ones(p)), Diagonal(ones(p)), ones(p)', μ̂)
Bmat = vcat(Diagonal(-1*ones(2*p)), zeros(2*p)', zeros(2*p)')
Cmat = vcat(zeros(2*p), 1, μ₀)

X̃ = hcat(fill(-1, n), rmat)
K = hcat(zeros(p), Diagonal(ones(p)));

using StatsBase

ρ = 10; τ = 0.1

tmp_z = fill(0, 2*p)/p
tmp_u = fill(0, 2*p+2)
tmp_β = fill(1, p)/p
tmp_β₀ = 1e-5

tmp_β̃ = vcat(tmp_β₀, tmp_β);

pl = sqrt.(counts(group_inx)); iter_outer = 1000; iter_inner = 100;

λ₁ = 0.001; λ₂ = 0.001;

let i = 1

    while i <= iter_outer

        let j = 1

            while j <= iter_inner

                tmp_W = Diagonal(ifelse.(abs.((X̃ * tmp_β̃) * 4) .< 1e-8, 1e+8, abs.(1 ./ ((X̃ * tmp_β̃)*4))))

                tmp_β̃ = (-1/2) * inv(X̃' * tmp_W * X̃ + (ρ/2) * (Amat * K)' * (Amat * K)) * ((τ - (1/2)) * X̃' * ones(n) +
                    ρ * (Amat*K)' * ((1/ρ)*tmp_u + Bmat*tmp_z - Cmat))
                j += 1
            end
        end

        tmp_β = tmp_β̃[2:end]

        v = (Amat * tmp_β - Cmat)
        v₁ = v[1:p]
        tmp_z₁ = v₁ + tmp_u[1:p]/ρ

        ### Proximal operator
        tmp_z₁ = ifelse.(abs.(tmp_z₁) .<= λ₁, 0, tmp_z₁ - λ₁ * sign.(tmp_z₁))

        v₂ = v[(p+1):(2*p)]
        u₂ = tmp_u[(p+1):(2*p)]

        tmp_z₂ = v₂ + u₂/ρ

        for l in sort(unique(group_inx))
            plₗ = pl[l]
            v₂ₗ = v₂[group_inx .== l]
            u₂ₗ = u₂[group_inx .== l]

            if(norm(tmp_z₂[group_inx .== l], 2) .> λ₂ * plₗ)
                tmp_z₂[group_inx .== l] = tmp_z₂[group_inx .== l] - λ₂ * plₗ * tmp_z₂[group_inx .== l] / norm(tmp_z₂[group_inx .== l],2)
            else
                tmp_z₂[group_inx .== l] .= 0
            end
        end

        tmp_z = vcat(tmp_z₁, tmp_z₂)

        tmp_u = tmp_u + ρ * (Amat * tmp_β + Bmat * tmp_z - Cmat)

        i += 1

    end

    return tmp_β̃
end

sum(tmp_β̃[2:end])

μ̂ * tmp_β̃[2:end]

max_stars_repo_name: chulhongsung/JuliaProject

max_issues_repo_name: chulhongsung/JuliaProject

max_forks_repo_name: chulhongsung/JuliaProject


text:
using BoysFunction
using Test

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


text:
#Hello World Julia Script
println("Hello World")

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
function parse_line_from_stdin()
    map(x->parse(Int64,x), convert(Array{String,1},split(readline(STDIN), " ")))
end

N,num_query = parse_line_from_stdin()

Arrs = Array{Array{Int64,1},1}(N)
for i in 1:N
    Arrs[i] = Array{Int64,1}(0)
end

lastAns = 0
for i in 1:num_query
    query_type,x,y = parse_line_from_stdin()

    if query_type == 0
        ind = mod(x$lastAns,N)+1
        push!(Arrs[ind],y)
    elseif query_type == 1
        ind = mod(x$lastAns,N)+1
        seq = Arrs[ind]
        size = length(seq)
        val = seq[y%size+1]
        println(val)
        lastAns = val
    else
        error("unkown query type")
    end
end


text:
# https://docs.julialang.org/en/v1/base/io-network/#Base.IOBuffer

# String -> IO
io = IOBuffer()

# # TODO: Is there a length limit?
write(io, "This is a string buffer") # Write string to the buffer

x = take!(io)      # x is a byte array as shown below
typeof(x)          # Array{UInt8,1}
println(String(x)) # Convert to String


text:
module EasyFITSTests

using Test
using FITSIO
using EasyFITS
using EasyFITS: find

function samevalues(A::AbstractArray, B::AbstractArray)
    @assert !Base.has_offset_axes(A, B)
    @assert size(A) == size(B)
    for i ∈ eachindex(A, B)
        A[i] != B[i] && return false
    end
    return true
end

function maxabsdif(A::AbstractArray, B::AbstractArray)
    @assert !Base.has_offset_axes(A, B)
    @assert size(A) == size(B)
    T = promote_type(eltype(A), eltype(B))
    res = zero(T)
    for i ∈ eachindex(A, B)
        res = max(res, abs(A[i] - B[i]))
    end
    return res
end

path = joinpath(@__DIR__, "test.fits")

dat1 = rand(Float32, 3, 4, 5)
hdr1 = FitsHeader()

dat2 = rand(Int32, 7, 8)
hdr2 = FitsHeader(
    "HDUNAME" => ("HDU-TWO", "Second custom HDU"),
    "GA"      => (-42,       "Some integer keyword"),
    "BU"      => ("Gibi",    "Some string keyword"),
    "ZO"      => (sqrt(2),   "Some real keyword"),
    "MEU"     => (false,     "Some boolean keyword"))

dat3 = rand(Int32, 11)
img3 = FitsImage(dat3,
                 HDUNAME = ("HDU-THREE", "Third custom HDU"),
                 GA      = (-42,         "Some integer keyword"),
                 BU      = ("Gibi",      "Some string keyword"),
                 ZO      = (sqrt(2),     "Some real keyword"),
                 MEU     = (false,       "Some boolean keyword"))

compute_bitpix(::Type{T}) where {T<:Integer} = 8*sizeof(T)
compute_bitpix(::Type{T}) where {T<:AbstractFloat} = -8*sizeof(T)

@testset "AbstractArray interface" begin
    @test typeof(similar(img3)) === typeof(img3)
    let dims = (21, 8, 7), T = Float32
        @test typeof(similar(img3, T)) === FitsImage{T,ndims(img3)}
        @test typeof(similar(img3, dims...)) === FitsImage{eltype(img3),length(dims)}
        @test size(similar(img3, dims...)) === dims
        @test typeof(similar(img3, dims)) === FitsImage{eltype(img3),length(dims)}
        @test size(similar(img3, dims)) === dims
        @test typeof(similar(img3, T, dims...)) === FitsImage{T,length(dims)}
        @test size(similar(img3, T, dims...)) === dims
        @test typeof(similar(img3, T, dims)) === FitsImage{T,length(dims)}
        @test size(similar(img3, T, dims)) === dims
    end
    @test FitsBitpix(img3) === FitsBitpix(compute_bitpix(eltype(img3)))
    @test FitsBitpix(img3) === FitsBitpix(eltype(img3))
    @test FitsBitpix(img3) === FitsBitpix(typeof(img3))
    @test eltype(FitsBitpix(img3)) === eltype(img3)
end

@testset "Create FITS File" begin
    @test length(hdr1) == 0
    hdr1["HDUNAME"] = ("HDU-ONE", "First custom HDU")
    hdr1["GA"]      = (42,        "Some integer keyword")
    hdr1["BU"]      = ("Shadok",  "Some string keyword")
    hdr1["ZO"]      = (3.1415,    "Some real keyword")
    hdr1["MEU"]     = (true,      "Some boolean keyword")

    #FitsIO(path, "w!") do io
    #    write(io, dat1, hdr1)
    #    write(io, dat2, hdr2)
    #    write(io, img3)
    #end
    write(FitsFile, path, (dat1, hdr1), (dat2, hdr2), img3; overwrite=true)

    # Check overwrite is forbidden.
    @test_throws ErrorException close(FitsIO(path, "w"))
end

hduname(obj::Union{FitsHDU,FitsIO}) =
    get(String, obj, "HDUNAME", nothing)

@testset "Low-level" begin
    str = "\'\'\' hello \'\'you   \'"
    @test tryparse(String, EasyFITS.FitsUnparsedValue(str)) == "\' hello \'you"
    @test tryparse(String, EasyFITS.FitsUnparsedValue(str[2:end])) === nothing
    @test tryparse(String, EasyFITS.FitsUnparsedValue(str[1:end-1])) === nothing
    FitsIO(path) do io
        @test find(hdu -> hduname(hdu) == "HDU-MISSING", io) === nothing
        @test findfirst(hdu -> hduname(hdu) == "HDU-MISSING", io) === nothing
        @test findlast(hdu -> hduname(hdu) == "HDU-MISSING", io) === nothing
        @test findnext(hdu -> hduname(hdu) == "HDU-MISSING", io, 2) === nothing
        @test findprev(hdu -> hduname(hdu) == "HDU-MISSING", io, 2) === nothing
        @test find(hdu -> hduname(hdu) == "HDU-ONE", io) == 1
        @test find(hdu -> hduname(hdu) == "HDU-TWO", io) == 2
        @test findfirst(hdu -> hduname(hdu) == "HDU-ONE", io) == 1
        @test findlast(hdu -> hduname(hdu) == "HDU-ONE", io) == 1
        @test findnext(hdu -> hduname(hdu) == "HDU-ONE", io, 2) === nothing
        @test findprev(hdu -> hduname(hdu) == "HDU-ONE", io, 2) === 1
        hdu = io[1]
        @test get(Int,     hdu, "GA") == 42
        @test get(String,  hdu, "BU") == "Shadok"
        @test get(Float64, hdu, "ZO") ≈ 3.1415
        @test get(Bool,    hdu, "MEU") == true
        @test FitsBitpix(hdu) === FitsBitpix(dat1)
        @test FitsBitpix(io[2]) === FitsBitpix(dat2)
        @test FitsBitpix(io[3]) === FitsBitpix(img3)
        @test eltype(hdu) === eltype(dat1)
        @test eltype(io[2]) === eltype(dat2)
        @test eltype(io[3]) === eltype(dat3)
        dat = read(hdu)
        @test eltype(dat) == eltype(dat1)
        @test size(dat) == size(dat1)
        @test samevalues(dat, dat1)
        hdu = io[2]
        @test get(Int,     hdu, "GA") == -42
        @test get(String,  hdu, "BU") == "Gibi"
        @test get(Float64, hdu, "ZO") ≈ sqrt(2)
        @test get(Bool,    hdu, "MEU") == false
        dat = read(hdu)
        @test eltype(dat) == eltype(dat2)
        @test size(dat) == size(dat2)
        @test samevalues(dat, dat2)
    end
end

@testset "High-level" begin
    # Read array+header data and check indexing of array and header.
    A1 = read(FitsImage, path)
    @test isa(A1, FitsImage)
    @test isa(A1["GA"],  Int)     && A1["GA"] == 42
    @test isa(A1["BU"],  String)  && A1["BU"] == "Shadok"
    @test isa(A1["ZO"],  Float64) && A1["ZO"] ≈ 3.1415
    @test isa(A1["MEU"], Bool)    && A1["MEU"] == true
    @test A1.GA  == A1["GA"]
    @test A1.BU  == A1["BU"]
    @test A1.ZO  == A1["ZO"]
    @test A1.MEU == A1["MEU"]
    A1["ZO"] = 15
    @test isa(A1["ZO"],  Int) && A1["ZO"] == 15
    @test A1.ZO  == A1["ZO"]
    A1.ZO = 172
    @test isa(A1["ZO"],  Int) && A1["ZO"] == 172
    @test A1.ZO  == A1["ZO"]
    @test get(FitsComment, A1, "ZO") == "Some real keyword"
    #
    A1["ZO"] = (21, "New comment")
    @test isa(A1["ZO"],  Int) && A1["ZO"] == 21
    @test get(FitsComment, A1, "ZO") == "New comment"
    #
    A1["GA"] = 127, "Another comment"
    @test A1["GA"] == 127 && get(FitsComment, A1, "GA") == "Another comment"
    A1["GA"] = 128, nothing
    @test A1["GA"] == 128 && get(FitsComment, A1, "GA") == ""
    A1["GA"] = (129, )
    @test A1["GA"] == 129 && get(FitsComment, A1, "GA") == ""
    #
    A1.ZO = 220.0, "Yet another comment"
    @test A1.ZO == 220 && get(FitsComment, A1, "ZO") == "Yet another comment"
    A1.ZO = 221.0, nothing
    @test A1.ZO == 221 && get(FitsComment, A1, "ZO") == ""
    A1.ZO = (223.0, )
    @test A1.ZO == 223 && get(FitsComment, A1, "ZO") == ""
    A1.DUMMY1 = (1, "Dummy keyword")
    A1.DUMMY2 = (2, "Other dummy keyword")
    A1.DUMMY3 = (3, "Yet another dummy keyword")
    @test haskey(A1, "DUMMY1") == true
    delete!(A1, "DUMMY1")
    @test haskey(A1, "DUMMY1") == false
    @test haskey(A1, "DUMMY2") == true
    @test pop!(A1, "DUMMY2") === 2
    @test haskey(A1, "DUMMY2") == false
    #
    @test eltype(A1) == eltype(dat1)
    @test ndims(A1) == ndims(dat1)
    @test size(A1) == size(dat1)
    @test samevalues(A1, dat1)
    A2 = read(FitsImage, path, ext=2)
    @test isa(A2, FitsImage)
    @test isa(A2["GA"],  Int)     && A2["GA"] == -42
    @test isa(A2["BU"],  String)  && A2["BU"] == "Gibi"
    @test isa(A2["ZO"],  Float64) && A2["ZO"] ≈ sqrt(2)
    @test isa(A2["MEU"], Bool)    && A2["MEU"] == false
    @test eltype(A2) == eltype(dat2)
    @test ndims(A2) == ndims(dat2)
    @test size(A2) == size(dat2)
    @test samevalues(A2, dat2)
    # Read headers.
    H1 = read(FitsHeader, path, ext=1)
    @test isa(H1, FitsHeader)
    @test isa(H1["GA"],  Int)     && H1["GA"] == 42
    @test isa(H1["BU"],  String)  && H1["BU"] == "Shadok"
    @test isa(H1["ZO"],  Float64) && H1["ZO"] ≈ 3.1415
    @test isa(H1["MEU"], Bool)    && H1["MEU"] == true
    @test H1.GA  == H1["GA"]
    @test H1.BU  == H1["BU"]
    @test H1.ZO  == H1["ZO"]
    @test H1.MEU == H1["MEU"]
    @test FitsBitpix(H1) === FitsBitpix(dat1)
    @test eltype(H1) === eltype(dat1)
    H1["ZO"] = 15
    @test isa(H1["ZO"],  Int) && H1["ZO"] == 15
    @test H1.ZO  == H1["ZO"]
    H1.ZO = 172
    @test isa(H1["ZO"],  Int) && H1["ZO"] == 172
    @test H1.ZO  == H1["ZO"]
    @test get(FitsComment, H1, "ZO") == "Some real keyword"
    #
    H2 = read(FitsHeader, path, ext=2)
    @test isa(H2, FitsHeader)
    @test FitsBitpix(H2) === FitsBitpix(dat2)
    @test eltype(H2) === eltype(dat2)
    # Read array data with contraints.
    A3 = read(FitsArray, path, ext=1)
    @test isa(A3, Array{eltype(dat1),ndims(dat1)})
    @test size(A3) == size(dat1)
    @test samevalues(A3, dat1)
    A4 = read(FitsArray{Int}, path, ext=2)
    @test isa(A4, Array{Int,ndims(dat2)})
    @test size(A4) == size(dat2)
    @test samevalues(A4, dat2)
    @test_throws MethodError read(FitsArray{Int,ndims(dat2)+1}, path, ext=2)
    A5 = read(FitsArray{Int,ndims(dat2)}, path, ext=2)
    @test isa(A5, Array{Int,ndims(dat2)})
    @test size(A5) == size(dat2)
    @test samevalues(A5, dat2)
    # Read array+header data with contraints.
    A6 = read(FitsImage, path, ext=1)
    @test isa(A6, FitsImage{eltype(dat1),ndims(dat1)})
    @test size(A6) == size(dat1)
    @test samevalues(A6, dat1)
    A7 = read(FitsImage{Float64}, path, ext=1)
    @test isa(A7, FitsImage{Float64,ndims(dat1)})
    @test size(A7) == size(dat1)
    @test samevalues(A7, dat1)
    @test_throws MethodError read(FitsImage{Int,ndims(dat2)+1}, path, ext=2)
    A8 = read(FitsImage{Int,ndims(dat2)}, path, ext=2)
    @test isa(A8, FitsImage{Int,ndims(dat2)})
    @test size(A8) == size(dat2)
    @test samevalues(A8, dat2)
    @test read(FitsArray, path, :, 2, 2:3, ext=1) ≈ dat1[:, 2, 2:3]
end

#rm(path)

end # module


text:
######################## Relevance Matching Pursuit ###########################
# equivalent to consecutive runs of forward and backward regression
# δ is the largest marginal increase/decrease in residual norm before
# rmp with outer iterations
function rmp(A::AbstractMatOrFac, b::AbstractVector, δ::Real,
                                    maxiter::Int = 1, x = spzeros(size(A, 2)))
    n = size(A, 1)
    xt = copy(x) # to keep track of changes in x
    P = StepwiseRegression(A, b, x.nzind)
    for i in 1:maxiter
        # foward stage
        for _ in 1:n
            forward_step!(P, x, 0, δ) || break # breaks if no change occured
        end
        !(xt ≈ x) || break # break outer loop if x is not changing anymore
        xt .= x
        # backward stage
        for _ in nnz(x):-1:1
            backward_step!(P, x, Inf, δ) || break # breaks if no change occured
        end
        !(xt ≈ x) || break # break outer loop if x is not changing anymore
        xt .= x
    end
    return x
end

# max_ε is the residual norm tolerance
# δ is the largest marginal increase/decrease in residual norm before
# either forward or backward algorithm terminates
# k is the desired sparsity of the solution
function rmp(A::AbstractMatOrFac, b::AbstractVector,
            k::Int, x = spzeros(size(A, 2)))
    n = size(A, 1)
    P = StepwiseRegression(A, b, x.nzind)
    for _ in 1:n
        forward_step!(P, x, 0, 0) || break
    end
    for _ in nnz(x):-1:k+1
        backward_step!(P, x, Inf, Inf) || break
    end
    return x
end

# see Adaptive Forward-Backward Greedy Algorithm for Sparse Learning with Linear Models
# terminates if algorithm can't decrease residual norm by more than δ
function foba(A::AbstractMatOrFac, b::AbstractVector, δ::Real, x = spzeros(size(A, 2));
                                                            isfast::Val = Val(true))
    n = size(A, 1)
    P = StepwiseRegression(A, b, x.nzind)
    for i in 1:n
        forward_step!(P, x, 0, δ) || break
        max_δ = sqrt(maximum(P.δ²)) # this is the change in residual norm of last forward step
        while backward_step!(P, x, Inf, max_δ/2, isfast) end # FoBa only takes backward steps if they increase the norm by half of decrease of forward algorithm
    end
    return x
end


text:
#=
xmlparser:
- Julia version: 1.0.3
- Author: rhdekker
- Date: 2019-01-09
=#

using LibExpat
using DataStructures

# variant graph types and structs
abstract type VariantGraphNode end
struct StartNode <: VariantGraphNode end
struct TextNode <: VariantGraphNode
    content::String
end

# xml types and structs
mutable struct XMLBlock
    tag::String
    content::String
    tail::String
end

# type to group xml_blocks together
# The groups are not really necessary... Since they only contain the array
# they are a bit redundant.
mutable struct Group
    blocks::Array{XMLBlock}
end

# traverse all nodes of the XML tree in depth first fashion
# without recursion
# Initialize the deque with the children of the root
# daarna halen we er steeds 1 op
# Mocht dat een element zijn
# dan zetten we de kinderen daarvan weer op de deck.
# NOTE: maybe a stack is better here?
# for now we create an array with the result
# implementing it as an iterator would be nicer
function create_an_array_of_the_xml_nodes(root)
    result = Union{AbstractString, ETree}[]
    deck = deque(Union{AbstractString, ETree})
    push!(deck, root)

    while !isempty(deck)
        a = pop!(deck)
        push!(result, a)
        if typeof(a) == ETree
            for element in reverse(a.elements)
                push!(deck, element)
            end
        end
    end
    return result
end

function convert_to_xml_blocks(all_nodes)
    xml_blocks::Array{XMLBlock} = []
    for node in all_nodes
        if typeof(node) == ETree
            xml_block = XMLBlock(node.name, "", "")
            push!(xml_blocks, xml_block)
        else
            # node is a String
            xml_block = last(xml_blocks)
            if xml_block.content == ""
                xml_block.content = node
            else
                xml_block.tail = node
            end
        end
    end
    return xml_blocks
end

function is_xml_block_textual_variation_or_not(xml_block)
    # TODO: this should really be a constant
    # TODO: maybe also better to make it a set, depending on the size
    variation_tags = ["del", "add"]
    return xml_block.tag in variation_tags
end

# we hebben een fuctie nodig die bij elk xml block ervoor zorgt
# dat er nodes worden aangemaakt.
# dat is dus per xml block 1 of 2 text nodes
function create_nodes_for_each_xml_block(xml_blocks::Array{XMLBlock})
    nodes::Array{TextNode} = []
    for block in xml_blocks
        !isempty(block.content) && push!(nodes, TextNode(block.content))
        !isempty(block.tail) && push!(nodes, TextNode(block.tail))
    end
    return nodes
end

function transform_block_into_text_nodes_named_tuple(block::XMLBlock)
    isempty(block.content) && return (tag=block.tag)
    isempty(block.tail) && return (tag=block.tag, content=TextNode(block.content))
    return (tag=block.tag, content=TextNode(block.content), tail=TextNode(block.tail))
end

# function transform_block_into_text_nodes(block::XMLBlock)::Array{TextNode}
#     nodes = []
#     !isempty(block.content) && push!(nodes, TextNode(block.content))
#     !isempty(block.tail) && push!(nodes, TextNode(block.tail))
#     return nodes
# end

# NOTE: this should probably be done with a fold
# Basically what happens here is that I turn an array into an array of arrays, also
# known as a tree :-).
# NOTE: There is a group by function in JuliaCollections / IterTools.jl
# but this function only looks at one element at the time
# there must be another function that groups elements into (previous, current) tuples.
# Or is it a take until problem?
function partition_block_into_groups(all_block)
    groups = []
    # we willen over alle nodes lopen steeds per twee, dus we houden een previous bij..
    previous = undef
    for block in all_block
        if previous == undef
            previous = block
            # create a new group for the first block
            group = Group([])
            # and add the block to the group
            push!(group.blocks, block)
            # add the group to the result
            push!(groups, group)
            continue
        end
        # hmm dit doe ik niet goed
        # compare an aspect of the current block with the previous block
        # if the aspect is different, but it in a new block
        # otherwise
        a = is_xml_block_textual_variation_or_not(previous)
        b = is_xml_block_textual_variation_or_not(block)

        # TODO: here we also need to check whether there is trialing text
        # if there is trailing text in a then they can not be placed in the same group
        if (a == b)
            # voeg toe aan de vorige group
            previous_group = last(groups)
            push!(previous_group.blocks, block)
        else
            group = Group([])
            push!(group.blocks, block)
            push!(groups, group)
        end
        previous = block
    end
    return groups
end

function main()
    xml = "<xml>Mondays are <del>well good</del><add>def bad</add>!</xml>"
    root = xp_parse(xml)
    # dump(root)

    all_nodes = create_an_array_of_the_xml_nodes(root)
    blocks = convert_to_xml_blocks(all_nodes)
    println(blocks)

    # We mappen de XML blocks naar een combo van tag en nul of meer text nodes.

    tuples_of_tag_and_text_nodes::Array{NamedTuple} = []
    for block in blocks
        named_tuple = transform_block_into_text_nodes_named_tuple(block)
        push!(tuples_of_tag_and_text_nodes, named_tuple)
    end

    println(tuples_of_tag_and_text_nodes)

#     textNodes = create_nodes_for_each_xml_block(blocks)
#     println(textNodes)

    # de partities doen we later wel.
    #partitions = partition_block_into_groups(blocks)
    #println(partitions)
end

main()


text:
using LogDensityProblems, CSV, Optim, DataFrames
using Distributions, Parameters
using Calculus, StatsPlots, PositiveFactorizations

include(string(pwd(),"/src/simulation.jl"))
include(string(pwd(),"/mcmc/priors_carlos.jl"))
include(string(pwd(),"/mcmc/foos_carlos.jl"))

dados = DataFrame(CSV.file("data/canada_data.csv"))

dados = Array(dados)[:,2:9]

par_cal = [0.67; 0.989; 1.66; # alfa, beta, phi
      #-0.34; -7.12; -8.98; # A0
       -40.66; -0.4467; -84.6557;
       0.872; 1.34e-3; 4.911e-05; -29.87; -0.212738; 0.002511; -84.1538; 0.31; 0.9410; # A1
       0.005; 0.798; 1.36; # shocks_star
       0.68; 1.27; 0.58; 0.07; 0.82; 0.13; 0.92; 0.02; 0.79; 0.60; 0.12; 0.88; 5.57; 10.34; 2.08; 1.10; 0.54; 0.30; 1.99; 0.70; 0.31; 0.79]

prob = DSGE_Model(dados,par_cal[1],par_cal[2],par_cal[3],par_cal[4:18])

t = problem_transform(prob)

P = TransformedLogDensity(t,prob)

unif = Uniform(0,1)

num_iter = 1_000_000

pars_aceitos = zeros(num_iter,22)

par_chute = randn(22)

nll(x) = -1*LogDensityProblems.logdensity(P,x)

mode = optimize(nll,par_chute, NelderMead(),Optim.Options(iterations =50_000))

start_par = mode.minimizer

#[rand(prior_bet),rand(prior_epsilon),rand(prior_theta),rand(prior_sig),rand(prior_s2),rand(prior_phi),rand(prior_phi_pi),rand(prior_phi_y),rand(prior_rho_v)]

npar = 22

coef_escala = 0.0075

hes = -Calculus.hessian(x->LogDensityProblems.logdensity(P,x),start_par)

hes_inv = pinv(hes)

hes_inv = (hes_inv + hes_inv')/2

hes_inv = cholesky(Positive,hes_inv)

hes_inv = hes_inv.L*hes_inv.L'

#hes = round.(hes;digits=3)

isposdef(hes_inv)

pars_aceitos[1,:] = pars_inicio

scale_matrix = cov(pars_aceitos[1:4000,:])

j = 2
rejec = 0

while j <= num_iter
    kernel_velho = MvNormal(pars_aceitos[j-1,:],coef_escala*scale_matrix)
    novo_par = rand(kernel_velho)
    #pars_aceitos[j,2:10] = novo_par
    kernel_novo = MvNormal(novo_par,coef_escala*scale_matrix)

    teste = LogDensityProblems.logdensity(P,novo_par)
    if isnan(teste)
        #pars_aceitos[j,2:10] = pars_aceitos[j-1,2:10]
        #global rejec += 1
        #global j += 1
        continue
    else
        num = teste + logpdf(kernel_novo,pars_aceitos[j-1,:])
        #println(string(teste," ",logpdf(kernel_novo, pars_aceitos[j-1,2:10])))
        dem = LogDensityProblems.logdensity(P,pars_aceitos[j-1,:]) + logpdf(kernel_velho,novo_par)
        alpha = exp(min(0,num - dem))
        p = rand(unif)
        if alpha < p
            pars_aceitos[j,:] = pars_aceitos[j-1,:]
            global rejec += 1
        else
            pars_aceitos[j,:] = novo_par
        end
        acc = 1 - rejec/j
        if j % 50 == 0
            println("Iteração ", j, " taxa de aceitação ", acc)
            sleep(0.4)
        end
        global j += 1
    end
end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using DrWatson, Test

# Test commit function
com = current_commit(dirname(@__DIR__))
@test com !== nothing
@test typeof(com) == String

# Test dictionary expansion
c = Dict(:a => [1, 2], :b => 4);
c1 = [ Dict(:a=>1,:b=>4)
 Dict(:a=>2,:b=>4)]

v1 = dict_list(c)
for el in c1
    @test el ∈ v1
end

c[:c] = "test"; c[:d] = ["lala", "lulu"];
c2 = [ Dict(:a=>1,:b=>4,:d=>"lala",:c=>"test")
 Dict(:a=>2,:b=>4,:d=>"lala",:c=>"test")
 Dict(:a=>1,:b=>4,:d=>"lulu",:c=>"test")
 Dict(:a=>2,:b=>4,:d=>"lulu",:c=>"test")]

v2 = dict_list(c)
for el in c2
    @test el ∈ v2
end

c[:e] = [[1, 2], [3, 5]];
c3 = [
Dict(:a=>1,:b=>4,:d=>"lala",:e=>[1, 2],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lala",:e=>[1, 2],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lulu",:e=>[1, 2],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lulu",:e=>[1, 2],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lala",:e=>[3, 5],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lala",:e=>[3, 5],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lulu",:e=>[3, 5],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lulu",:e=>[3, 5],:c=>"test")
]

v3 = dict_list(c)
@test dict_list_count(c) == length(c3) == length(v3)
for el in c3
    @test el ∈ v3
end

# tag!
d1 = Dict(:x => 3, :y => 4)
d2 = Dict("x" => 3, "y" => 4)
for d in (d1, d2)
    d = tag!(d, dirname(@__DIR__))

    @test haskey(d, keytype(d)(:commit))
    @test d[keytype(d)(:commit)] |> typeof == String
end


text:
using memo
using Test

@testset "memo" begin
    # Write your tests here.
end


text:
Test.@testset "assert.jl" begin
    Test.@test_nowarn CompatHelper.always_assert(true)
    Test.@test CompatHelper.always_assert(true) isa Nothing
    Test.@test CompatHelper.always_assert(true) == nothing
    Test.@test Test.@test_nowarn CompatHelper.always_assert(true) isa Nothing
    Test.@test Test.@test_nowarn CompatHelper.always_assert(true) == nothing
    Test.@test_throws CompatHelper.AlwaysAssertionError CompatHelper.always_assert(false)
end

Test.@testset "ci_service.jl" begin
    withenv("GITHUB_REPOSITORY" => "foo/bar") do
        Test.@test CompatHelper.auto_detect_ci_service() isa CompatHelper.CIService
        Test.@test CompatHelper.auto_detect_ci_service() isa CompatHelper.GitHubActions
        Test.@test_throws ArgumentError CompatHelper.main(; keep_existing_compat = false, drop_existing_compat = false)
    end
    withenv("GITHUB_REPOSITORY" => nothing) do
        Test.@test_throws ErrorException CompatHelper.auto_detect_ci_service()
    end
end

Test.@testset "git.jl" begin
    Test.@test CompatHelper.git_decide_master_branch(CompatHelper.DefaultBranch(), "  abcdefg  ") == "abcdefg"
    Test.@test CompatHelper.git_decide_master_branch(" abc ", " xyz ") == "abc"
end

Test.@testset "new_versions.jl" begin
    Test.@test_throws ArgumentError CompatHelper.old_compat_to_new_compat("", "", :abc)
end


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 75225d8a-ad5e-11eb-131a-b76726982c01
using CSV, DataFrames, Query, Statistics, NaNMath

# ╔═╡ 91525f7c-027d-41b2-a2eb-33f9e1eddf46
using Plots

# ╔═╡ 630e63c6-d61e-46ad-b895-05db5902ad27
begin
	data = Dict{String,Any}()
	path = "/mnt/data_simons/NEID/DRPv0.7-fixedflatfielding2/output/manifests"
	for (root, dirs, files) in walkdir(path)
	    #=
		println("Directories in $root")
		for dir in dirs
	        println(joinpath(root, dir)) # path to directories
	    end
		=#
	    println("Files in $root")
	    for file in files
			m = match(r"(\d+)/manifest\.csv",joinpath(root, file))
			if m == nothing continue end
			datestr = m.captures[1]
	        #println(joinpath(root, file)) # path to files
			println(datestr) # path to files
			df = CSV.read(joinpath(root, file),DataFrame)
			data[datestr] = df
	    end
		#break
	end
end

# ╔═╡ faf36f02-c231-4cb0-a095-d981762a2f99
begin
	order_to_use = 60
	for (i,k) in enumerate(sort(collect(keys(data))))
		nobs = length(data[k].order_snrs)
		if nobs == 0
			pop!(data,k)
		end
	end
	ndays = length(keys(data))
	days_snrs = fill(Vector{Float64}(undef,0),ndays)
	days_airmass = fill(Vector{Float64}(undef,0),ndays)
	days_sha = fill(Vector{Float64}(undef,0),ndays)
	for (i,k) in enumerate(sort(collect(keys(data))))
		nobs = length(data[k].order_snrs)
		days_snrs[i] = map(i->parse.(Float64,split(data[k].order_snrs[i][2:end-1],','))[order_to_use],1:length(data[k].order_snrs))
		days_airmass[i] = map(i->data[k].airmass[i],1:length(data[k].airmass))
		days_sha[i] = map(i->data[k].bjd[i].-floor(data[k].bjd[i]),1:length(data[k].bjd))
	end
end

# ╔═╡ f6545aab-0865-44aa-973a-644101b65e20
scatter(mean.(days_snrs))

# ╔═╡ 10da0778-c51c-4f98-90e6-b1d647274b98
scatter(maximum.(days_snrs))

# ╔═╡ aa3aa153-51e1-4649-af89-07246ed0c00d
scatter(std.(days_snrs))

# ╔═╡ 62e42b0f-07e9-47b7-bc2c-d52c86c65e38
threshold = maximum(maximum.(days_snrs))

# ╔═╡ 755fbea9-adcf-4c24-b837-118aaba14ff6
begin
	plt = plot(legend=:none)
	for obsid in 1:ndays
		idx_keep = days_snrs[obsid] .>= 0.80*threshold
		scatter!(plt,days_sha[obsid][idx_keep].*24,days_snrs[obsid][idx_keep]./threshold, ms=1)
		idx_keep = days_snrs[obsid] .>= 0.92*threshold
		scatter!(plt,days_sha[obsid][idx_keep].*24,days_snrs[obsid][idx_keep]./threshold, ms=2)
	end
	plt
end

# ╔═╡ 9db1b522-f069-4319-8563-762e33db910f
4/24

# ╔═╡ Cell order:
# ╠═75225d8a-ad5e-11eb-131a-b76726982c01
# ╠═630e63c6-d61e-46ad-b895-05db5902ad27
# ╠═faf36f02-c231-4cb0-a095-d981762a2f99
# ╠═f6545aab-0865-44aa-973a-644101b65e20
# ╠═10da0778-c51c-4f98-90e6-b1d647274b98
# ╠═aa3aa153-51e1-4649-af89-07246ed0c00d
# ╠═91525f7c-027d-41b2-a2eb-33f9e1eddf46
# ╠═755fbea9-adcf-4c24-b837-118aaba14ff6
# ╠═62e42b0f-07e9-47b7-bc2c-d52c86c65e38
# ╠═9db1b522-f069-4319-8563-762e33db910f


text:
# This file contains functions that are related to the AMR capabilities of the DG solver

# Refine elements in the DG solver based on a list of cell_ids that should be refined
function refine!(dg::Dg3D{Eqn, MeshType, NVARS, POLYDEG}, mesh::TreeMesh,
                 cells_to_refine::AbstractArray{Int}) where {Eqn, MeshType, NVARS, POLYDEG}
  # Return early if there is nothing to do
  if isempty(cells_to_refine)
    return
  end

  # Determine for each existing element whether it needs to be refined
  needs_refinement = falses(nelements(dg.elements))
  tree = mesh.tree
  # The "Ref(...)" is such that we can vectorize the search but not the array that is searched
  elements_to_refine = searchsortedfirst.(Ref(dg.elements.cell_ids[1:nelements(dg.elements)]),
                                          cells_to_refine)
  needs_refinement[elements_to_refine] .= true

  # Retain current solution data
  old_n_elements = nelements(dg.elements)
  old_u = dg.elements.u

  # Get new list of leaf cells
  leaf_cell_ids = leaf_cells(tree)

  # Initialize new elements container
  elements = init_elements(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG))
  n_elements = nelements(elements)

  # Loop over all elements in old container and either copy them or refine them
  u_tmp1 = Array{Float64, 4}(undef, nvariables(dg), nnodes(dg), nnodes(dg), nnodes(dg))
  u_tmp2 = Array{Float64, 4}(undef, nvariables(dg), nnodes(dg), nnodes(dg), nnodes(dg))
  element_id = 1
  for old_element_id in 1:old_n_elements
    if needs_refinement[old_element_id]
      # Refine element and store solution directly in new data structure
      refine_element!(elements.u, element_id, old_u, old_element_id,
                      dg.mortar_forward_upper, dg.mortar_forward_lower, dg, u_tmp1, u_tmp2)
      element_id += 2^ndims(dg)
    else
      # Copy old element data to new element container
      @views elements.u[:, :, :, :, element_id] .= old_u[:, :, :, :, old_element_id]
      element_id += 1
    end
  end

  # Initialize new interfaces container
  interfaces = init_interfaces(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG), elements)
  n_interfaces = ninterfaces(interfaces)

  # Initialize boundaries
  boundaries, n_boundaries_per_direction = init_boundaries(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG), elements)
  n_boundaries = nboundaries(boundaries)

  # Initialize new mortar containers
  l2mortars = init_mortars(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG), elements, dg.mortar_type)
  n_l2mortars = nmortars(l2mortars)

  # Sanity check
  if isperiodic(mesh.tree) && n_l2mortars == 0
    @assert n_interfaces == 3*n_elements ("For 3D and periodic domains and conforming elements, "
                                        * "n_surf must be the same as 3*n_elem")
  end

  # Update DG instance with new data
  dg.elements = elements
  dg.n_elements = n_elements
  dg.n_elements_global = n_elements
  dg.interfaces = interfaces
  dg.n_interfaces = n_interfaces
  dg.boundaries = boundaries
  dg.n_boundaries = n_boundaries
  dg.n_boundaries_per_direction = n_boundaries_per_direction
  dg.l2mortars = l2mortars
  dg.n_l2mortars = n_l2mortars
end

# Refine solution data u for an element, using L2 projection (interpolation)
function refine_element!(u, element_id, old_u, old_element_id,
                         forward_upper, forward_lower, dg::Dg3D, u_tmp1, u_tmp2)
  # Store new element ids
  bottom_lower_left_id  = element_id
  bottom_lower_right_id = element_id + 1
  bottom_upper_left_id  = element_id + 2
  bottom_upper_right_id = element_id + 3
  top_lower_left_id     = element_id + 4
  top_lower_right_id    = element_id + 5
  top_upper_left_id     = element_id + 6
  top_upper_right_id    = element_id + 7

  # Interpolate to bottom lower left element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, bottom_lower_left_id), forward_lower, forward_lower, forward_lower,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to bottom lower right element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, bottom_lower_right_id), forward_upper, forward_lower, forward_lower,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to bottom upper left element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, bottom_upper_left_id), forward_lower, forward_upper, forward_lower,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to bottom upper right element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, bottom_upper_right_id), forward_upper, forward_upper, forward_lower,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to top lower left element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, top_lower_left_id), forward_lower, forward_lower, forward_upper,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to top lower right element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, top_lower_right_id), forward_upper, forward_lower, forward_upper,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to top upper left element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, top_upper_left_id), forward_lower, forward_upper, forward_upper,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)

  # Interpolate to top upper right element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, top_upper_right_id), forward_upper, forward_upper, forward_upper,
    view(old_u, :, :, :, :, old_element_id), u_tmp1, u_tmp2)
end

# Coarsen elements in the DG solver based on a list of cell_ids that should be removed
function coarsen!(dg::Dg3D{Eqn, MeshType, NVARS, POLYDEG}, mesh::TreeMesh,
                  child_cells_to_coarsen::AbstractArray{Int}) where {Eqn, MeshType, NVARS, POLYDEG}
  # Return early if there is nothing to do
  if isempty(child_cells_to_coarsen)
    return
  end

  # Determine for each old element whether it needs to be removed
  to_be_removed = falses(nelements(dg.elements))
  # The "Ref(...)" is such that we can vectorize the search but not the array that is searched
  elements_to_remove = searchsortedfirst.(Ref(dg.elements.cell_ids[1:nelements(dg.elements)]),
                                          child_cells_to_coarsen)
  to_be_removed[elements_to_remove] .= true

  # Retain current solution data
  old_n_elements = nelements(dg.elements)
  old_u = dg.elements.u

  # Get new list of leaf cells
  leaf_cell_ids = leaf_cells(mesh.tree)

  # Initialize new elements container
  elements = init_elements(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG))
  n_elements = nelements(elements)

  # Loop over all elements in old container and either copy them or coarsen them
  u_tmp1 = Array{Float64, 4}(undef, nvariables(dg), nnodes(dg), nnodes(dg), nnodes(dg))
  u_tmp2 = Array{Float64, 4}(undef, nvariables(dg), nnodes(dg), nnodes(dg), nnodes(dg))
  skip = 0
  element_id = 1
  for old_element_id in 1:old_n_elements
    # If skip is non-zero, we just coarsened 2^ndims elements and need to omit the following elements
    if skip > 0
      skip -= 1
      continue
    end

    if to_be_removed[old_element_id]
      # If an element is to be removed, sanity check if the following elements
      # are also marked - otherwise there would be an error in the way the
      # cells/elements are sorted
      @assert all(to_be_removed[old_element_id:(old_element_id+2^ndims(dg)-1)]) "bad cell/element order"

      # Coarsen elements and store solution directly in new data structure
      coarsen_elements!(elements.u, element_id, old_u, old_element_id,
                        dg.l2mortar_reverse_upper, dg.l2mortar_reverse_lower, dg, u_tmp1, u_tmp2)
      element_id += 1
      skip = 2^ndims(dg) - 1
    else
      # Copy old element data to new element container
      @views elements.u[:, :, :, :, element_id] .= old_u[:, :, :, :, old_element_id]
      element_id += 1
    end
  end

  # Initialize new interfaces container
  interfaces = init_interfaces(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG), elements)
  n_interfaces = ninterfaces(interfaces)

  # Initialize boundaries
  boundaries, n_boundaries_per_direction = init_boundaries(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG), elements)
  n_boundaries = nboundaries(boundaries)

  # Initialize new mortar containers
  l2mortars = init_mortars(leaf_cell_ids, mesh, Val(NVARS), Val(POLYDEG), elements, dg.mortar_type)
  n_l2mortars = nmortars(l2mortars)

  # Sanity check
  if isperiodic(mesh.tree) && n_l2mortars == 0
    @assert n_interfaces == 3*n_elements ("For 3D and periodic domains and conforming elements, "
                                        * "n_surf must be the same as 3*n_elem")
  end

  # Update DG instance with new data
  dg.elements = elements
  dg.n_elements = n_elements
  dg.n_elements_global = n_elements
  dg.interfaces = interfaces
  dg.n_interfaces = n_interfaces
  dg.boundaries = boundaries
  dg.n_boundaries = n_boundaries
  dg.n_boundaries_per_direction = n_boundaries_per_direction
  dg.l2mortars = l2mortars
  dg.n_l2mortars = n_l2mortars
end

# Coarsen solution data u for four elements, using L2 projection
function coarsen_elements!(u, element_id, old_u, old_element_id,
                           reverse_upper, reverse_lower, dg::Dg3D, u_tmp1, u_tmp2)
  # Store old element ids
  bottom_lower_left_id  = old_element_id
  bottom_lower_right_id = old_element_id + 1
  bottom_upper_left_id  = old_element_id + 2
  bottom_upper_right_id = old_element_id + 3
  top_lower_left_id     = old_element_id + 4
  top_lower_right_id    = old_element_id + 5
  top_upper_left_id     = old_element_id + 6
  top_upper_right_id    = old_element_id + 7

  # Project from bottom lower left element
  multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_lower, reverse_lower, reverse_lower,
    view(old_u, :, :, :, :, bottom_lower_left_id), u_tmp1, u_tmp2)

  # Project from bottom lower right element_variables
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_upper, reverse_lower, reverse_lower,
    view(old_u, :, :, :, :, bottom_lower_right_id), u_tmp1, u_tmp2)

  # Project from bottom upper left element
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_lower, reverse_upper, reverse_lower,
    view(old_u, :, :, :, :, bottom_upper_left_id), u_tmp1, u_tmp2)

  # Project from bottom upper right element
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_upper, reverse_upper, reverse_lower,
    view(old_u, :, :, :, :, bottom_upper_right_id), u_tmp1, u_tmp2)

  # Project from top lower left element
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_lower, reverse_lower, reverse_upper,
    view(old_u, :, :, :, :, top_lower_left_id), u_tmp1, u_tmp2)

  # Project from top lower right element
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_upper, reverse_lower, reverse_upper,
    view(old_u, :, :, :, :, top_lower_right_id), u_tmp1, u_tmp2)

  # Project from top upper left element
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_lower, reverse_upper, reverse_upper,
    view(old_u, :, :, :, :, top_upper_left_id), u_tmp1, u_tmp2)

  # Project from top upper right element
  add_multiply_dimensionwise!(
    view(u,     :, :, :, :, element_id), reverse_upper, reverse_upper, reverse_upper,
    view(old_u, :, :, :, :, top_upper_right_id), u_tmp1, u_tmp2)
end

# Calculate an AMR indicator value for each element/leaf cell
#
# The indicator value λ ∈ [-1,1] is ≈ -1 for cells that should be coarsened, ≈
# 0 for cells that should remain as-is, and ≈ 1 for cells that should be
# refined.
#
# Note: The implementation here implicitly assumes that we have an element for
# each leaf cell and that they are in the same order.
#
# FIXME: This is currently implemented for each test case - we need something
# appropriate that is both equation and test case independent
function calc_amr_indicator(dg::Dg3D, mesh::TreeMesh, time::Float64)
  lambda = zeros(dg.n_elements)

  if dg.amr_indicator === :gauss
    base_level = 4
    max_level = 6
    threshold_high = 0.6
    threshold_low = 0.1

    # Iterate over all elements
    for element_id in 1:dg.n_elements
      # Determine target level from peak value
      peak = maximum(dg.elements.u[:, :, :, :, element_id])
      if peak > threshold_high
        target_level = max_level
      elseif peak > threshold_low
        target_level = max_level - 1
      else
        target_level = base_level
      end

      # Compare target level with actual level to set indicator
      cell_id = dg.elements.cell_ids[element_id]
      actual_level = mesh.tree.levels[cell_id]
      if actual_level < target_level
        lambda[element_id] = 1.0
      elseif actual_level > target_level
        lambda[element_id] = -1.0
      else
        lambda[element_id] = 0.0
      end
    end
  elseif dg.amr_indicator === :blob
    base_level = 1
    max_level = 6 # originally 6; set to 7 (or even 8) to get increased resolution
    blending_factor_threshold1 = 0.3 # Löhner original choice
    blending_factor_threshold2 = 0.1 # Löhner original choice

    # (Re-)initialize element variable storage for blending factor
    if (!haskey(dg.element_variables, :amr_indicator_values) ||
        length(dg.element_variables[:amr_indicator_values]) != dg.n_elements)
      dg.element_variables[:amr_indicator_values] = Vector{Float64}(undef, dg.n_elements)
    end

    alpha = dg.element_variables[:amr_indicator_values]
    calc_loehner_indicator!(alpha, dg.elements.u, density, dg.thread_cache, dg)

    # OPTIONAL: use max level were shock capturing is maxed
    # (Re-)initialize element variable storage for blending factor
    #if (!haskey(dg.element_variables, :blending_factor) ||
    #    length(dg.element_variables[:blending_factor]) != dg.n_elements)
    #  dg.element_variables[:blending_factor] = Vector{Float64}(undef, dg.n_elements)
    #end
    #if (!haskey(dg.element_variables, :blending_factor_tmp) ||
    #    length(dg.element_variables[:blending_factor_tmp]) != dg.n_elements)
    #  dg.element_variables[:blending_factor_tmp] = Vector{Float64}(undef, dg.n_elements)
    #end

    #alpha1     = dg.element_variables[:blending_factor]
    #alpha1_tmp = dg.element_variables[:blending_factor_tmp]
    #calc_blending_factors!(alpha1, alpha1_tmp, dg.elements.u, dg.shock_alpha_max, dg.shock_alpha_min, true,
    #                       dg.shock_indicator_variable, dg.thread_cache, dg)

    # Iterate over all elements
    for element_id in 1:dg.n_elements
      cell_id = dg.elements.cell_ids[element_id]
      actual_level = mesh.tree.levels[cell_id]
      target_level = actual_level
      # adapt for the amr indicator
      if alpha[element_id] >= blending_factor_threshold1
        target_level = max_level
      elseif alpha[element_id] <= blending_factor_threshold2
        target_level = base_level
      end
      # make sure that a highly troubled shock cell is not coarsened
      #if isapprox.(dg.shock_alpha_max, alpha1[element_id], atol=1e-12)
      #  target_level = max_level
      #end
      # Compare target level with actual level to set indicator
      if actual_level < target_level
        lambda[element_id] = 1.0
      elseif actual_level > target_level
        lambda[element_id] = -1.0
      else
        lambda[element_id] = 0.0
      end

    end
  elseif dg.amr_indicator === :density_pulse
    # Works with initial_conditions_density_pulse for compressible Euler equations
    base_level = 4
    max_level = 6
    threshold_high = 1.3
    threshold_low = 1.05

    # Iterate over all elements
    for element_id in 1:dg.n_elements
      # Determine target level from peak value
      peak = maximum(dg.elements.u[1, :, :, :, element_id])
      if peak > threshold_high
        target_level = max_level
      elseif peak > threshold_low
        target_level = max_level - 1
      else
        target_level = base_level
      end

      # Compare target level with actual level to set indicator
      cell_id = dg.elements.cell_ids[element_id]
      actual_level = mesh.tree.levels[cell_id]
      if actual_level < target_level
        lambda[element_id] = 1.0
      elseif actual_level > target_level
        lambda[element_id] = -1.0
      else
        lambda[element_id] = 0.0
      end
    end
  elseif dg.amr_indicator === :sedov_self_gravity
    base_level = 2
    max_level = parameter("max_refinement_level", 7)::Int # TODO: This is just so we can test this, but we should have a global switch...
    blending_factor_threshold = 0.0003

    # (Re-)initialize element variable storage for blending factor
    if (!haskey(dg.element_variables, :amr_indicator_values) ||
        length(dg.element_variables[:amr_indicator_values]) != dg.n_elements)
      dg.element_variables[:amr_indicator_values] = Vector{Float64}(undef, dg.n_elements)
    end
    if (!haskey(dg.element_variables, :amr_indicator_values_tmp) ||
        length(dg.element_variables[:amr_indicator_values_tmp]) != dg.n_elements)
      dg.element_variables[:amr_indicator_values_tmp] = Vector{Float64}(undef, dg.n_elements)
    end

    alpha     = dg.element_variables[:amr_indicator_values]
    alpha_tmp = dg.element_variables[:amr_indicator_values_tmp]
    calc_blending_factors!(alpha, alpha_tmp, dg.elements.u, dg.amr_alpha_max, dg.amr_alpha_min, dg.amr_alpha_smooth,
                           density_pressure, dg.thread_cache, dg)

    # Iterate over all elements
    for element_id in 1:dg.n_elements
      if alpha[element_id] > blending_factor_threshold
        target_level = max_level
      else
        target_level = base_level
      end

      # Compare target level with actual level to set indicator
      cell_id = dg.elements.cell_ids[element_id]
      actual_level = mesh.tree.levels[cell_id]
      if actual_level < target_level
        lambda[element_id] = 1.0
      elseif actual_level > target_level
        lambda[element_id] = -1.0
      else
        lambda[element_id] = 0.0
      end
    end
  else
    error("unknown AMR indicator '$(dg.amr_indicator)'")
  end

  return lambda
end


text:
@testset "Sensibility checks" begin
    @testset "Sanity checks on partitioning" begin
        st = Date(2018, 1, 1)
        ed = Date(2019, 1, 1)

        datesets = (
            st:Day(1):ed,
            st..ed,
            AnchoredInterval{Day(365), Date}(st),
        )

        selectors = (
            NoneSelector(),
            PeriodicSelector(Week(2), Week(1)),
            RandomSelector(42),
        )

        @testset "$(repr(selector))" for selector in selectors, dateset in datesets
            a, b = partition(dateset, selector)
            @test all(in(dateset), a)
            @test all(in(dateset), b)
            @test isempty(intersect(a, b))

            if dateset isa AbstractInterval
                 _dateset = DateSelectors._interval2daterange(dateset)
                 @test sort(union(a, b)) == collect(_dateset)
            else
                @test sort(union(a, b)) == collect(dateset)
            end

        end
    end

    @testset "Vector of days is not allowed" begin
        date_range = collect(Date(2019, 1, 1):Day(1):Date(2019, 2, 1))
        @test_throws MethodError partition(date_range, NoneSelector())
    end

    @testset "Weekly intervals are not allowed" begin
        st = Date(2019, 1, 1)
        ed = Date(2019, 2, 1)
        weekly_dates = st:Week(1):ed

        for selector in (
            NoneSelector(),
            PeriodicSelector(Week(2), Week(1)),
            RandomSelector(42),
        )
            @test_throws MethodError partition(weekly_dates, selector)
        end
    end
end


text:
module Bridges

using VecMathOptInterface
const MOI = VecMathOptInterface
const MOIU = MOI.Utilities

const CI = MOI.ConstraintIndex

include("bridge.jl")
include("set_map.jl")
include("bridge_optimizer.jl")

include("Variable/Variable.jl")
include("Constraint/Constraint.jl")
include("Objective/Objective.jl")

include("flip_sign.jl")
include("soc_rsoc.jl")

include("lazy_bridge_optimizer.jl")
include("debug.jl")

"""
    full_bridge_optimizer(model::MOI.ModelLike, ::Type{T}) where {T}

Returns a [`LazyBridgeOptimizer`](@ref) bridging `model` for every bridge
defined in this package (see below for the few exceptions) and for the
coefficient type `T` in addition to the bridges in the list returned by
`MOI.get(model, MOI.Bridges.ListOfNonstandardBridges{T}())`.

See also [`ListOfNonstandardBridges`](@ref).

!!! note
    The following bridges are not added by `full_bridge_optimizer` except if
    they are in the list returned by `MOI.get(model, MOI.Bridges.ListOfNonstandardBridges{T}())`
    (see the docstrings of the corresponding bridge for the reason they are not
    added):
    * [`Constraint.SOCtoNonConvexQuadBridge`](@ref),
      [`Constraint.RSOCtoNonConvexQuadBridge`](@ref) and
      [`Constraint.SOCtoPSDBridge`](@ref).
    * The subtypes of [`Constraint.AbstractToIntervalBridge`](@ref) (i.e.
      [`Constraint.GreaterToIntervalBridge`](@ref) and
      [`Constraint.LessToIntervalBridge`](@ref)) if `T` is not a subtype of
      `AbstractFloat`.
"""
function full_bridge_optimizer(model::MOI.ModelLike, ::Type{T}) where {T}
    bridged_model = LazyBridgeOptimizer(model)
    for BT in MOI.get(model, ListOfNonstandardBridges{T}())
        add_bridge(bridged_model, BT)
    end
    Variable.add_all_bridges(bridged_model, T)
    Constraint.add_all_bridges(bridged_model, T)
    Objective.add_all_bridges(bridged_model, T)
    return bridged_model
end

"""
    ListOfNonstandardBridges{T}() <: MOI.AbstractOptimizerAttribute

Any optimizer can be wrapped in a [`LazyBridgeOptimizer`](@ref) using
[`full_bridge_optimizer`](@ref). However, by default [`LazyBridgeOptimizer`](@ref)
uses a limited set of bridges that are:

  1. implemented in `MOI.Bridges`
  2. generally applicable for all optimizers.

For some optimizers however, it is useful to add additional bridges, such as
those that are implemented in external packages (e.g., within the solver package
itself) or only apply in certain circumstances (e.g.,
[`Constraint.SOCtoNonConvexQuadBridge`](@ref)).

Such optimizers should implement the `ListOfNonstandardBridges` attribute to
return a vector of bridge types that are added by [`full_bridge_optimizer`](@ref)
in addition to the list of default bridges.

Note that optimizers implementing `ListOfNonstandardBridges` may require
package-specific functions or sets to be used if the non-standard bridges
are not added. Therefore, you are recommended to use
`model = MOI.instantiate(Package.Optimizer; with_bridge_type = T)` instead of
`model = MOI.instantiate(Package.Optimizer)`. See
[`VecMathOptInterface.instantiate`](@ref).

## Examples

### An optimizer using a non-default bridge in `MOI.Bridges`

Solvers supporting [`MOI.ScalarQuadraticFunction`](@ref) can support
[`MOI.SecondOrderCone`](@ref) and [`MOI.RotatedSecondOrderCone`](@ref) by
defining:
```julia
function MOI.get(::MyQuadraticOptimizer, ::ListOfNonstandardBridges{Float64})
    return Type[
        MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64},
        MOI.Bridges.Constraint.RSOCtoNonConvexQuadBridge{Float64},
    ]
end
```

### An optimizer defining an internal bridge

Suppose an optimizer can exploit specific structure of a constraint, e.g., it
can exploit the structure of the matrix `A` in the linear system of equations
`A * x = b`.

The optimizer can define the function:
```julia
struct MatrixAffineFunction{T} <: MOI.AbstractVectorFunction
    A::SomeStructuredMatrixType{T}
    b::Vector{T}
end
```
and then a bridge
```julia
struct MatrixAffineFunctionBridge{T} <: MOI.Constraint.AbstractBridge
    # ...
end
# ...
```
from `VectorAffineFunction{T}` to the `MatrixAffineFunction`. Finally, it
defines:
```julia
function MOI.get(::Optimizer{T}, ::ListOfNonstandardBridges{T}) where {T}
    return Type[MatrixAffineFunctionBridge{T}]
end
```
"""
struct ListOfNonstandardBridges{T} <: MOI.AbstractOptimizerAttribute end

attribute_value_type(::ListOfNonstandardBridges) = Vector{Type}

MOI.is_copyable(::ListOfNonstandardBridges) = false

MOI.get_fallback(model::MOI.ModelLike, ::ListOfNonstandardBridges) = Type[]

print_num_bridges(io::IO, ::Variable.EmptyMap) = nothing

print_num_bridges(io::IO, ::Constraint.EmptyMap) = nothing

print_num_bridges(io::IO, ::Objective.EmptyMap) = nothing

function print_num_bridges(io::IO, B::Variable.Map)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    n = length(B)
    print(io, "\n$(indent)with $(n) variable bridge$(s(n))")
    return
end

function print_num_bridges(io::IO, B::Constraint.Map)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    n = length(B)
    print(io, "\n$(indent)with $(n) constraint bridge$(s(n))")
    return
end

function print_num_bridges(io::IO, B::Objective.Map)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    n = length(B)
    print(io, "\n$(indent)with $(n) objective bridge$(s(n))")
    return
end

function Base.show(io::IO, B::AbstractBridgeOptimizer)
    MOIU.print_with_acronym(io, summary(B))
    print_num_bridges(io, Variable.bridges(B))
    print_num_bridges(io, Constraint.bridges(B))
    print_num_bridges(io, Objective.bridges(B))
    if :model in propertynames(B)
        indent = " "^get(io, :indent, 0)
        print(io, "\n$(indent)with inner model ")
        show(IOContext(io, :indent => get(io, :indent, 0) + 2), B.model)
    end
    return
end

end


text:
using Turing: gradient_logp, ForwardDiffAD
using Test

"""
    test_reverse_mode_ad(forward, f, ȳ, x...; rtol=1e-6, atol=1e-6)

Check that the reverse-mode sensitivities produced by an AD library are correct for `f`
at `x...`, given sensitivity `ȳ` w.r.t. `y = f(x...)` up to `rtol` and `atol`.
"""
function test_reverse_mode_ad( f, ȳ, x...; rtol=1e-6, atol=1e-6)
    # Perform a regular forwards-pass.
    y = f(x...)

    # Use Tracker to compute reverse-mode sensitivities.
    y_tracker, back_tracker = Tracker.forward(f, x...)
    x̄s_tracker = back_tracker(ȳ)

    # Use Zygote to compute reverse-mode sensitivities.
    y_zygote, back_zygote = Zygote.pullback(f, x...)
    x̄s_zygote = back_zygote(ȳ)

    test_rd = length(x) == 1 && y isa Number
    if test_rd
        # Use ReverseDiff to compute reverse-mode sensitivities.
        if x[1] isa Array
            x̄s_rd = similar(x[1])
            tp = ReverseDiff.GradientTape(x -> f(x), x[1])
            ReverseDiff.gradient!(x̄s_rd, tp, x[1])
            x̄s_rd .*= ȳ
            y_rd = ReverseDiff.value(tp.output)
            @assert y_rd isa Number
        else
            x̄s_rd = [x[1]]
            tp = ReverseDiff.GradientTape(x -> f(x[1]), [x[1]])
            ReverseDiff.gradient!(x̄s_rd, tp, [x[1]])
            y_rd = ReverseDiff.value(tp.output)[1]
            x̄s_rd = x̄s_rd[1] * ȳ
            @assert y_rd isa Number
        end
    end

    # Use finite differencing to compute reverse-mode sensitivities.
    x̄s_fdm = FDM.j′vp(central_fdm(5, 1), f, ȳ, x...)

    # Check that Tracker forwards-pass produces the correct answer.
    @test isapprox(y, Tracker.data(y_tracker), atol=atol, rtol=rtol)

    # Check that Zygpte forwards-pass produces the correct answer.
    @test isapprox(y, y_zygote, atol=atol, rtol=rtol)

    if test_rd
        # Check that ReverseDiff forwards-pass produces the correct answer.
        @test isapprox(y, y_rd, atol=atol, rtol=rtol)
    end

    # Check that Tracker reverse-mode sensitivities are correct.
    @test all(zip(x̄s_tracker, x̄s_fdm)) do (x̄_tracker, x̄_fdm)
        isapprox(Tracker.data(x̄_tracker), x̄_fdm; atol=atol, rtol=rtol)
    end

    # Check that Zygote reverse-mode sensitivities are correct.
    @test all(zip(x̄s_zygote, x̄s_fdm)) do (x̄_zygote, x̄_fdm)
        isapprox(x̄_zygote, x̄_fdm; atol=atol, rtol=rtol)
    end

    if test_rd
        # Check that ReverseDiff reverse-mode sensitivities are correct.
        @test isapprox(x̄s_rd, x̄s_zygote[1]; atol=atol, rtol=rtol)
    end
end

function test_model_ad(model, f, syms::Vector{Symbol})
    # Set up VI.
    vi = Turing.VarInfo(model)

    # Collect symbols.
    vnms = Vector(undef, length(syms))
    vnvals = Vector{Float64}()
    for i in 1:length(syms)
        s = syms[i]
        vnms[i] = getfield(vi.metadata, s).vns[1]

        vals = getval(vi, vnms[i])
        for i in eachindex(vals)
            push!(vnvals, vals[i])
        end
    end

    spl = SampleFromPrior()
    _, ∇E = gradient_logp(ForwardDiffAD{1}(), vi[spl], vi, model)
    grad_Turing = sort(∇E)

    # Call ForwardDiff's AD
    grad_FWAD = sort(ForwardDiff.gradient(f, vec(vnvals)))

    # Compare result
    @test grad_Turing ≈ grad_FWAD atol=1e-9
end


text:
"Construct state from a list of terms (e.g. initial predicates and fluents)."
function State(terms::Vector{<:Term}, types::Vector{<:Term}=Term[])
    state = State(Set{Term}(types), Set{Term}(), Dict{Symbol,Any}())
    for t in terms
        if t.name == :(==)
            # Initialize fluents
            term, val = t.args[1], t.args[2]
            @assert !isa(term, Var) "Initial terms cannot be unbound variables."
            @assert isa(val, Const) "Terms must be initialized to constants."
            state[term] = val.name
        else
            push!(state.facts, t)
        end
    end
    return state
end

"Returns the list of all terms in a state."
function get_terms(state::State)
    return Term[get_types(state); get_facts(state); get_fluents(state)]
end

"Returns the list of all type declarations in a state."
function get_types(state::State)
    return collect(state.types)
end

"Returns the list of all facts in a state."
function get_facts(state::State)
    return collect(state.facts)
end

"Returns the list of all fluent terms (but not their values) in a state."
function get_fluents(state::State)
    terms = Term[]
    for (name, val) in state.fluents
        if isa(val, Dict)
            append!(terms, [Compound(name, Const.(collect(args)))
                            for args in keys(val)])
        else
            push!(terms, Const(name))
        end
    end
    return terms
end

"Returns the list of all fluent assignments in a state."
function get_assignments(state::State)
    terms = Term[]
    for (name, val) in state.fluents
        if isa(val, Dict)
            for (args, v) in val
                fluent = Compound(name, Const.(collect(args)))
                assignment = @julog(:fluent == $v)
                push!(terms, assignment)
            end
        else
            assignment = @julog($name == $val)
            push!(terms, assignment)
        end
    end
    return terms
end

Base.copy(s::State) =
    State(copy(s.types), copy(s.facts), deepcopy(s.fluents))
Base.:(==)(s1::State, s2::State) =
    s1.types == s2.types && s1.facts == s2.facts && s1.fluents == s2.fluents
Base.hash(s::State, h::UInt) =
    hash(s.fluents, hash(s.facts, hash(s.types, h)))
Base.issubset(s1::State, s2::State) =
    s1.types ⊆ s2.types && s1.facts ⊆ s2.facts &&
    (let f1 = get_fluents(s1), f2 = get_fluents(s2)
        all(f in f2 && s1[f] == s2[f] for f in f1)
    end)

"Access the value of a fluent or fact in a state."
Base.getindex(state::State, term::Term) =
    evaluate(term, state; as_const=false)
Base.getindex(state::State, term::String) =
    evaluate(Parser.parse_formula(term), state; as_const=false)
Base.getindex(state::State, term::Symbol) =
    evaluate(Const(term), state; as_const=false)
Base.getindex(state::State, term::Expr) =
    evaluate(eval(Julog.parse_term(term, identity)), state; as_const=false)
Base.getindex(state::State, domain::Domain, term::Term) =
    evaluate(term, state, domain; as_const=false)
Base.getindex(state::State, domain::Domain, term::String) =
    evaluate(Parser.parse_formula(term), state, domain; as_const=false)
Base.getindex(state::State, domain::Domain, term::Symbol) =
    evaluate(Const(term), state, domain; as_const=false)
Base.getindex(state::State, domain::Domain, term::Expr) =
    evaluate(eval(Julog.parse_term(term, identity)), state, domain; as_const=false)

"Set the value of a fluent or fact in a state."
Base.setindex!(state::State, val::Bool, term::Const) =
    (if val push!(state.facts, term) else delete!(state.facts, term) end; val)
Base.setindex!(state::State, val::Bool, term::Compound) =
    (if val push!(state.facts, term) else delete!(state.facts, term) end; val)
Base.setindex!(state::State, val::Any, term::Const) =
    (state.fluents[term.name] = val)
Base.setindex!(state::State, val::Any, term::Compound) =
    (d = get!(state.fluents, term.name, Dict());
     d[Tuple(a.name for a in term.args)] = val)
Base.setindex!(state::State, val::Any, term::Symbol) =
    Base.setindex!(state, val, Const(term))
Base.setindex!(state::State, val::Any, term::Expr) =
    Base.setindex!(state, val, eval(Julog.parse_term(term, identity)))