text:
mutable struct mcpop_data
    n # number of all variables
    nb # number of binary variables
    m # number of all constraints
    numeq # number of equality constraints
    supp # support data
    coe # coefficient data
    basis # monomial bases
    rlorder # relaxation order
    ksupp # extending support at the k-th step
    cql # number of cliques
    cliques # cliques of variables
    cliquesize # numbers of cliques
    J # constraints associated to each clique
    ncc # constraints associated to no clique
    sb # sizes of different blocks
    numb # numbers of different blocks
    blocks # block structure
    cl # numbers of blocks
    blocksize # sizes of blocks
    Mmatrix # Moment matrix
    solver # SDP solver
    tol # tolerance to certify global optimality
    flag # 0 if global optimality is certified; 1 otherwise
end

"""
    opt,sol,data = cs_tssos_first(pop, x, d; nb=0, numeq=0, foc=100, CS="MF", minimize=false,
    assign="first", TS="block", merge=false, md=3, solver="Mosek", QUIET=false, solve=true, solution=false,
    MomentOne=true, tol=1e-4)

Compute the first step of the CS-TSSOS hierarchy for constrained polynomial optimization with
relaxation order `d`. Return the optimum, the (near) optimal solution (if `solution=true`) and
other auxiliary data.

# Arguments
- `pop`: the vector of the objective function, inequality constraints, and equality constraints.
- `x`: the set of variables.
- `d`: the relaxation order of the moment-SOS hierarchy.
- `nb`: the number of binary variables in `x`.
- `md`: the tunable parameter for merging blocks.
- `numeq`: the number of equality constraints.
"""
function cs_tssos_first(pop, x, d; nb=0, numeq=0, foc=100, CS="MF", minimize=false, assign="first", TS="block",
    merge=false, md=3, solver="Mosek", tune=false, QUIET=false, solve=true, solution=false, MomentOne=true, Mommat=false, tol=1e-4)
    n,supp,coe = polys_info(pop, x)
    opt,sol,data = cs_tssos_first(supp, coe, n, d, numeq=numeq, nb=nb, foc=foc, CS=CS, minimize=minimize, assign=assign, TS=TS,
    merge=merge, md=md, QUIET=QUIET, solver=solver, tune=tune, solve=solve, solution=solution, MomentOne=MomentOne, Mommat=Mommat, tol=tol)
    return opt,sol,data
end

function polys_info(pop, x)
    n = length(x)
    m = length(pop)-1
    coe = Vector{Vector{Float64}}(undef, m+1)
    supp = Vector{Vector{Vector{UInt16}}}(undef, m+1)
    for k = 1:m+1
        mon = monomials(pop[k])
        coe[k] = coefficients(pop[k])
        lm = length(mon)
        supp[k] = [UInt16[] for i=1:lm]
        for i = 1:lm
            ind = mon[i].z .> 0
            vars = mon[i].vars[ind]
            exp = mon[i].z[ind]
            for j = 1:length(vars)
                l = ncbfind(x, n, vars[j])
                append!(supp[k][i], l*ones(UInt16, exp[j]))
            end
        end
    end
    return n,supp,coe
end

"""
    opt,sol,data = cs_tssos_first(supp::Vector{Vector{Vector{UInt16}}}, coe, n, d; numeq=0,
    nb=0, foc=100, CS="MF", minimize=false, assign="first", TS="block", merge=false, md=3,
    QUIET=false, solver="Mosek", solve=true, solution=false, MomentOne=true, tol=1e-4)

Compute the first step of the CS-TSSOS hierarchy for constrained polynomial optimization with
relaxation order `d`. Here the polynomial optimization problem is defined by `supp` and `coe`,
corresponding to the supports and coeffients of `pop` respectively.
Return the optimum, the (near) optimal solution (if `solution=true`) and other auxiliary data.

# Arguments
- `supp`: the supports of the polynomial optimization problem.
- `coe`: the coeffients of the polynomial optimization problem.
- `d`: the relaxation order of the moment-SOS hierarchy.
- `numeq`: the number of equality constraints.
"""
function cs_tssos_first(supp::Vector{Vector{Vector{UInt16}}}, coe, n, d; numeq=0, nb=0,
    foc=100, CS="MF", minimize=false, assign="first", TS="block", merge=false, md=3, QUIET=false,
    solver="Mosek", tune=false, solve=true, solution=false, MomentOne=true, Mommat=false, tol=1e-4)
    println("************************TSSOS************************")
    println("TSSOS is launching...")
    m = length(supp)-1
    supp[1],coe[1] = resort(supp[1], coe[1])
    dg = [maximum(length.(supp[i])) for i=2:m+1]
    time = @elapsed begin
    cliques,cql,cliquesize = clique_decomp(n, m, dg, supp, order=d, alg=CS, minimize=minimize)
    end
    if CS != false && QUIET == false
        mc = maximum(cliquesize)
        println("Obtained the variable cliques in $time seconds. The maximal size of cliques is $mc.")
    end
    J,ncc = assign_constraint(m, supp, cliques, cql, cliquesize, assign=assign)
    rlorder = init_order(dg, J, cliquesize, cql, foc=foc, order=d)
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,basis,status = get_cblocks_mix(dg, J, rlorder, m, supp, cliques, cql, cliquesize, nb=nb, TS=TS, merge=merge, md=md)
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp,moment = blockcpop_mix(n, m, supp, coe, basis, cliques, cql, cliquesize, J, ncc, blocks, cl, blocksize, numeq=numeq, nb=nb, QUIET=QUIET, TS=TS, solver=solver, tune=tune, solve=solve, solution=solution, MomentOne=MomentOne, Mommat=Mommat)
    data = mcpop_data(n, nb, m, numeq, supp, coe, basis, rlorder, ksupp, cql, cliques, cliquesize, J, ncc, sb, numb, blocks, cl, blocksize, moment, solver, tol, 1)
    sol = nothing
    if solution == true
        sol,data.flag = approx_sol(opt, moment, n, cliques, cql, cliquesize, supp, coe, numeq=numeq, tol=tol)
        if data.flag == 1
            sol,ub,gap = refine_sol(opt, sol, data, QUIET=true)
            if gap != nothing
                if gap < tol
                    data.flag = 0
                else
                    println("Found a local optimal solution giving an upper bound: $ub and a relative optimality gap: $gap.")
                end
            end
        end
    end
    return opt,sol,data
end

"""
    opt,sol,data = cs_tssos_higher!(data; TS="block", merge=false, md=3, QUIET=false, solve=true,
    solution=false, MomentOne=false)

Compute higher steps of the CS-TSSOS hierarchy.
Return the optimum, the (near) optimal solution (if `solution=true`) and other auxiliary data.
"""
function cs_tssos_higher!(data; TS="block", merge=false, md=3, QUIET=false, solve=true, tune=false, solution=false, ipart=true, MomentOne=false, Mommat=false)
    n = data.n
    nb = data.nb
    m = data.m
    numeq = data.numeq
    supp = data.supp
    coe = data.coe
    basis = data.basis
    rlorder = data.rlorder
    ksupp = data.ksupp
    cql = data.cql
    cliques = data.cliques
    cliquesize = data.cliquesize
    J = data.J
    ncc = data.ncc
    blocks = data.blocks
    cl = data.cl
    blocksize = data.blocksize
    sb = data.sb
    numb = data.numb
    solver = data.solver
    tol = data.tol
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,basis,status = get_cblocks_mix([], J, rlorder, m, supp, cliques, cql, cliquesize, tsupp=ksupp, basis=basis, blocks=blocks, cl=cl, blocksize=blocksize, sb=sb, numb=numb, nb=nb, TS=TS, merge=merge, md=md)
    end
    opt = nothing
    sol = nothing
    if status == 1
        if TS != false && QUIET == false
            mb = maximum(maximum.(sb))
            println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
        end
        opt,ksupp,moment = blockcpop_mix(n, m, supp, coe, basis, cliques, cql, cliquesize, J, ncc, blocks, cl, blocksize, numeq=numeq, nb=nb, QUIET=QUIET, solver=solver, solve=solve, tune=tune, solution=solution, ipart=ipart, MomentOne=MomentOne, Mommat=Mommat)
        if solution == true
            sol,data.flag = approx_sol(opt, moment, n, cliques, cql, cliquesize, supp, coe, numeq=numeq, tol=tol)
            if data.flag == 1
                sol,ub,gap = refine_sol(opt, sol, data, QUIET=true)
                if gap != nothing
                    if gap < tol
                        data.flag = 0
                    else
                        println("Found a local optimal solution giving an upper bound: $ub and a relative optimality gap: $gap.")
                    end
                end
            end
        end
        data.ksupp = ksupp
        data.blocks = blocks
        data.cl = cl
        data.blocksize = blocksize
        data.Mmatrix = moment
        data.sb = sb
        data.numb = numb
    else
        println("No higher CS-TSSOS hierarchy!")
    end
    return opt,sol,data
end

function blockcpop_mix(n, m, supp::Vector{Vector{Vector{UInt16}}}, coe, basis, cliques, cql, cliquesize,
    J, ncc, blocks, cl, blocksize; numeq=0, nb=0, QUIET=false, TS="block", solver="Mosek", tune=false,
    solve=true, solution=false, ipart=false, MomentOne=false, Mommat=false)
    tsupp = Vector{UInt16}[]
    for i = 1:cql, j = 1:cl[i][1], k = 1:blocksize[i][1][j], r = k:blocksize[i][1][j]
        @inbounds bi = sadd(basis[i][1][blocks[i][1][j][k]], basis[i][1][blocks[i][1][j][r]], nb=nb)
        push!(tsupp, bi)
    end
    if (MomentOne == true || solution == true) && TS != false
        ksupp = copy(tsupp)
    end
    for i = 1:cql, (j, w) in enumerate(J[i])
        for l = 1:cl[i][j+1], t = 1:blocksize[i][j+1][l], r = t:blocksize[i][j+1][l], s = 1:length(supp[w+1])
            ind1 = blocks[i][j+1][l][t]
            ind2 = blocks[i][j+1][l][r]
            @inbounds bi = sadd(sadd(basis[i][j+1][ind1], supp[w+1][s], nb=nb), basis[i][j+1][ind2], nb=nb)
            push!(tsupp, bi)
        end
    end
    for i ∈ ncc, j = 1:length(supp[i+1])
        push!(tsupp, supp[i+1][j])
    end
    if (MomentOne == true || solution == true) && TS != false
        for i = 1:cql, j = 1:cliquesize[i]
            push!(tsupp, [cliques[i][j]])
            for k = j+1:cliquesize[i]
                push!(tsupp, [cliques[i][j];cliques[i][k]])
            end
        end
    end
    sort!(tsupp)
    unique!(tsupp)
    if (MomentOne == true || solution == true) && TS != false
        sort!(ksupp)
        unique!(ksupp)
    else
        ksupp = tsupp
    end
    objv = nothing
    moment = nothing
    if solve == true
        ltsupp = length(tsupp)
        if QUIET == false
            println("Assembling the SDP...")
            println("There are $ltsupp affine constraints.")
        end
        if solver == "Mosek"
            model = Model(optimizer_with_attributes(Mosek.Optimizer))
            if tune == true
                set_optimizer_attributes(model,
                "MSK_DPAR_INTPNT_CO_TOL_MU_RED" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_INFEAS" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_REL_GAP" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_DFEAS" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_PFEAS" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_NEAR_REL" => 1e6,
                "MSK_IPAR_BI_IGNORE_NUM_ERROR" => 1,
                "MSK_DPAR_BASIS_TOL_X" => 1e-3,
                "MSK_DPAR_BASIS_TOL_S" => 1e-3,
                "MSK_DPAR_BASIS_REL_TOL_S" => 1e-5)
            end
        elseif solver == "COSMO"
            model = Model(optimizer_with_attributes(COSMO.Optimizer, "max_iter" => 10000))
        elseif solver == "SDPT3"
            model = Model(optimizer_with_attributes(SDPT3.Optimizer))
        else
            @error "The solver is currently not supported!"
            return nothing,nothing,nothing
        end
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        time = @elapsed begin
        cons = [AffExpr(0) for i=1:ltsupp]
        for i = 1:cql
            if (MomentOne == true || solution == true) && TS != false
                bs = cliquesize[i]+1
                pos = @variable(model, [1:bs, 1:bs], PSD)
                for t = 1:bs, r = t:bs
                    if t == 1 && r == 1
                        bi = UInt16[]
                    elseif t == 1 && r > 1
                        bi = [cliques[i][r-1]]
                    else
                        bi = sadd(cliques[i][t-1], cliques[i][r-1], nb=nb)
                    end
                    Locb = bfind(tsupp, ltsupp, bi)
                    if t == r
                        @inbounds add_to_expression!(cons[Locb], pos[t,r])
                    else
                        @inbounds add_to_expression!(cons[Locb], 2, pos[t,r])
                    end
                end
            end
            for l = 1:cl[i][1]
                @inbounds bs = blocksize[i][1][l]
                if bs == 1
                    pos = @variable(model, lower_bound=0)
                    @inbounds bi = sadd(basis[i][1][blocks[i][1][l][1]], basis[i][1][blocks[i][1][l][1]], nb=nb)
                    Locb = bfind(tsupp, ltsupp, bi)
                    @inbounds add_to_expression!(cons[Locb], pos)
                else
                    pos = @variable(model, [1:bs, 1:bs], PSD)
                    for t = 1:bs, r = t:bs
                        @inbounds ind1 = blocks[i][1][l][t]
                        @inbounds ind2 = blocks[i][1][l][r]
                        @inbounds bi = sadd(basis[i][1][ind1], basis[i][1][ind2], nb=nb)
                        Locb = bfind(tsupp, ltsupp, bi)
                        if t == r
                            @inbounds add_to_expression!(cons[Locb], pos[t,r])
                        else
                            @inbounds add_to_expression!(cons[Locb], 2, pos[t,r])
                        end
                    end
                end
            end
        end
        for i in ncc
            if i <= m-numeq
                pos = @variable(model, lower_bound=0)
            else
                pos = @variable(model)
            end
            for j = 1:length(supp[i+1])
                Locb = bfind(tsupp, ltsupp, supp[i+1][j])
                @inbounds add_to_expression!(cons[Locb], coe[i+1][j], pos)
            end
        end
        for i = 1:cql, (j, w) in enumerate(J[i]), l = 1:cl[i][j+1]
            bs = blocksize[i][j+1][l]
            if bs == 1
                if w <= m-numeq
                    pos = @variable(model, lower_bound=0)
                else
                    pos = @variable(model)
                end
                ind = blocks[i][j+1][l][1]
                for s = 1:length(supp[w+1])
                    @inbounds bi = sadd(sadd(basis[i][j+1][ind], supp[w+1][s], nb=nb), basis[i][j+1][ind], nb=nb)
                    Locb = bfind(tsupp, ltsupp, bi)
                    @inbounds add_to_expression!(cons[Locb], coe[w+1][s], pos)
                end
            else
                if w <= m-numeq
                    pos = @variable(model, [1:bs, 1:bs], PSD)
                else
                    pos = @variable(model, [1:bs, 1:bs], Symmetric)
                end
                for t = 1:bs, r = t:bs
                    ind1 = blocks[i][j+1][l][t]
                    ind2 = blocks[i][j+1][l][r]
                    for s = 1:length(supp[w+1])
                        @inbounds bi = sadd(sadd(basis[i][j+1][ind1], supp[w+1][s], nb=nb), basis[i][j+1][ind2], nb=nb)
                        Locb = bfind(tsupp, ltsupp, bi)
                        if t == r
                            @inbounds add_to_expression!(cons[Locb], coe[w+1][s], pos[t,r])
                        else
                            @inbounds add_to_expression!(cons[Locb], 2*coe[w+1][s], pos[t,r])
                        end
                    end
                end
            end
        end
        bc = zeros(ltsupp)
        for i = 1:length(supp[1])
            Locb = bfind(tsupp, ltsupp, supp[1][i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,nothing,nothing
            else
               bc[Locb] = coe[1][i]
            end
        end
        @variable(model, lower)
        if solution == true || Mommat == true
            cons[1] += lower
            @constraint(model, con[i=1:ltsupp], cons[i]==bc[i])
        else
            @constraint(model, cons[2:end].==bc[2:end])
            @constraint(model, cons[1]+lower==bc[1])
        end
        @objective(model, Max, lower)
        end
        if QUIET == false
            println("SDP assembling time: $time seconds.")
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
        if solution == true
            measure = -dual.(con)
            moment = get_moment(measure, tsupp, cliques, cql, cliquesize, nb=nb)
        end
        if Mommat == true
            measure = -dual.(con)
            moment = get_moment(measure, tsupp, cliques, cql, cliquesize, basis=basis, nb=nb)
        end
    end
    return objv,ksupp,moment
end

function init_order(dg, J, cliquesize, cql; foc=100, order="min")
    rlorder = ones(Int, cql)
    if order == "min"
        for i = 1:cql
            if !isempty(J[i])
                rlorder[i] = ceil(Int, maximum(dg[J[i]])/2)
            end
        end
    else
        for i = 1:cql
            if cliquesize[i] <= foc
                rlorder[i] = order
            end
        end
    end
    return rlorder
end

function get_cblocks_mix(dg, J, rlorder, m, supp::Vector{Vector{Vector{UInt16}}}, cliques, cql,
    cliquesize; tsupp=[], basis=[], blocks=[], cl=[], blocksize=[], sb=[], numb=[], TS="block",
    nb=0, merge=false, md=3)
    if isempty(tsupp)
        blocks = Vector{Vector{Vector{Vector{UInt16}}}}(undef, cql)
        cl = Vector{Vector{UInt16}}(undef, cql)
        blocksize = Vector{Vector{Vector{UInt16}}}(undef, cql)
        sb = Vector{Vector{UInt16}}(undef, cql)
        numb = Vector{Vector{UInt16}}(undef, cql)
        basis = Vector{Vector{Vector{Vector{UInt16}}}}(undef, cql)
        tsupp = copy(supp[1])
        for i = 2:m+1, j = 1:length(supp[i])
            push!(tsupp, supp[i][j])
        end
        sort!(tsupp)
        unique!(tsupp)
        flag = 1
    else
        flag = 0
    end
    status = ones(Int, cql)
    for i = 1:cql
        lc = length(J[i])
        nvar = cliquesize[i]
        ind = [issubset(tsupp[j], cliques[i]) for j=1:length(tsupp)]
        fsupp = copy(tsupp[ind])
        if flag == 1
            basis[i] = Vector{Vector{Vector{UInt16}}}(undef, lc+1)
            basis[i][1] = get_sbasis(cliques[i], rlorder[i], nb=nb)
            for j = 1:length(basis[i][1])
                push!(fsupp, sadd(basis[i][1][j], basis[i][1][j], nb=nb))
            end
            sort!(fsupp)
            unique!(fsupp)
            for s = 1:lc
                basis[i][s+1] = get_sbasis(cliques[i], rlorder[i]-ceil(Int, dg[J[i][s]]/2), nb=nb)
            end
            blocks[i] = Vector{Vector{Vector{UInt16}}}(undef, lc+1)
            cl[i] = Vector{UInt16}(undef, lc+1)
            blocksize[i] = Vector{Vector{UInt16}}(undef, lc+1)
            sb[i] = Vector{UInt16}(undef, lc+1)
            numb[i] = Vector{UInt16}(undef, lc+1)
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_cblocks(lc, fsupp, supp[J[i].+1], basis[i], TS=TS, nb=nb, QUIET=true, merge=merge, md=md)
        else
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_cblocks(lc, fsupp, supp[J[i].+1], basis[i], blocks=blocks[i], cl=cl[i], blocksize=blocksize[i], sb=sb[i], numb=numb[i], TS=TS, nb=nb, QUIET=true, merge=merge, md=md)
        end
    end
    return blocks,cl,blocksize,sb,numb,basis,maximum(status)
end

function assign_constraint(m, supp::Vector{Vector{Vector{UInt16}}}, cliques, cql, cliquesize; assign="first")
    J = [UInt32[] for i=1:cql]
    ncc = UInt32[]
    for i = 2:m+1
        rind = copy(supp[i][1])
        for j = 2:length(supp[i])
            append!(rind, supp[i][j])
        end
        unique!(rind)
        if assign == "first"
            ind = findfirst(k->issubset(rind, cliques[k]), 1:cql)
            if ind != nothing
                push!(J[ind], i-1)
            else
                push!(ncc, i-1)
            end
        else
            temp = UInt32[]
            for j = 1:cql
                if issubset(rind, cliques[j])
                    push!(temp, j)
                end
            end
            if !isempty(temp)
                if assign == "min"
                    push!(J[temp[argmin(cliquesize[temp])]], i-1)
                else
                    push!(J[temp[argmax(cliquesize[temp])]], i-1)
                end
            else
                push!(ncc, i-1)
            end
        end
    end
    return J,ncc
end

function get_sbasis(var, d; nb=0)
    n = length(var)
    lb = binomial(n+d, d)
    basis = Vector{Vector{UInt16}}(undef, lb)
    basis[1] = UInt16[]
    i = 0
    t = 1
    while i < d+1
        t += 1
        if sum(basis[t-1]) == var[n]*i
           if i < d
               basis[t] = var[1]*ones(UInt16, i+1)
           end
           i += 1
        else
            j = bfind(var, n, basis[t-1][1])
            basis[t] = copy(basis[t-1])
            ind = findfirst(x->basis[t][x]!=var[j], 1:length(basis[t]))
            if ind == nothing
                ind = length(basis[t])+1
            end
            if j != 1
                basis[t][1:ind-2] = var[1]*ones(UInt16, ind-2)
            end
            basis[t][ind-1] = var[j+1]
        end
    end
    if nb > 0
        ind = [!any([basis[i][j]==basis[i][j+1]&&basis[i][j]<=nb for j=1:length(basis[i])-1]) for i=1:lb]
        basis = basis[ind]
    end
    return basis
end

function get_graph(tsupp::Vector{Vector{UInt16}}, basis::Vector{Vector{UInt16}}; nb=0)
    lb = length(basis)
    G = SimpleGraph(lb)
    ltsupp = length(tsupp)
    for i = 1:lb, j = i+1:lb
        bi = sadd(basis[i], basis[j], nb=nb)
        if bfind(tsupp, ltsupp, bi)!=0
            add_edge!(G, i, j)
        end
    end
    return G
end

function get_cgraph(tsupp::Vector{Vector{UInt16}}, supp::Vector{Vector{UInt16}}, basis::Vector{Vector{UInt16}}; nb=0)
    lb = length(basis)
    ltsupp = length(tsupp)
    G = SimpleGraph(lb)
    for i = 1:lb, j = i+1:lb
        r = 1
        while r <= length(supp)
            bi = sadd(sadd(basis[i], supp[r], nb=nb), basis[j], nb=nb)
            if bfind(tsupp, ltsupp, bi)!=0
               break
            else
                r += 1
            end
        end
        if r <= length(supp)
            add_edge!(G, i, j)
        end
    end
    return G
end

function clique_decomp(n, m, dg, supp::Vector{Vector{Vector{UInt16}}}; order="min", alg="MF", minimize=false)
    if alg == false
        cliques = [UInt16[i for i=1:n]]
        cql = 1
        cliquesize = [n]
    else
        G = SimpleGraph(n)
        for i = 1:m+1
            if order == "min" || i == 1 || order == ceil(Int, dg[i-1]/2)
                for j = 1:length(supp[i])
                    add_clique!(G, unique(supp[i][j]))
                end
            else
                temp = copy(supp[i][1])
                for j = 2:length(supp[i])
                    append!(temp, supp[i][j])
                end
                add_clique!(G, unique(temp))
            end
        end
        if alg == "NC"
            cliques,cql,cliquesize = max_cliques(G)
        else
            cliques,cql,cliquesize = chordal_cliques!(G, method=alg, minimize=minimize)
        end
    end
    uc = unique(cliquesize)
    sizes=[sum(cliquesize.== i) for i in uc]
    println("------------------------------------------------------")
    println("The clique sizes of varibles:\n$uc\n$sizes")
    println("------------------------------------------------------")
    return cliques,cql,cliquesize
end

function sadd(a, b; nb=0)
    c = [a;b]
    sort!(c)
    if nb > 0
        i = 1
        while i < length(c)
            if c[i] <= nb
                if c[i] == c[i+1]
                    deleteat!(c, i:i+1)
                else
                    i += 1
                end
            else
                break
            end
        end
    end
    return c
end

function approx_sol(opt, moment, n, cliques, cql, cliquesize, supp, coe; numeq=0, tol=1e-4)
    qsol = Float64[]
    lcq = sum(cliquesize)
    A = zeros(lcq,n)
    q = 1
    for k = 1:cql
        cqs = cliquesize[k]
        F = eigen(moment[k], cqs+1:cqs+1)
        temp = sqrt(F.values[1])*F.vectors[:,1]
        if temp[1] == 0
            temp = zeros(cqs)
        else
            temp = temp[2:cqs+1]./temp[1]
        end
        append!(qsol, temp)
        for j = 1:cqs
            A[q,cliques[k][j]] = 1
            q += 1
        end
    end
    sol = (A'*A)\(A'*qsol)
    flag = 0
    upper_bound = seval(supp[1], coe[1], sol)
    gap = abs(upper_bound)>1 ? abs((opt-upper_bound)/upper_bound) : abs(opt-upper_bound)
    if gap >= tol
        flag = 1
    end
    m = length(supp)-1
    for i = 1:m-numeq
        if seval(supp[i+1], coe[i+1], sol) <= -tol
            flag = 1
        end
    end
    for i = m-numeq+1:m
        if abs(seval(supp[i+1], coe[i+1], sol)) >= tol
            flag = 1
        end
    end
    if flag == 0
        println("Global optimality certified!")
    end
    return sol,flag
end

function get_moment(measure, tsupp, cliques, cql, cliquesize; basis=[], nb=0)
    moment = Vector{Union{Float64, Symmetric{Float64}, Array{Float64,2}}}(undef, cql)
    ltsupp = length(tsupp)
    for i = 1:cql
        if basis == []
            lb = cliquesize[i]+1
        else
            lb = length(basis[i][1])
        end
        moment[i] = zeros(Float64, lb, lb)
        if basis == []
            for j = 1:lb, k = j:lb
                if j == 1
                    if k == 1
                        bi = UInt16[]
                    else
                        bi = [cliques[i][k-1]]
                    end
                else
                    bi = sadd(cliques[i][j-1], cliques[i][k-1], nb=nb)
                end
                Locb = bfind(tsupp, ltsupp, bi)
                moment[i][j,k] = measure[Locb]
            end
        else
            for j = 1:lb, k = j:lb
                bi = sadd(basis[i][1][j], basis[i][1][k], nb=nb)
                Locb = bfind(tsupp, ltsupp, bi)
                if Locb != 0
                    moment[i][j,k] = measure[Locb]
                end
            end
        end
        moment[i] = Symmetric(moment[i],:U)
    end
    return moment
end

function ncbfind(A, l, a)
    low = 1
    high = l
    while low <= high
        mid = Int(ceil(1/2*(low+high)))
        if A[mid] == a
           return mid
        elseif A[mid] < a
            high = mid-1
        else
            low = mid+1
        end
    end
    return 0
end

function seval(supp, coe, x)
    val = 0
    for i = 1:length(supp)
        if isempty(supp[i])
            temp = 1
        else
            temp = prod(x[supp[i][j]] for j=1:length(supp[i]))
        end
        val += coe[i]*temp
    end
    return val
end


text:
@testset "Cycles" begin
    completedg = CompleteDiGraph(4)
    pathdg = PathDiGraph(5)
    triangle = random_regular_graph(3, 2)
    quadrangle = random_regular_graph(4, 2)
    pentagon = random_regular_graph(5, 2)

    @test maxsimplecycles(pathdg) == 0
    @test maxsimplecycles(completedg) == 20
    @test maxsimplecycles(4) == 20
    @test maxsimplecycles(pathdg, false) == 84

    @test length(simplecycles(completedg)) == 20
    @test simplecycles(completedg) == @inferred(simplecycles_iter(completedg))
    @test simplecyclescount(completedg) == 20
    @test simplecycleslength(completedg) == ([0, 6, 8, 6], 20)

    @test simplecyclescount(pathdg) == 0
    @test length(simplecycles(pathdg)) == 0
    @test isempty(simplecycles(pathdg)) == true
    @test isempty(simplecycles_iter(pathdg)) == true
    @test simplecycleslength(pathdg) == (zeros(5), 0)

    @test simplecyclescount(completedg, 10) == 10
    @test simplecycleslength(completedg, 10)[2] == 10

    @test simplecyclescount(pathdg, 10) == 0
    @test isempty(simplecycles_iter(pathdg, 10)) == true
    @test simplecycleslength(pathdg, 10) == (zeros(5), 0)

    trianglelengths, triangletotal = simplecycleslength(DiGraph(triangle))
    @test sum(trianglelengths) == triangletotal

    quadranglelengths, quadrangletotal = simplecycleslength(DiGraph(quadrangle))
    @test sum(quadranglelengths) == quadrangletotal
    @test simplecycles(DiGraph(quadrangle)) == @inferred(simplecycles_iter(DiGraph(quadrangle)))

    pentagonlengths, pentagontotal = simplecycleslength(DiGraph(pentagon))
    @test sum(pentagonlengths) == pentagontotal

    selfloopg = DiGraph([
        0 1 0 0;
        0 0 1 0;
        1 0 1 0;
        0 0 0 1;
    ])
    cycles = simplecycles(selfloopg)
    @test [3] in cycles
    @test [4] in cycles
    @test [1, 2, 3] in cycles || [2, 3, 1] in cycles || [3, 1, 2] in cycles
    @test length(cycles) == 3

    cycles2 = simplecycles_iter(selfloopg)
    @test [3] in cycles2
    @test [4] in cycles2
    @test [1, 2, 3] in cycles2 || [2, 3, 1] in cycles2 || [3, 1, 2] in cycles2
    @test length(cycles2) == 3
end

ext: jl


text:
abstract type Aircraft end

# Mass properties
"""
    get_mass(ac::Aircraft)

Get aircraft mass (kg).
"""
function get_mass(ac::Aircraft) end

"""
    get_inertia_tensor(ac::Aircraft)

Get aircraft inertia tensor (kg·m²).
"""
function get_inertia_tensor(ac::Aircraft) end

"""
    get_xcg_mac(ac::Aircraft)

Get aircraft center of gravity position (units of mean aerodynamic chord).
"""
function get_xcg_mac(ac::Aircraft) end

# Geometric properties
"""
    get_surface(ac:Aircraft)

Get aircraft wing surface, ``S`` (m²).
"""
function get_surface(ac::Aircraft) end

"""
    get_chord(ac::Aircraft)

Get aircraft mean aerodynamic chord, ``c``, (m).
"""
function get_chord(ac::Aircraft) end

"""
    get_wing_span(ac::Aircraft)

Get aircraft wing span, ``b``, (m).
"""
function get_wing_span(ac::Aircraft) end

"""
    get_aspect_ratio(ac::Aircraft)

Get aspect ratio, ``b^2/S``.
"""
function get_aspect_ratio(ac::Aircraft)
    b = get_wing_span(ac)
    S = get_surface(ac)
    return b^2 / S
end

# Forces and moments
# TODO: document once headers are standarized. Right now they are too specific to F16Stevens
function calculate_prop_forces_moments end
function calculate_prop_gyro_effects end
function calculate_pdot end

function calculate_aero_forces_moments end

# Flight Control system
function tgear end


text:
using LearnConvection
using Plots

D=16
N=8

problem  = Sequential("T")

train = ["general_strat_4_profiles.jld2", "general_strat_32_profiles.jld2"]
validate = ["general_strat_12_profiles.jld2", "general_strat_24_profiles.jld2"]
test = ["general_strat_8_profiles.jld2", "general_strat_16_profiles.jld2", "general_strat_20_profiles.jld2", "general_strat_28_profiles.jld2"]

𝒟_train = LearnConvection.Data.data(train, problem; D=D, N=N)
𝒟_validate  = LearnConvection.Data.data(validate, problem; D=D, N=N)
𝒟_test = LearnConvection.Data.data(test, problem; D=D, N=N)

train_validate_test(𝒟_train, 𝒟_validate, 𝒟_test, problem; log_γs=-0.4:0.1:0.4)

𝒢 = LearnConvection.GaussianProcess.model(𝒟_train; kernel=get_kernel(2,0.0,0.0,euclidean_distance))
predict(𝒢, 𝒟_train; postprocessed=true)
LearnConvection.GaussianProcess.get_me_true_check(𝒢, 𝒟_validate)

k=5
distance=euclidean_distance
get_min_gamma(k, distance, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=-0.4:0.1:0.4)
get_min_gamma_alpha(5, distance, 𝒟_train, 𝒟_validate, 𝒟_test; log_γs=-0.4:0.1:0.4)

##

using LearnConvection
using Plots

# Construct a ProfileData object, 𝒟, consisting of the data to train on,
# and a GP object, 𝒢, with which to perform the regression.

# data
filename = "general_strat_32_profiles.jld2"
problem  = Residual("KPP", KPP.Parameters()) # v: "T" for temperature profile; "wT" for temperature flux profile
problem  = Residual("TKE", TKEMassFlux.TKEParameters())
D        = 32       # collapse profile data down to 16 gridpoints
N        = 4        # collect every 4 timesteps' data for training

# Now let's define our kernel. We don't know what the best long(length-scale) parameter value is yet so let's guess 0.0 (γ=1) and see what happens.
# Let's use an exponential kernel.

# kernel
k        = 3        # kernel function ID
logγ     = -0.6      # log(length-scale parameter)
logσ     = 0.0      # log(signal variance parameter)
distance = derivative_distance  # distance metric to use in the kernel
kernel   = get_kernel(k, logγ, logσ, distance)

# data
𝒟_train = LearnConvection.Data.data(["general_strat_16_profiles.jld2","general_strat_32_profiles.jld2"], problem; D=D, N=N)
𝒟_test = LearnConvection.Data.data("general_strat_32_profiles.jld2", problem; D=D, N=N)

# model
𝒢 = LearnConvection.GaussianProcess.model(𝒟_train; kernel = kernel)

# Animate the mean GP prediction.
# anim = animate_profile(𝒢, 𝒟)
# gif(anim, "animated_profile.gif")

anim = animate_profile_and_model_output(𝒢, 𝒟_test)
gif(anim, "animated_profile_and_model_output_KPP32_train_gs16_gs32_test_gs32.gif")

##

# Not great. Let's try to optimize our γ value.
# We can use get_min_gamma function to search the range -2:0.1:2 for the log(γ) value that minimizes the mean error on the true check
min_log_param, min_error = get_min_gamma(k, 𝒟, distance, -2:0.01:2)

# Let's create a kernel with the better log(γ) value, put the new kernel into a new model, and see how it does.
new_kernel = get_kernel(k, min_log_param, logσ, distance)
𝒢 = LearnConvection.GaussianProcess.model(𝒟; kernel = new_kernel)
anim = animate_profile(𝒢, 𝒟)
gif(anim, "animated_profile_optimized.gif")

# There's a lot of flexibility in how we design our kernel function.
# We can modify the hyperparameter values, the distance metric, and the type of kernel function.

p1 = plot([1 2 3], [4 5 6], title="hello", ylabel="yes")
p2 = plot([1 2 3], [4 5 6])
p3 = plot([1 2 3], [4 5 6])
p4 = plot([1 2 3], [4 5 6])
p5 = plot([1 2 3], [4 5 6])
p6 = plot([1 2 3], [4 5 6])

layout = @layout [a b c; d e f]
titles = ["Sq Exp / l²-norm" "Matérn 1/2 / l²-norm" "Matérn 3/2 / l²-norm" "Sq Exp / H¹-norm" "Sq Exp / H¹-norm" "Sq Exp / H¹-norm"]

plot(p1, p2, p3, p4, p5, p6; layout = layout, title=titles)

# data
𝒟 = OceanConvect.ModelData.data("general_strat_16_profiles.jld2", Sequential("dT"); D=16, N=4)
logγs = -3.0:0.1:3.0
p1  = plot_landscapes_compare_error_metrics(1, 𝒟, euclidean_distance, log_γs)
p2  = plot_landscapes_compare_error_metrics(2, 𝒟, euclidean_distance, log_γs)
p3  = plot_landscapes_compare_error_metrics(3, 𝒟, euclidean_distance, log_γs)
p4  = plot_landscapes_compare_error_metrics(4, 𝒟, euclidean_distance, log_γs)
p5  = plot_landscapes_compare_error_metrics(5, 𝒟, euclidean_distance, log_γs)
p6  = plot_landscapes_compare_error_metrics(1, 𝒟, derivative_distance, log_γs)
p7  = plot_landscapes_compare_error_metrics(2, 𝒟, derivative_distance, log_γs)
p8  = plot_landscapes_compare_error_metrics(3, 𝒟, derivative_distance, log_γs)
p9  = plot_landscapes_compare_error_metrics(4, 𝒟, derivative_distance, log_γs)
p10 = plot_landscapes_compare_error_metrics(5, 𝒟, derivative_distance, log_γs)
p11 = plot_landscapes_compare_error_metrics(1, 𝒟, antiderivative_distance, log_γs)
p12 = plot_landscapes_compare_error_metrics(2, 𝒟, antiderivative_distance, log_γs)
p13 = plot_landscapes_compare_error_metrics(3, 𝒟, antiderivative_distance, log_γs)
p14 = plot_landscapes_compare_error_metrics(4, 𝒟, antiderivative_distance, log_γs)
p15 = plot_landscapes_compare_error_metrics(5, 𝒟, antiderivative_distance, log_γs)

##

# problems
params   = KPP.Parameters( CSL = 1.0, CNL = 1.0, Cb_T = 1.0, CKE = 1.0)
params   = KPP.Parameters( )

problem  = Residual("KPP", params) # v: "T" for temperature profile; "wT" for temperature flux profile

problem  = Sequential("dT") # v: "T" for temperature profile; "wT" for temperature flux profile

# data
filename = "general_strat_32_profiles.jld2"
D        = 16
N        = 4

# kernel
k        = 1
logγ     = -3.0
logσ     = 0.0
distance = euclidean_distance
kernel   = get_kernel(k, logγ, logσ, distance)

# data
𝒟 = LearnConvection.Data.data(filename, problem; D=D, N=N)

# model
𝒢 = LearnConvection.GaussianProcess.model(𝒟; kernel = kernel)

# Animate the mean GP prediction.
anim = animate_profile(𝒢, 𝒟)
gif(anim, "animated_profile.gif")

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using Bridge, ConjugatePriors, Plots, Distributions, LaTeXStrings

srand(5);
pgfplots()

simname = :gammachain
simid = 1

#joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)paths.jld2")

Nf = 800001 # divisible by 25, 50, 100
step = [25, 50, 100, 200][4] # subsampling stepsize
H = [20,40,80,160,320][4]

skipp = 100

s1(t) = 3/2 + sin(2*(4t-2)) + 2exp(-16((4t-2)).^2) # called s_1 in the paper
b(x) = -10x + 20. # called b_1 in the paper

ss = [0.1,0.13,0.15,0.23,0.25,0.4,0.44,0.65,0.76,0.78,0.81]
h = [4,-5,3,-4,5,-4.2,2.1,4.3,-3.1,2.1,-4.2]
s2(t) = sum(h[])

s = s1
dtf = 1/(Nf-1)
ttf = 0:dtf:1

N = (Nf-1) ÷ step + 1 # number of observations: n
td = (0:dtf:1)[1:step:end]
dt = td[2] - td[1]

M = (length(td)-1) ÷ H # observed intervals per bin
pri = InverseGamma(0.1, 0.1)

r = length(td)-1 - H*M
assert(H*M == length(td)-1)

a1 = 1/1000
a, aζ = 1000.0/sqrt(H), 1000.0/sqrt(H)

import Bridge.piecewise

function piecewise(Y::SamplePath, tend = Y.tt[end])
    tt = [Y.tt[1]]
    n = length(Y.yy)
    append!(tt, repeat(Y.tt[2:n], inner=2))
    push!(tt, tend)
    tt, repeat(Y.yy, inner=2)
end
function piecewise(X::Bridge.EstSamplePath, tend = X.tt[end])
    Bridge.EstSamplePath(Bridge.piecewise(SamplePath(X.tt, X.yy), tend)...,
        Bridge.piecewise(SamplePath(X.tt, X.vv), tend)[2])
end

@recipe function plot(X::Bridge.EstSamplePath{Float64}; quantile = 0.95)
    q = Distributions.quantile(Normal(), 1 - (1 - quantile)/2)
    if Plots.backend() != Plots.PGFPlotsBackend()
        ribbon --> q*sqrt.(X.vv)
        X.tt, X.yy
    else
        ribbon --> (NaN,NaN)
        linewidth --> 0
        vcat(X.tt, reverse(X.tt)), vcat(X.yy + q*sqrt.(X.vv), reverse(X.yy - q*sqrt.(X.vv)))
    end
end

@recipe function plot(X::SamplePathBand)
    if Plots.backend() != Plots.PGFPlotsBackend()
        error("todo")
    else
        ribbon --> (NaN,NaN)
        linewidth --> 0
        vcat(X.tt, reverse(X.tt)), vcat(X.uu, reverse(X.ll))
    end
end

"""
White noise on the grid tt
"""
dW(tt) = randn(length(tt)-1).*sqrt.(diff(tt))

"""
Finite difference operator
"""
function diff0(x)
    y = similar(x)
    y[1] = zero(x)
    for i in 2:length(x)
        y[i] = x[i]-x[i-1]
    end
    y
end

"""
Euler scheme with start in y1, space dependent drift b and time dependent dispersion coefficient
"""
function euler1(tt, y1, b, sigma, yy = zeros(tt), dW = dW(tt))
    yy[1] = y1
    for i in 2:length(yy)
        yy[i] = yy[i-1] + b(yy[i-1]) * (tt[i]-tt[i-1]) + sigma(tt[i-1])*dW[i-1]
    end
    yy
end

Xf = euler1(ttf, 0., b, s, zeros(ttf), dW(ttf))
;

Xd = Xf[1:step:end]

ttplot = ttf[1:skipp:end]
Xplot = Xf[1:skipp:end]
plot(ttplot, Xplot, label=L"X_t")

post = Vector(H)
Z = Vector(H)

for i in 1:H
    ii = 1+(i-1)*M:(i)*M+1
    tx = td[ii]
    x = diff(Xd[ii])./sqrt.(diff(tx)) # normalize obs
    Z[i] = sum(diff(Xd[ii]).^2)
    post[i] = ConjugatePriors.posterior((0., pri), Normal, x)
end

@show rate(post[end])
@show rate(InverseGamma(0.1 + M/2, 0.1 + Z[end]/2/dt))

xhat = sqrt.(median.(post));

plot(Bridge.piecewise(SamplePath(td[1:M:end],xhat)))
plot!(ttplot, s.(ttplot))

v = zeros(H)
v[:] = xhat
#ζ = zeros(H + 1);
ζ = zeros(H);

v[1] = rand(InverseGamma(a1, a1))

for i in 2:H
    ζ[i] = rand(InverseGamma(aζ, inv(v[i-1]/aζ)))
    v[i] = rand(InverseGamma(a, inv(ζ[i]/a)))
end

plot(v[1:end])

iterations = 10000
state = Bridge.mcstart(v)
for iter in 1:iterations
    # sample z
    ζ[1] = rand(InverseGamma(aζ + a + 1, aζ/bζ + a/v[1]))
    for h in 2:H
        ζ[h] = rand(InverseGamma(aζ + a, (aζ/v[h-1] + a/v[h])))
    end
  #  ζ[1] = 1/rand(Gamma(aζ + a, aζ/v[H] + a/v[1]))
    for h in 1:H-1
        v[h] = rand(InverseGamma(aζ + a + M/2, (a/ζ[h] + aζ/ζ[h+1] + Z[h]/2/dt)))
    end
    v[H] = rand(InverseGamma(aζ + (M+r)/2, a/ζ[H] + Z[H]/2/dt))
    state = Bridge.mcnext(state, v)
end
v;

plot(Bridge.EstSamplePath(td[1:M:end-1], Bridge.mcstats(state)...), size = (900,300),
    yaxis =((0,20),0:5:20), color=colorant"#0044FF", label = "marginal 95\\% cov.~band", layout=2, subplot=1)
plot!(ttplot, (s.(ttplot)).^2, linewidth=2, color="red", label=L"s_2", subplot=1)
plot!(Bridge.EstSamplePath(td[1:M:end-1], mean.(post), var.(post) ),
    yaxis = ((0,20),0:5:20), color=colorant"#0044FF", label = "marginal 95\\% cov.~band", layout=2, subplot=2)
plot!(ttplot, (s.(ttplot)).^2, linewidth=2, color="red", label = L"s_2", subplot=2)

#savefig("chain.pdf")

plot(piecewise( Bridge.EstSamplePath(td[1:M:end], Bridge.mcstats(state)...) ), size = (900,300),
    yaxis =((0,20),0:5:20), color=colorant"#0044FF", label = "marginal 95\\% cov.~band", layout=2, subplot=1)
plot!(ttplot, (s.(ttplot)).^2, linewidth=2, color="red", label=L"s_2", subplot=1)
plot!(piecewise( Bridge.EstSamplePath(td[1:M:end], mean.(post), var.(post)) ),
    yaxis = ((0,20),0:5:20), color=colorant"#0044FF", label = "marginal 95\\% cov.~band", layout=2, subplot=2)
plot!(ttplot, (s.(ttplot)).^2, linewidth=2, color="red", label = L"s_2", subplot=2)

diff(td[1:M:M+1]),
1/H

#length(td[1:M:end-1]), length(Bridge.mcstats(state)[1])

#plot(piecewise(Bridge.EstSamplePath(td[1:M:end], Bridge.mcstats(state)...)))


text:
module ToeplitzMatrices
using StatsBase

import Base: convert, *, \, getindex, print_matrix, size, Matrix, +, -, copy, similar, sqrt, copyto!,
    adjoint, transpose
import LinearAlgebra: BlasReal, Cholesky, DimensionMismatch, cholesky, cholesky!, eigvals, inv, ldiv!,
    mul!, pinv, rmul!, tril, triu

using FFTW
using FFTW: Plan
flipdim(A, d) = reverse(A, dims=d)

export Toeplitz, SymmetricToeplitz, Circulant, TriangularToeplitz, Hankel,
       chan, strang

include("iterativeLinearSolvers.jl")

# Abstract
abstract type AbstractToeplitz{T<:Number} <: AbstractMatrix{T} end

size(A::AbstractToeplitz) = (size(A, 1), size(A, 2))
function getindex(A::AbstractToeplitz, i::Integer)
    return A[mod(i - 1, size(A, 1)) + 1, div(i - 1, size(A, 1)) + 1]
end

convert(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S)
convert(::Type{AbstractArray{T}}, S::AbstractToeplitz) where {T} = convert(AbstractToeplitz{T}, S)

# Convert an abstract Toeplitz matrix to a full matrix
function Matrix(A::AbstractToeplitz{T}) where T
    m, n = size(A)
    Af = Matrix{T}(undef, m, n)
    for j = 1:n
        for i = 1:m
            Af[i,j] = A[i,j]
        end
    end
    return Af
end

convert(::Type{Matrix}, A::AbstractToeplitz) = Matrix(A)

# Fast application of a general Toeplitz matrix to a column vector via FFT
function mul!(y::StridedVector{T}, A::AbstractToeplitz{T}, x::StridedVector, α::T, β::T) where T
    m = size(A,1)
    n = size(A,2)
    N = length(A.vcvr_dft)
    if m != length(y)
        throw(DimensionMismatch(""))
    end
    if n != length(x)
        throw(DimensionMismatch(""))
    end

    # In any case, scale/initialize y
    if iszero(β)
        fill!(y, 0)
    else
        rmul!(y, β)
    end

    @inbounds begin
        # Small case: don't use FFT
        if N < 512
            for j in 1:n
                tmp = α * x[j]
                for i in 1:m
                    y[i] += tmp*A[i,j]
                end
            end
            return y
        end

        # Large case: use FFT
        for i in 1:n
            A.tmp[i] = x[i]
        end
        for i in n+1:N
            A.tmp[i] = 0
        end
        mul!(A.tmp, A.dft, A.tmp)
        for i = 1:N
            A.tmp[i] *= A.vcvr_dft[i]
        end
        A.dft \ A.tmp
        for i in 1:m
            y[i] += α * (T <: Real ? real(A.tmp[i]) : A.tmp[i])
        end
        return y
    end
end
# Avoid ambiguity error
mul!(y::StridedVector{T}, A::AbstractToeplitz{T}, x::StridedVector, α::T, β::T) where {T<:Number} =
    invoke(mul!, Tuple{StridedVector{T},ToeplitzMatrices.AbstractToeplitz{T},StridedVector,T,T} where T,
        y, A, x, α, β)

# Application of a general Toeplitz matrix to a general matrix
function mul!(C::StridedMatrix{T}, A::AbstractToeplitz{T}, B::StridedMatrix, α::T, β::T) where T
    l = size(B, 2)
    if size(C, 2) != l
        throw(DimensionMismatch("input and output matrices must have same number of columns"))
    end
    for j = 1:l
        mul!(view(C, :, j), A, view(B, :, j), α, β)
    end
    return C
end
# Avoid ambiguity error
mul!(C::StridedMatrix{T}, A::AbstractToeplitz{T}, B::StridedMatrix, α::T, β::T) where {T<:Number} =
    invoke(mul!, Tuple{StridedMatrix{T},ToeplitzMatrices.AbstractToeplitz{T},StridedMatrix,T,T} where T,
        C, A, B, α, β)

# Translate three to five argument mul!
mul!(y::StridedVecOrMat, A::AbstractToeplitz, x::StridedVecOrMat) =
    mul!(y, A, x, one(eltype(A)), zero(eltype(A)))

# Left division of a general matrix B by a general Toeplitz matrix A, i.e. the solution x of Ax=B.
function ldiv!(A::AbstractToeplitz, B::StridedMatrix)
    if size(A, 1) != size(A, 2)
        error("Division: Rectangular case is not supported.")
    end
    for j = 1:size(B, 2)
        ldiv!(A, view(B, :, j))
    end
    return B
end

function (\)(A::AbstractToeplitz, b::AbstractVector)
    T = promote_type(eltype(A), eltype(b))
    if T != eltype(A)
        throw(ArgumentError("promotion of Toeplitz matrices not handled yet"))
    end
    bb = similar(b, T)
    copyto!(bb, b)
    ldiv!(A, bb)
end

# General Toeplitz matrix
mutable struct Toeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
    vc::Vector{T}
    vr::Vector{T}
    vcvr_dft::Vector{S}
    tmp::Vector{S}
    dft::Plan{S}
end

# Ctor
function Toeplitz{T}(vc::Vector, vr::Vector) where {T}
    m, n = length(vc), length(vr)
    if vc[1] != vr[1]
        error("First element of the vectors must be the same")
    end

    vcp, vrp = Vector{T}(vc), Vector{T}(vr)

    tmp = Vector{promote_type(T, Complex{Float32})}(undef, m + n - 1)
    copyto!(tmp, vcp)
    for i = 1:n - 1
        tmp[i + m] = vrp[n - i + 1]
    end
    dft = plan_fft!(tmp)
    return Toeplitz(vcp, vrp, dft*tmp, similar(tmp), dft)
end

Toeplitz(vc::Vector, vr::Vector) =
    Toeplitz{promote_type(eltype(vc), eltype(vr), Float32)}(vc, vr)

# Toeplitz(A::AbstractMatrix) projects onto Toeplitz part using the first row/col
Toeplitz(A::AbstractMatrix) = Toeplitz(A[:,1], A[1,:])
Toeplitz{T}(A::AbstractMatrix) where T = Toeplitz{T}(A[:,1], A[1,:])

convert(::Type{AbstractToeplitz{T}}, A::Toeplitz) where {T} = convert(Toeplitz{T}, A)
convert(::Type{Toeplitz{T}}, A::Toeplitz) where {T} = Toeplitz(convert(Vector{T}, A.vc),
                                                               convert(Vector{T}, A.vr))

adjoint(A::Toeplitz) = Toeplitz(conj(A.vr), conj(A.vc))
adjoint(A::Toeplitz{<:Real}) = transpose(A)
transpose(A::Toeplitz) = Toeplitz(A.vr, A.vc)

# Size of a general Toeplitz matrix
function size(A::Toeplitz, dim::Int)
    if dim == 1
        return length(A.vc)
    elseif dim == 2
        return length(A.vr)
    elseif dim > 2
        return 1
    else
        error("arraysize: dimension out of range")
    end
end

# Retrieve an entry
function getindex(A::Toeplitz, i::Integer, j::Integer)
    m = size(A,1)
    n = size(A,2)
    if i > m || j > n
        error("BoundsError()")
    end

    if i >= j
        return A.vc[i - j + 1]
    else
        return A.vr[1 - i + j]
    end
end

# Form a lower triangular Toeplitz matrix by annihilating all entries above the k-th diaganal
function tril(A::Toeplitz, k = 0)
    if k > 0
        error("Second argument cannot be positive")
    end
    Al = TriangularToeplitz(copy(A.vc), 'L', length(A.vr))
    if k < 0
      for i in -1:-1:k
          Al.ve[-i] = zero(eltype(A))
      end
    end
    return Al
end

# Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diaganal
function triu(A::Toeplitz, k = 0)
    if k < 0
        error("Second argument cannot be negative")
    end
    Al = TriangularToeplitz(copy(A.vr), 'U', length(A.vc))
    if k > 0
      for i in 1:k
          Al.ve[i] = zero(eltype(A))
      end
    end
    return Al
end

ldiv!(A::Toeplitz, b::StridedVector) =
    copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, strang(A), 1000, 100eps())[1])

# Symmetric
mutable struct SymmetricToeplitz{T<:BlasReal} <: AbstractToeplitz{T}
    vc::Vector{T}
    vcvr_dft::Vector{Complex{T}}
    tmp::Vector{Complex{T}}
    dft::Plan
end

function SymmetricToeplitz{T}(vc::Vector{T}) where T<:BlasReal
	tmp = convert(Array{Complex{T}}, [vc; zero(T); reverse(vc[2:end])])
	dft = plan_fft!(tmp)
	return SymmetricToeplitz{T}(vc, dft*tmp, similar(tmp), dft)
end

SymmetricToeplitz{T}(vc::AbstractVector) where T<:BlasReal = SymmetricToeplitz{T}(convert(Vector{T}, vc))
SymmetricToeplitz{T}(vc::AbstractVector{T}) where T = SymmetricToeplitz{promote_type(Float32, T)}(vc)
SymmetricToeplitz{T}(vc::AbstractVector) where T = SymmetricToeplitz{T}(convert(Vector{T}, vc))
SymmetricToeplitz(vc::AbstractVector{T}) where T = SymmetricToeplitz{T}(vc)

SymmetricToeplitz{T}(A::AbstractMatrix) where T = SymmetricToeplitz{T}(A[1, :])
SymmetricToeplitz(A::AbstractMatrix) = SymmetricToeplitz{eltype(A)}(A)

convert(::Type{AbstractToeplitz{T}}, A::SymmetricToeplitz) where {T} = convert(SymmetricToeplitz{T},A)
convert(::Type{SymmetricToeplitz{T}}, A::SymmetricToeplitz) where {T} = SymmetricToeplitz(convert(Vector{T},A.vc))

adjoint(A::SymmetricToeplitz) = SymmetricToeplitz(conj(A.vr), conj(A.vc))
adjoint(A::SymmetricToeplitz{<:Real}) = A
transpose(A::SymmetricToeplitz) = A

function size(A::SymmetricToeplitz, dim::Int)
    if 1 <= dim <= 2
        return length(A.vc)
    else
        error("arraysize: dimension out of range")
    end
end

getindex(A::SymmetricToeplitz, i::Integer, j::Integer) = A.vc[abs(i - j) + 1]

ldiv!(A::SymmetricToeplitz, b::StridedVector) =
    copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])

# Circulant
mutable struct Circulant{T<:Number,S<:Number} <: AbstractToeplitz{T}
    vc::Vector{T}
    vcvr_dft::Vector{S}
    tmp::Vector{S}
    dft::Plan
end

function Circulant{T}(vc::Vector{T}) where T
    tmp = zeros(promote_type(T, Complex{Float32}), length(vc))
    return Circulant(vc, fft(vc), tmp, plan_fft!(tmp))
end

Circulant{T}(vc::AbstractVector) where T = Circulant{T}(convert(Vector{T}, vc))
Circulant(vc::AbstractVector) = Circulant{promote_type(eltype(vc), Float32)}(vc)
Circulant{T}(A::AbstractMatrix) where T = Circulant{T}(A[:,1])
Circulant(A::AbstractMatrix) = Circulant(A[:,1])

convert(::Type{AbstractToeplitz{T}}, A::Circulant) where {T} = convert(Circulant{T}, A)
convert(::Type{Circulant{T}}, A::Circulant) where {T} = Circulant(convert(Vector{T}, A.vc))

function size(C::Circulant, dim::Integer)
    if 1 <= dim <= 2
        return length(C.vc)
    else
        error("arraysize: dimension out of range")
    end
end

function getindex(C::Circulant, i::Integer, j::Integer)
    n = size(C, 1)
    if i > n || j > n
        error("BoundsError()")
    end
    return C.vc[mod(i - j, length(C.vc)) + 1]
end

function Ac_mul_B(A::Circulant{T}, B::Circulant{T}) where T<:Real
    tmp = similar(A.vcvr_dft)
    for i = 1:length(tmp)
        tmp[i] = conj(A.vcvr_dft[i]) * B.vcvr_dft[i]
    end
    c_tmp = copy(tmp)
    return Circulant(real(A.dft \ tmp), c_tmp, A.tmp, A.dft)
end
function Ac_mul_B(A::Circulant, B::Circulant)
    T = promote_type(eltype(A), eltype(B))
    tmp = similar(A.vcvr_dft, T)
    for i = 1:length(tmp)
        tmp[i] = conj(A.vcvr_dft[i]) * B.vcvr_dft[i]
    end
    c_tmp = copy(tmp)
    tmp2 = A.dft \ tmp
    return Circulant(Vector{T}(tmp2), c_tmp, Vector{T}(A.tmp), eltype(A) == T ? A.dft : plan_fft!(tmp2))
end

function ldiv!(C::Circulant{T}, b::AbstractVector{T}) where T
    n = length(b)
    size(C, 1) == n || throw(DimensionMismatch(""))
    for i = 1:n
        C.tmp[i] = b[i]
    end
    C.dft * C.tmp
    for i = 1:n
        C.tmp[i] /= C.vcvr_dft[i]
    end
    C.dft \ C.tmp
    for i = 1:n
        b[i] = (T <: Real ? real(C.tmp[i]) : C.tmp[i])
    end
    return b
end

function inv(C::Circulant{T}) where T<:Real
    vdft = 1 ./ C.vcvr_dft
    # Copy the vector because the dft below will overwrite it
    c_vdft = copy(vdft)
    return Circulant(real(C.dft \ vdft), c_vdft, similar(vdft), C.dft)
end
function inv(C::Circulant)
    vdft = 1 ./ C.vcvr_dft
    c_vdft = copy(vdft)
    return Circulant(C.dft \ vdft, c_vdft, similar(vdft), C.dft)
end

function strang(A::AbstractMatrix{T}) where T
    n = size(A, 1)
    v = Vector{T}(undef, n)
    n2 = div(n, 2)
    for i = 1:n
        if i <= n2 + 1
            v[i] = A[i,1]
        else
            v[i] = A[1, n - i + 2]
        end
    end
    return Circulant(v)
end
function chan(A::AbstractMatrix{T}) where T
    n = size(A, 1)
    v = Vector{T}(undef, n)
    for i = 1:n
        v[i] = ((n - i + 1) * A[i, 1] + (i - 1) * A[1, min(n - i + 2, n)]) / n
    end
    return Circulant(v)
end

function pinv(C::Circulant{T}, tolerance::T = eps(T)) where T<:Real
    vdft = copy(C.vcvr_dft)
    vdft[abs.(vdft).<tolerance] .= Inf
    vdft .= 1 ./ vdft
    c_vdft = copy(vdft)
    return Circulant(real(C.dft \ vdft), c_vdft, similar(vdft), C.dft)
end

function pinv(C::Circulant{T}, tolerance::Real = eps(real(T))) where T<:Number
    vdft = copy(C.vcvr_dft)
    vdft[abs.(vdft).<tolerance] .= Inf
    vdft .= 1 ./ vdft
    c_vdft = copy(vdft)
    return Circulant(C.dft \ vdft, c_vdft, similar(vdft), C.dft)
end

eigvals(C::Circulant) = copy(C.vcvr_dft)
sqrt(C::Circulant{T}) where T<:Real = Circulant(real(ifft(sqrt.(C.vcvr_dft))))
sqrt(C::Circulant) = Circulant(ifft(sqrt.(C.vcvr_dft)))
copy(C::Circulant) = Circulant(copy(C.vc))
similar(C::Circulant) = Circulant(similar(C.vc))
function copyto!(dest::Circulant{U,S}, src::Circulant{V,S}) where {U,V,S}
    copyto!(dest.vc, src.vc)
    copyto!(dest.vcvr_dft, src.vcvr_dft)
end

function (+)(C1::Circulant, C2::Circulant)
    @boundscheck (size(C1)==size(C2)) || throw(BoundsError())
    Circulant(C1.vc+C2.vc)
end

function (-)(C1::Circulant, C2::Circulant)
    @boundscheck (size(C1)==size(C2)) || throw(BoundsError())
    Circulant(C1.vc-C2.vc)
end

(-)(C::Circulant) = Circulant(-C.vc)

function (*)(C1::Circulant{<:Real}, C2::Circulant{<:Real})
    @boundscheck (size(C1)==size(C2)) || throw(BoundsError())
    Circulant(real(ifft(C1.vcvr_dft.*C2.vcvr_dft)))
end

function (*)(C1::Circulant, C2::Circulant)
    @boundscheck (size(C1)==size(C2)) || throw(BoundsError())
    Circulant(ifft(C1.vcvr_dft.*C2.vcvr_dft))
end

(*)(scalar::Number, C::Circulant) = Circulant(scalar*C.vc)
(*)(C::Circulant,scalar::Number) = Circulant(scalar*C.vc)

# Triangular
mutable struct TriangularToeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
    ve::Vector{T}
    uplo::Char
    vcvr_dft::Vector{S}
    tmp::Vector{S}
    dft::Plan
end

function TriangularToeplitz{T}(vep::Vector{T}, uplo::Symbol) where T
    n = length(vep)

    tmp = zeros(promote_type(T, Complex{Float32}), 2n - 1)
    if uplo == :L
        copyto!(tmp, vep)
    else
        tmp[1] = vep[1]
        for i = 1:n - 1
            tmp[n + i] = vep[n - i + 1]
        end
    end
    dft = plan_fft!(tmp)
    return TriangularToeplitz(vep, string(uplo)[1], dft * tmp, similar(tmp), dft)
end

TriangularToeplitz{T}(ve::AbstractVector, uplo::Symbol) where T =
    TriangularToeplitz(convert(Vector{T}, ve), uplo)

TriangularToeplitz(ve::AbstractVector, uplo::Symbol) =
    TriangularToeplitz{promote_type(eltype(ve), Float32)}(ve, uplo)

TriangularToeplitz{T}(A::AbstractMatrix, uplo::Symbol) where T =
    TriangularToeplitz{T}(uplo == :U ? A[1,:] : A[:,1], uplo)

TriangularToeplitz(A::AbstractMatrix, uplo::Symbol) =
    TriangularToeplitz(uplo == :U ? A[1,:] : A[:,1], uplo)

function convert(::Type{Toeplitz}, A::TriangularToeplitz)
    if A.uplo == 'L'
        Toeplitz(A.ve, [A.ve[1]; zeros(length(A.ve) - 1)])
    else
        @assert A.uplo == 'U'
        Toeplitz([A.ve[1]; zeros(length(A.ve) - 1)], A.ve)
    end
end

convert(::Type{AbstractToeplitz{T}}, A::TriangularToeplitz) where {T} = convert(TriangularToeplitz{T},A)
convert(::Type{TriangularToeplitz{T}}, A::TriangularToeplitz) where {T} =
    TriangularToeplitz(convert(Vector{T},A.ve),A.uplo=='U' ? (:U) : (:L))

function size(A::TriangularToeplitz, dim::Int)
    if dim == 1 || dim == 2
        return length(A.ve)
    elseif dim > 2
        return 1
    else
        error("arraysize: dimension out of range")
    end
end

function getindex(A::TriangularToeplitz{T}, i::Integer, j::Integer) where T
    if A.uplo == 'L'
        return i >= j ? A.ve[i - j + 1] : zero(T)
    else
        return i <= j ? A.ve[j - i + 1] : zero(T)
    end
end

function (*)(A::TriangularToeplitz, B::TriangularToeplitz)
    n = size(A, 1)
    if n != size(B, 1)
        throw(DimensionMismatch(""))
    end
    if A.uplo == B.uplo
        return TriangularToeplitz(conv(A.ve, B.ve)[1:n], A.uplo)
    end
    return Triangular(Matrix(A), A.uplo) * Triangular(Matrix(B), B.uplo)
end

Ac_mul_B(A::TriangularToeplitz, b::AbstractVector) =
    TriangularToeplitz(A.ve, A.uplo == 'U' ? :L : :U) * b

# NB! only valid for lower triangular
function smallinv(A::TriangularToeplitz{T}) where T
    n = size(A, 1)
    b = zeros(T, n)
    b[1] = 1 ./ A.ve[1]
    for k = 2:n
        tmp = zero(T)
        for i = 1:k-1
            tmp += A.uplo == 'L' ? A.ve[k - i + 1]*b[i] : A.ve[i + 1] * b[k - i]
        end
        b[k] = -tmp/A.ve[1]
    end
    return TriangularToeplitz(b, symbol(A.uplo))
end

function inv(A::TriangularToeplitz{T}) where T
    n = size(A, 1)
    if n <= 64
        return smallinv(A)
    end
    np2 = nextpow2(n)
    if n != np2
        return TriangularToeplitz(inv(TriangularToeplitz([A.ve, zeros(T, np2 - n)],
            symbol(A.uplo))).ve[1:n], symbol(A.uplo))
    end
    nd2 = div(n, 2)
    a1 = inv(TriangularToeplitz(A.ve[1:nd2], symbol(A.uplo))).ve
    return TriangularToeplitz([a1, -(TriangularToeplitz(a1, symbol(A.uplo)) *
        (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))], symbol(A.uplo))
end

# ldiv!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b
ldiv!(A::TriangularToeplitz, b::StridedVector) =
    copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, chan(A), 1000, 100eps())[1])

# extend levinson
StatsBase.levinson!(x::StridedVector, A::SymmetricToeplitz, b::StridedVector) =
    StatsBase.levinson!(A.vc, b, x)
function StatsBase.levinson!(C::StridedMatrix, A::SymmetricToeplitz, B::StridedMatrix)
    n = size(B, 2)
    if n != size(C, 2)
        throw(DimensionMismatch("input and output matrices must have same number of columns"))
    end
    for j = 1:n
        StatsBase.levinson!(view(C, :, j), A, view(B, :, j))
    end
    C
end
StatsBase.levinson(A::AbstractToeplitz, B::StridedVecOrMat) =
    StatsBase.levinson!(zeros(size(B)), A, copy(B))

# BlockTriangular
# type BlockTriangularToeplitz{T<:BlasReal} <: AbstractMatrix{T}
#     Mc::Array{T,3}
#     uplo::Char
#     Mc_dft::Array{Complex{T},3}
#     tmp::Vector{Complex{T}}
#     dft::Plan
# end
# function BlockTriangularToeplitz{T<:BlasReal}(Mc::Array{T,3}, uplo::Symbol)
#     n, p, _ = size(Mc)
#     tmp = zeros(Complex{T}, 2n)
#     dft = plan_fft!(tmp)
#     Mc_dft = Array{Complex{T}}(2n, p, p)
#     for j = 1:p
#         for i = 1:p
#             Mc_dft[1,i,j] = complex(Mc[1,i,j])
#             for t = 2:n
#                 Mc_dft[t,i,j] = uplo == :L ? complex(Mc[t,i,j]) : zero(Complex{T})
#             end
#             Mc_dft[n+1,i,j] = zero(Complex{T})
#             for t = n+2:2n
#                 Mc_dft[t,i,j] = uplo == :L ? zero(Complex{T}) : complex(Mc[2n-t+2,i,j])
#             end
#             dft(view(Mc_dft, 1:2n, 1:p, 1:p))
#         end
#     end
#     return BlockTriangularToeplitz(Mc, string(uplo)[1], Mc_dft, tmp, dft, idft)
# end

#= Hankel Matrix
 A Hankel matrix is a matrix that is constant across the anti-diagonals:

  [a_0 a_1 a_2 a_3 a_4
   a_1 a_2 a_3 a_4 a_5
   a_2 a_3 a_4 a_5 a_6]

 This is precisely a Toeplitz matrix with the columns reversed:
                             [0 0 0 0 1
  [a_4 a_3 a_2 a_1 a_0        0 0 0 1 0
   a_5 a_4 a_3 a_2 a_1   *    0 0 1 0 0
   a_6 a_5 a_4 a_3 a_2]       0 1 0 0 0
                              1 0 0 0 0]
 We represent the Hankel matrix by wrapping the corresponding Toeplitz matrix.=#

# Hankel Matrix, use _Hankel as Hankel(::Toeplitz) should project to Hankel
function _Hankel end

# Hankel Matrix
mutable struct Hankel{TT<:Number} <: AbstractMatrix{TT}
    T::Toeplitz{TT}
    global _Hankel(T::Toeplitz{TT}) where TT<:Number = new{TT}(T)
end

# Ctor: vc is the leftmost column and vr is the bottom row.
function Hankel{T}(vc::AbstractVector, vr::AbstractVector) where T
    if vc[end] != vr[1]
        error("First element of rows must equal last element of columns")
    end
    n = length(vr)
    p = [vc; vr[2:end]]
    _Hankel(Toeplitz{T}(p[n:end],p[n:-1:1]))
end

Hankel(vc::AbstractVector, vr::AbstractVector) =
    Hankel{promote_type(eltype(vc), eltype(vr))}(vc, vr)

Hankel{T}(A::AbstractMatrix) where T = Hankel{T}(A[:,1], A[end,:])
Hankel(A::AbstractMatrix) = Hankel(A[:,1], A[end,:])

convert(::Type{Array}, A::Hankel) = convert(Matrix, A)
convert(::Type{Matrix}, A::Hankel{T}) where T = convert(Matrix{T}, A)
function convert(::Type{Matrix{T}}, A::Hankel) where T
    m, n = size(A)
    Af = Matrix{T}(undef, m, n)
    for j = 1:n
        for i = 1:m
            Af[i,j] = A[i,j]
        end
    end
    return Af
end

convert(::Type{AbstractArray{T}}, A::Hankel{T}) where {T<:Number} = A
convert(::Type{AbstractArray{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A)
convert(::Type{AbstractMatrix{T}}, A::Hankel{T}) where {T<:Number} = A
convert(::Type{AbstractMatrix{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A)
convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = _Hankel(convert(Toeplitz{T}, A.T))

# Size
size(H::Hankel,k...) = size(H.T,k...)

# Retrieve an entry by two indices
getindex(A::Hankel, i::Integer, j::Integer) = A.T[i,end-j+1]

# Fast application of a general Hankel matrix to a general vector
*(A::Hankel, b::AbstractVector) = A.T * reverse(b)

# Fast application of a general Hankel matrix to a general matrix
*(A::Hankel, B::AbstractMatrix) = A.T * flipdim(B, 1)
## BigFloat support

(*)(A::Toeplitz{T}, b::Vector) where {T<:BigFloat} = irfft(
    rfft([
        A.vc;
        reverse(A.vr[2:end])]
    ) .* rfft([
        b;
        zeros(length(b) - 1)
    ]),
    2 * length(b) - 1
)[1:length(b)]

function cholesky!(L::AbstractMatrix, T::SymmetricToeplitz)

    L[:, 1] .= T.vc ./ sqrt(T.vc[1])
    v = copy(L[:, 1])
    N = size(T, 1)

    @inbounds for n in 1:N-1
        sinθn = v[n + 1] / L[n, n]
        cosθn = sqrt(1 - sinθn^2)

        for n′ in n+1:N
            v[n′] = (v[n′] - sinθn * L[n′ - 1, n]) / cosθn
            L[n′, n + 1] = -sinθn * v[n′] + cosθn * L[n′ - 1, n]
        end
    end
    return Cholesky(L, 'L', 0)
end

"""
    cholesky(T::SymmetricToeplitz)

Implementation of the Bareiss Algorhithm, adapted from "On the stability of the Bareiss and
related Toeplitz factorization algorithms", Bojanczyk et al, 1993.
"""
function cholesky(T::SymmetricToeplitz)
    return cholesky!(Matrix{eltype(T)}(undef, size(T, 1), size(T, 1)), T)
end

end #module


text:
using HTSFileFormats
using Base.Test
using Bio.Seq
using BGZFStreams

function testfile(filename)
    return Pkg.dir("HTSFileFormats", "test", filename)
end

function rand_interval(intval)
    x = rand(intval)
    y = rand(intval)
    if x < y
        return x:y
    else
        return y:x
    end
end

@testset "AuxDataDict" begin
    AuxDataDict = HTSFileFormats.AuxDataDict

    dict = AuxDataDict()
    @test length(dict) == 0
    @test isempty(dict)
    @test_throws KeyError dict["NM"]

    dict = AuxDataDict(
        "X1" => UInt8(1),
        "X2" => UInt16(2),
        "X3" => UInt32(3),
        "X4" => Int8(4),
        "X5" => Int16(5),
        "X6" => Int32(6),
        "X7" => Float32(7),
        "X8" => "eight",
        "X9" => Int32[9])
    @test length(dict) == 9
    @test !isempty(dict)
    @test dict["X1"] === UInt8(1)
    @test dict["X2"] === UInt16(2)
    @test dict["X3"] === UInt32(3)
    @test dict["X4"] === Int8(4)
    @test dict["X5"] === Int16(5)
    @test dict["X6"] === Int32(6)
    @test dict["X7"] === Float32(7)
    @test dict["X8"] == "eight"
    @test typeof(dict["X8"]) == String
    @test dict["X9"] == Int32[9]
    @test typeof(dict["X9"]) == Vector{Int32}

    dict = AuxDataDict("NM" => 0x01, "XY" => Int32(100), "XZ" => [0x11, 0x23])
    @test length(dict) == 3
    @test dict["NM"] === 0x01
    @test dict["XY"] === Int32(100)
    @test dict["XZ"] == [0x11, 0x23]
    @test eltype(dict["XZ"]) == UInt8

    dict = AuxDataDict("NM" => 0x01, "MD" => "8T1T39")
    @test length(dict) == 2
    @test dict["NM"] === 0x01
    @test dict["MD"] == "8T1T39"
    dict["NM"] = 0x00
    @test dict["NM"] === 0x00
    dict["MD"] = "50"
    @test dict["MD"] == "50"
    @test collect(dict) == ["NM" => 0x00, "MD" => "50"]
    dict["XY"] = "foobar"
    @test dict["XY"] == "foobar"
    @test collect(dict) == ["NM" => 0x00, "MD" => "50", "XY" => "foobar"]
    delete!(dict, "NM")
    @test length(dict) == 2
    @test collect(dict) == ["MD" => "50", "XY" => "foobar"]
end

@testset "SAM" begin
    @testset "SAMHeader" begin
        h = SAMHeader()
        @test isa(h, Associative)
        @test isempty(h)
        h["HD"] = Dict("VN" => "100.100", "SO" => "unknown")
        @test length(h) == 1
        @test h["HD"]["VN"] == "100.100"
        h["CO"] = ["comment1", "comment2"]
        @test length(h) == 2
        @test h["CO"] == ["comment1", "comment2"]
        delete!(h, "CO")
        @test length(h) == 1
    end

    @testset "Record" begin
        rec = SAMRecord()
        @test !ismapped(rec)
        # default values
        @test seqname(rec) == ""
        @test flag(rec) == 0x0000
        @test refname(rec) == "*"
        @test nextrefname(rec) == "*"
        @test position(rec) == 0
        @test nextposition(rec) == 0
        @test templatelength(rec) == 0
        @test cigar(rec) == "*"
        @test sequence(rec) == dna""
        @test qualities(rec) == UInt8[]
    end

    @testset "Reader" begin
        reader = open(testfile("sam1.sam"), SAM)
        @test isa(reader, HTSFileFormats.SAMReader)
        n = 0
        aln = SAMRecord()
        read!(reader, aln)
        n += 1
        read!(reader, aln)
        n += 1
        @test aln["XT"] === 'U'
        @test aln["NM"] === Int32(5)
        while !eof(reader)
            read!(reader, aln)
            n += 1
        end
        close(reader)
        @test n == 200
    end

    @testset "Round trip" begin
        mktemp() do path, io
            # copy
            reader = open(testfile("sam1.sam"), SAM)
            writer = HTSFileFormats.SAMWriter(io)
            write(writer, header(reader))
            aln = SAMRecord()
            while !eof(reader)
                read!(reader, aln)
                write(writer, aln)
            end
            close(reader)
            close(writer)

            # read again
            reader = open(path, SAM)
            n = 0
            aln = SAMRecord()
            while !eof(reader)
                read!(reader, aln)
                n += 1
            end
            close(reader)
            @test n == 200
        end
    end
end

@testset "BAM" begin
    @testset "Record" begin
        rec = BAMRecord()
        @test !ismapped(rec)
        # default values
        @test refname(rec) == "*"
        @test refid(rec) == 0
        @test position(rec) == 0
        @test mappingquality(rec) == 0
        @test flag(rec) == 0
        @test nextrefname(rec) == "*"
        @test nextrefid(rec) == 0
        @test nextposition(rec) == 0
        @test templatelength(rec) == 0
        @test seqname(rec) == ""
        @test cigar(rec) == ""
        @test sequence(rec) == dna""
        @test qualities(rec) == UInt8[]

        @test HTSFileFormats.rightmost_position(rec) === Int32(-1)
        @test HTSFileFormats.alignment_length(rec) === 0

        rec = BAMRecord()
        @test !haskey(rec, "MN")
        rec["MN"] = 0x01
        @test rec["MN"] === 0x01
        @test haskey(rec, "MN")
        delete!(rec, "MN")
        @test !haskey(rec, "MN")
    end

    @testset "Reader" begin
        reader = open(testfile("bam1.bam"), BAM)
        h = header(reader)
        @test h["SQ"] == [Dict("SN" => "1", "LN" => "239940")]
        @test h["PG"] == [Dict("ID" => "bwa", "PN" => "bwa", "VN" => "0.6.2-r126")]
        @test isa(reader, HTSFileFormats.BAMReader)
        n = 0
        aln = BAMRecord()
        read!(reader, aln)
        n += 1
        read!(reader, aln)
        n += 1
        @test aln["XT"] === 'U'
        # TODO: check type
        @test aln["NM"] == 5
        while !eof(reader)
            read!(reader, aln)
            n += 1
        end
        close(reader)
        @test n == 200
    end

    @testset "Round trip" begin
        mktemp() do path, _
            # copy
            reader = open(testfile("bam1.bam"), BAM)
            writer = HTSFileFormats.BAMWriter(
                BGZFStream(path, "w"),
                header(reader, true))
            aln = BAMRecord()
            while !eof(reader)
                read!(reader, aln)
                write(writer, aln)
            end
            close(reader)
            close(writer)

            # read again
            reader = open(path, BAM)
            n = 0
            aln = BAMRecord()
            while !eof(reader)
                read!(reader, aln)
                n += 1
            end
            close(reader)
            @test n == 200
        end
    end

    @testset "random access" begin
        reader = open(testfile("test.bam"), BAM)
        for (name, interval, expected) in [
                ("chr1", 1000:10000, 21),
                ("chr1", 8000:10000, 20),
                ("chr1", 766_000:800_000, 142),
                ("chr1", 786_000:800_000, 1),
                ("chr1", 796_000:800_000, 0)]
            n = 0
            for rec in intersect(reader, name, interval)
                n += 1
            end
            @test n == expected
        end

        for n in 1:50
            name = "chr1"
            refid = 1
            interval = rand_interval(1:1_000_000)

            expected = BAMRecord[]
            seekstart(reader)
            for rec in reader
                if HTSFileFormats.isoverlapping(rec, refid, interval)
                    push!(expected, rec)
                end
            end

            actual = BAMRecord[]
            for rec in intersect(reader, name, interval)
                push!(actual, rec)
            end

            @test map(seqname, actual) == map(seqname, expected)
        end
    end
end

@testset "Tabix" begin
    bedfile = testfile("knownGene.part.bed.gz")
    indexfile = bedfile * ".tbi"

    @testset "load" begin
        index = Tabix(indexfile)
        @test isa(index, Tabix)
        @test index.format === Int32(0x10000)
        @test index.columns === (1, 2, 3)
        @test index.meta === '#'
        @test index.skip === 0
        @test index.names == ["chr1", "chr2", "chr3"]
        @test length(index.indexes) == 3
    end

    @testset "chunks" begin
        stream = BGZFStream(bedfile)
        index = Tabix(indexfile)

        # Expected values were counted using the tabix tool as follows:
        #   $ tabix test/knownGene.part.bed.gz chr1:5,000,000-10,000,000 | wc -l
        #   $     64
        for (seqname, interval, expected) in [
                ("chr1", 5_000_000:10_000_000, 64),
                ("chr1", 6_000_000:10_000_000, 54),
                ("chr2", 5_000_000:10_000_000, 81),
                ("chr3", 9_000_000:10_000_000, 15),
                ("chr3", 9_800_000:10_000_000,  3),
                ("chr3", 9_900_000:10_000_000,  0)]
            n = 0
            for chunk in overlapchunks(index, seqname, interval)
                seek(stream, chunk)
                while virtualoffset(stream) in chunk
                    line = readline(stream)
                    values = split(chomp(line), '\t')
                    nm = values[index.columns[1]]
                    @test nm == seqname
                    int = parse(Int, values[index.columns[2]]):parse(Int, values[index.columns[3]])-1
                    if !isempty(intersect(int, interval))
                        n += 1
                    elseif first(int) > last(interval)
                        break
                    end
                end
            end
            @test n == expected
        end
    end
end


text:
function parseinput(lines)
    idx = findall(==(""), lines)[1]
    rules = lines[1:idx-1]
    messages = lines[idx+1:end]
    return (rules, messages)
end

function recpattern(v, mapping)
    vals = split(v, " ")
    pattern = join(mapping[a] for a in vals)
    return pattern
end

function makemapping(rules)
    starting_rules = filter(x -> '"' ∈ x, rules)
    recursive_rules= filter(x -> !occursin("\"", x), rules)
    mapping = Dict{String, String}()
    for r in starting_rules
        k, v = split(r, ": \"")
        mapping[k] = v[1:1]
    end
    while length(keys(mapping)) != length(rules)
        nextrules(recursive_rules, mapping)
    end
    return mapping
end

function nextrules(rules, mapping)
    possiblevalues = [keys(mapping)...,"|"]
    for r in rules
        k, v = split(r, ": ")
        if k in keys(mapping)
	    continue
	end
       if all(split(v, " ") .∈ Ref(possiblevalues))
           if occursin("|", v)
               v1, v2  = split(v, " | ")
	       p1 = recpattern(v1, mapping)
               p2 = recpattern(v2, mapping)
	       pattern = "($p1|$p2)"
           else
               pattern = recpattern(v, mapping)
	   end
           mapping[k] = pattern
        end
    end
    return mapping
end

function part1(lines)
    rules, messages = parseinput(lines)
    mapping = makemapping(rules)
    regex = Regex("^$(mapping["0"])\$")
    total = sum(occursin(regex, m) for m in messages)
    return total
end

data = readlines("day19.txt")

println("Part 1: $(part1(data))")

# Horrible solution but works
# 0 is 8 11
# the new 8 is 8+ (regex 8 once or more)
# the new 11 is 42 {i} 31 {i}
# making the regex once is not possible as it
# will be too big so I make a lot of regexes

function makeregex(mapping, i)
    d = "($(mapping["8"]))+"
    a = mapping["42"]
    b = mapping["31"]
    c = "($a){$i}($b){$i}"
    m = "^$d$c\$"
    return Regex(m)
end

function part2(lines)
    rules, messages = parseinput(lines)
    mapping = makemapping(rules)
    regexes =[makeregex(mapping, i) for i in 1:10]
    total = 0
    for m in messages
        if any(occursin.(regexes, m))
	    total += 1
	end
    end
    return total
end

println("Part2 : $(part2(data))")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Deployment <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, DeploymentProperties } # spec name: properties

    function Deployment(;properties=nothing)
        o = new()
        validate_property(Deployment, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type Deployment

const _property_map_Deployment = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"))
const _property_types_Deployment = Dict{Symbol,String}(Symbol("properties")=>"DeploymentProperties")
Base.propertynames(::Type{ Deployment }) = collect(keys(_property_map_Deployment))
Swagger.property_type(::Type{ Deployment }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_Deployment[name]))}
Swagger.field_name(::Type{ Deployment }, property_name::Symbol) =  _property_map_Deployment[property_name]

function check_required(o::Deployment)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

function validate_property(::Type{ Deployment }, name::Symbol, val)
end


text:
function HMMforward(v,phghm,ph1,pvgh;UseLogArray=true)
    #function [logalpha,loglik]=HMMforward(v,phghm,ph1,pvgh)
    #%HMMFORWARD HMM Forward Pass
    #% [alpha,loglik]=HMMforward(v,phghm,ph1,pvgh)
    #%
    #% Inputs:
    #% v : visible (observation) sequence being a vector v=[2 1 3 3 1 ...]
    #% phghm : homogeneous transition distribution phghm(i,j)=p(h(t)=i|h(t-1)=j)
    #% ph1 : initial distribution
    #% pvgh : homogeneous emission distribution pvgh(i,j)=p(v(t)=i|h(t)=j)
    #%
    #% Outputs:
    #% logalpha : log alpha messages: log p(h(t),v(1:t))
    #% loglik : sequence log likelihood log p(v(1:T))
    #% See also HMMbackward.m, HMMviterbi.m,  HMMsmooth.m, demoHMMinference.m

    T=length(v); H=length(ph1);

    if !UseLogArray # alpha recursion (not recommended due to numerical underflow)
        z=zeros(1,T)
        alpha=zeros(H,T)
        alpha[:,1] = pvgh[v[1],:]'.*ph1
        z[1]=sum(alpha[:,1])
        alpha[:,1]=condp(alpha[:,1])
        for t=2:T
            alpha[:,t]=pvgh[v[t],:]'.*(phghm*alpha[:,t-1])
            z[t]=sum(alpha[:,t])
            alpha[:,t]=condp(alpha[:,t])
        end
        loglik = sum(log(z)) # log likelihood
    end

    if UseLogArray
        logpvgh=LogArray(log(pvgh))
        logphghm=LogArray(log(phghm))
        logph1=LogArray(log(ph1))
        logalpha=Dict{Integer,LogArray}()
        logalphaoutput=zeros(H,T)
        logalpha[1] = logpvgh[v[1],:].'.*logph1
        logalphaoutput[:,1]=(logalpha[1]).content
        for t=2:T
            logalpha[t]=logpvgh[v[t],:].'.*(logphghm*logalpha[t-1])
            #if any(isnan(logalpha[t].content))
            #    return t,v,logpvgh,logphghm,logalpha
            #end

            logalphaoutput[:,t]=(logalpha[t]).content
        end
        loglik = (sum(logalpha[T])).content # log likelihood
    end

    return logalphaoutput, loglik

end


text:
import GreyDecision.Utility: makeminmax
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Marcos: marcos

@testset "Marcos" begin
    tol = 0.0001

    weights = [0.19019, 0.15915, 0.19819, 0.19019, 0.15115, 0.11111]

    fns = [maximum, maximum, maximum, maximum, maximum, maximum]
    Fns = convert(Array{Function,1}, fns)

    decmat = [
        8.675 8.433 8.000 7.800 8.025 8.043
        8.825 8.600 7.420 7.463 7.825 8.229
        8.325 7.600 8.040 7.700 7.925 7.600
        8.525 8.667 7.180 7.375 7.750 8.071
    ]

    nrow, ncol = size(decmat)
    # Converting white decision matrix to
    # grey decision matrix
    greydecmat = Array{GreyNumber{Float64}, 2}(undef, nrow, ncol)

    for i in 1:nrow
        for j in 1:ncol
            greydecmat[i, j] = GreyNumber(decmat[i, j])
        end
    end

    result = marcos(greydecmat, weights, Fns)

    @test result[:bestIndex] == 1
    @test isapprox(result[:scores][1].a, 0.684865943528, atol = tol)
    @test isapprox(result[:scores][2].a, 0.672767106696, atol = tol)
    @test isapprox(result[:scores][3].a, 0.662596906139, atol = tol)
    @test isapprox(result[:scores][4].a, 0.661103207660, atol = tol)

end


text:
struct LittleInt <: Integer
    val::Int8
    function LittleInt(n::Real)
        1 ≤ n ≤ 10 || throw(ArgumentError("LittleInt number must be in [1, 10]"))
        return new(Int8(n))
    end
end
Base.show(io::IO, x::LittleInt) = print(io, x.val)
Base.convert(::Type{T}, x::LittleInt) where T<:Number = convert(T, x.val)
Base.promote_rule(::Type{LittleInt}, ::Type{T}) where T<:Number = T

for op in (:+, :*, :÷, :-, :&, :|, :$, :<, :>, :(==))
    @eval (Base.$op)(a::LittleInt, b::LittleInt) = LittleInt(($op)(a.val, b.val))
end

# Test
a = LittleInt(3)
b = LittleInt(4.0)
@show a b
@show a + b
@show b - a
@show a * LittleInt(2)
@show b ÷ LittleInt(2)
@show a * b

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Random
using Test

using DataStructuresAlgorithms

include("search.jl")
include("sort.jl")
include("stack.jl")
include("dynamicarray.jl")
include("list.jl")
include("queue.jl")
include("tree.jl")
include("heap.jl")
include("hash.jl")
include("graph.jl")


text:
using OceanTurb, Printf

using OceanTurb.ModularKPP: DiagnosticPlumeModel

@use_pyplot_utils

usecmbright()

modelsetup = (N=16, L=128, stepper=:BackwardEuler)

Fb = 1e-7
Fu = 0.0
N² = 1e-5
Δt = 2minute
t_plot = (1hour, 2hour, 3hour)

name = "Free convection
    \\small{with \$ \\overline{w b} |_{z=0} = 10^{-7} \\, \\mathrm{m^2 \\, s^{-3}}\$}"

# The standard setup except with a plume model rather than a counter-gradient flux model.
model = ModularKPP.Model(; N=1024, L=64, stepper=:BackwardEuler,
                         nonlocalflux = DiagnosticPlumeModel())

# Initial condition and fluxes
dTdz = N² / (model.constants.α * model.constants.g)
T₀(z) = 20 + dTdz * z

model.solution.T = T₀

Fθ = Fb / (model.constants.α * model.constants.g)
model.bcs.U.top = FluxBoundaryCondition(Fu)

model.bcs.T.top = FluxBoundaryCondition(Fθ)
model.bcs.T.bottom = GradientBoundaryCondition(dTdz)

fig, axs = subplots(ncols=3)

sca(axs[1])
removespines("top", "right")
xlabel("Temperature \$ \\, {}^\\circ \\mathrm{C} \$")
ylabel(L"z \, \mathrm{(m)}")

sca(axs[3])
removespines("top", "right", "left")
axs[3].tick_params(left=false, labelleft=false)
xlabel("Vertical momentum")

sca(axs[2])
removespines("top", "right", "left")
axs[2].tick_params(left=false, labelleft=false)
xlabel("Plume temperature excess")

function iterate_and_plot(Δt, n=1)
    iterate!(model, Δt=Δt, Nt=1)

    sca(axs[1])
    plot(model.solution.T, "-", alpha=0.8, markersize=1.5)
    plot(model.state.plume.T, "--", alpha=0.8, markersize=1.5)

    ΔT = model.state.plume.T - model.solution.T
    sca(axs[2])
    plot(ΔT, "-", alpha=0.8, markersize=1.5)

    sca(axs[3])
    plot(model.state.plume.W², "-")

    return nothing
end


text:
using MacroTools
using Swizzles

macro sugar(ex)
    op = nothing
    reg = nothing
    vars = nothing
    init = nothing
    elem = nothing

        if @capture(ex, [vars__] = elem_)
    elseif @capture(ex, [vars__] += elem_) op = :+
    elseif @capture(ex, [vars__] *= elem_) op = :*
    elseif @capture(ex, [vars__] |= elem_) op = :|
    elseif @capture(ex, [vars__] &= elem_) op = :&
    elseif @capture(ex, [vars__] <op_>= elem_)
    elseif @capture(ex, [vars__] = init_ += elem_) op = :+
    elseif @capture(ex, [vars__] = init_ *= elem_) op = :*
    elseif @capture(ex, [vars__] = init_ |= elem_) op = :|
    elseif @capture(ex, [vars__] = init_ &= elem_) op = :&
    elseif @capture(ex, [vars__] = init_ <op_>= elem_)
    elseif @capture(ex, reg_[vars__] = elem_) init = :($reg[$(vars...)])
    elseif @capture(ex, reg_[vars__] += elem_) init = :($reg[$(vars...)]); op = :+
    elseif @capture(ex, reg_[vars__] *= elem_) init = :($reg[$(vars...)]); op = :*
    elseif @capture(ex, reg_[vars__] |= elem_) init = :($reg[$(vars...)]); op = :|
    elseif @capture(ex, reg_[vars__] &= elem_) init = :($reg[$(vars...)]); op = :&
    elseif @capture(ex, reg_[vars__] <op_>= elem_) init = :($reg[$(vars...)])
    elseif @capture(ex, reg_[vars__] = init_ += elem_) op = :+
    elseif @capture(ex, reg_[vars__] = init_ *= elem_) op = :*
    elseif @capture(ex, reg_[vars__] = init_ |= elem_) op = :|
    elseif @capture(ex, reg_[vars__] = init_ &= elem_) op = :&
    elseif @capture(ex, reg_[vars__] = init_ <op_>= elem_)
    else
        return :(error("TODO"))
    end

    function dim(var, env)
        if var == :(:)
            return env[gensym()] = length(env) + 1
        elseif var isa Symbol
            if !haskey(dims, var)
                env[var] = length(env) + 1
            end
            return universe[var]
        elseif var isa Integer
            return nil
        end
        error("TODO")
    end

    N = length(vars)
    Dict(vars[n] => n for n in 1:N)

    function beamify(ex, env)
        if @capture(ex, arr_[inds__])
            return :($(Beam(map(ind -> dim(ind, env), inds)...)).($arr))
        elseif ex isa Expr
	    return Expr(ex.head, map(arg->beamify(arg, env), ex.args)...)
        else
            return ex
        end
    end

    init_env = Dict(vars[n] => n for n in 1:N)
    if init === nothing
        init = ()
    else
        init = (Base.Broadcast.__dot__(beamify(init, init_env)), )
    end

    elem_env = Dict()
    mask = Val((map(var-> dim(var, elem_env), vars)...,))
    elem = Base.Broadcast.__dot__(beamify(elem, elem_env))

    if reg === nothing
        return esc(:(Swizzle($op, $mask).($(init...), $elem)))
    else
        return esc(:($reg .= Swizzle($op, $mask).($(init...), $elem)))
    end
end


text:
# compare performance of graph traversal

using Graphs

# construct graphs

nv = 10000
deg = 100

g_adj = simple_adjlist(nv)
g_inc = simple_inclist(nv)
g_gra = simple_graph(nv)

shuff = [1:nv]

for u = 1 : nv
    for j in 1 :  deg

        v = u

        while v == u
            k = rand(j+1:nv)
            shuff[j], shuff[k] = shuff[k], shuff[j]
            v = shuff[j]
        end
        add_edge!(g_adj, u, v)
        add_edge!(g_inc, u, v)
        add_edge!(g_gra, u, v)
    end
end

eweights = rand(num_edges(g_inc))

# test tasks

function neighbor_scan(g::AbstractGraph)
    for v in vertices(g)
        for u in out_neighbors(v, g)
            u
        end
    end
end

function outedge_scan(g::AbstractGraph)
    for v in vertices(g)
        for u in out_edges(v, g)
            e
        end
    end
end

function bf_traverse(g::AbstractGraph)
    traverse_graph(g, BreadthFirst(), 1, TrivialGraphVisitor())
end

function df_traverse(g::AbstractGraph)
    traverse_graph(g, DepthFirst(), 1, TrivialGraphVisitor())
end

function run_dijkstra(g::AbstractGraph)
    dijkstra_shortest_paths(g, eweights, 1)
end

# generic benchmark macro

macro graph_bench(fun, name, g, repeats)
    quote
        $(fun)($g)  # warming
        et = @elapsed for i = 1 : $repeats
            $(fun)($g)
        end
        @printf("    On %-25s: avg elapsed = %.4f ms\n", $(string(name)), 1000 * et / $repeats)
    end
end

# benchmarks

println("Benchmark on a graph with $(nv) vertices and $(nv * deg) edges")

println("Benchmark of neighbor scan")
@graph_bench neighbor_scan SimpleAdjacencyList g_adj 50
@graph_bench neighbor_scan SimpleIncidenceList g_inc 50
@graph_bench neighbor_scan SimpleGraph g_gra 50

println("Benchmark of out-edge scan")
@graph_bench outedge_scan SimpleIncidenceList g_inc 50
@graph_bench outedge_scan SimpleGraph g_gra 50

println("Benchmark of breadth-first traversal")
@graph_bench bf_traverse SimpleAdjacencyList g_adj 10
@graph_bench bf_traverse SimpleIncidenceList g_inc 10
@graph_bench bf_traverse SimpleGraph g_gra 10

println("Benchmark of depth-first traversal")
@graph_bench df_traverse SimpleAdjacencyList g_adj 5
@graph_bench df_traverse SimpleIncidenceList g_inc 2
@graph_bench df_traverse SimpleGraph g_gra 5

println("Benchmark of Dijkstra shortest paths")
@graph_bench run_dijkstra SimpleIncidenceList g_inc 5
@graph_bench run_dijkstra SimpleGraph g_gra 5


text:
# Autogenerated wrapper script for Pikchr_jll for armv7l-linux-musleabihf
export libpikchr

JLLWrappers.@generate_wrapper_header("Pikchr")
JLLWrappers.@declare_library_product(libpikchr, "libpikchr.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpikchr,
        "lib/libpikchr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const require = ["Atom", "Gadfly"]

cachedir = Pkg.dir("..", "lib", "v$(VERSION.major).$(VERSION.minor)")

isdir(cachedir) && rm(cachedir, recursive = true)

Pkg.init()
Pkg.update()
for pkg in require
  Pkg.add(pkg)
  Base.require(symbol(pkg))
end

# Cache Preprocessing

function startswith!(io::IO, s)
  pos = position(io)
  for c in s
    if eof(io) || c ≠ read(io, UInt8)
      seek(io, pos)
      return false
    end
  end
  return true
end

startswith!(io::IO, s::UTF8String) = startswith!(io, s.data)

prefix = Pkg.dir()
dummy = "i sure hope this isn't in any packages"

function upto(io::IO)
  len = position(io)
  seek(io, 0)
  return readbytes(io, len)
end

cd(cachedir) do
  for cache in readdir()
    f = open(cache)
    Base.isvalid_cache_header(f) || continue
    Base.cache_dependencies(f)

    out = IOBuffer()
    write(out, upto(f))
    io = IOBuffer(readbytes(f))
    close(f)

    while !eof(io)
      if startswith!(io, prefix)
        write(out, rpad(dummy, length(prefix)))
      else
        write(out, read(io, UInt8))
      end
    end
    open(io -> write(io, takebuf_array(out)), cache, "w")
  end
end

exit()

max_stars_repo_name: JunoLab/atom

max_issues_repo_name: JunoLab/atom

max_forks_repo_name: JunoLab/atom


text:
abstract type AbstractPenalization end
abstract type AbstractPenalizer end
# Penalization API #
# ------------------------------------------------------------
penalty(ph::AbstractProgressiveHedgingSolver) = penalty(ph, ph.penalization)
init_penalty!(ph::AbstractProgressiveHedgingSolver) = init_penalty!(ph, ph.penalization)
update_penalty!(ph::AbstractProgressiveHedgingSolver) = update_penalty!(ph, ph.penalization)
# ------------------------------------------------------------
include("fixed.jl")
include("adaptive.jl")


text:
using Test
using LinearAlgebra, StaticArrays
import ClimaCore.Spaces: Quadratures

f(x) = x^3 + 2x^2 + 3x + 4
fd(x) = 3x^2 + 4x + 3

@testset "test LGL quadrature" begin
    quad = Quadratures.GLL{3}()
    @test quad isa Quadratures.QuadratureStyle
    @test Quadratures.polynomial_degree(quad) == 2
    @test Quadratures.degrees_of_freedom(quad) == 3

    points, weights = Quadratures.quadrature_points(Float64, quad)

    @test eltype(points) === eltype(weights) === Float64
    @test length(points) == length(weights) == 3
    @test dot(f.(points), weights) ≈ 28 / 3
end

@testset "test GL quadrature" begin
    quad = Quadratures.GL{4}()
    @test quad isa Quadratures.QuadratureStyle
    @test Quadratures.polynomial_degree(quad) == 3
    @test Quadratures.degrees_of_freedom(quad) == 4

    points, weights = Quadratures.quadrature_points(Float32, quad)

    @test eltype(points) === eltype(weights) === Float32
    @test length(points) == length(weights) == 4
    @test dot(f.(points), weights) ≈ 28 / 3
end

@testset "differentiation_matrix" begin
    quad = Quadratures.GL{4}()
    points, weights = Quadratures.quadrature_points(Float64, quad)

    D = Quadratures.differentiation_matrix(Float64, quad)
    @test D isa SMatrix
    @test size(D) == (4, 4)

    @test D * f.(points) ≈ fd.(points)

    quad = Quadratures.GLL{4}()
    points, weights = Quadratures.quadrature_points(Float64, quad)

    D = Quadratures.differentiation_matrix(Float64, quad)
    @test D isa SMatrix
    @test size(D) == (4, 4)

    @test D * f.(points) ≈ fd.(points)
end

@testset "interpolation matrix" begin
    quad1 = Quadratures.GL{4}()
    quad2 = Quadratures.GLL{5}()
    points1, weights1 = Quadratures.quadrature_points(Float64, quad1)
    points2, weights2 = Quadratures.quadrature_points(Float64, quad2)
    I = Quadratures.interpolation_matrix(Float64, quad2, quad1)
    @test size(I) == (5, 4)

    @test I * f.(points1) ≈ f.(points2)

end


text:
using Test
using Metatheory
using Metatheory.Library
using Metatheory.Schedulers
using TermInterface

mult_t = @commutative_monoid (*) 1
plus_t = @commutative_monoid (+) 0

minus_t = @theory a b begin
  # TODO Jacques Carette's post in zulip chat
  a - a --> 0
  a - b --> a + (-1 * b)
  -a --> -1 * a
  a + (-b) --> a + (-1 * b)
end

mulplus_t = @theory a b c begin
  # TODO FIXME this rules improves performance and avoids commutative
  # explosion of the egraph
  a + a --> 2 * a
  0 * a --> 0
  a * 0 --> 0
  a * (b + c) == ((a * b) + (a * c))
  a + (b * a) --> ((b + 1) * a)
end

pow_t = @theory x y z n m p q begin
  (y^n) * y --> y^(n + 1)
  x^n * x^m == x^(n + m)
  (x * y)^z == x^z * y^z
  (x^p)^q == x^(p * q)
  x^0 --> 1
  0^x --> 0
  1^x --> 1
  x^1 --> x
  x * x --> x^2
  inv(x) == x^(-1)
end

div_t = @theory x y z begin
  x / 1 --> x
  # x / x => 1 TODO SIGN ANALYSIS
  x / (x / y) --> y
  x * (y / x) --> y
  x * (y / z) == (x * y) / z
  x^(-1) == 1 / x
end

trig_t = @theory θ begin
  sin(θ)^2 + cos(θ)^2 --> 1
  sin(θ)^2 - 1 --> cos(θ)^2
  cos(θ)^2 - 1 --> sin(θ)^2
  tan(θ)^2 - sec(θ)^2 --> 1
  tan(θ)^2 + 1 --> sec(θ)^2
  sec(θ)^2 - 1 --> tan(θ)^2
  cot(θ)^2 - csc(θ)^2 --> 1
  cot(θ)^2 + 1 --> csc(θ)^2
  csc(θ)^2 - 1 --> cot(θ)^2
end

# Dynamic rules
fold_t = @theory a b begin
  -(a::Number) => -a
  a::Number + b::Number => a + b
  a::Number * b::Number => a * b
  a::Number^b::Number => begin
    b < 0 && a isa Int && (a = float(a))
    a^b
  end
  a::Number / b::Number => a / b
end

using Calculus: differentiate
diff_t = @theory x y begin
  ∂(y, x::Symbol) => begin
    z = extract!(_egraph, simplcost; root = y.id)
    @show z
    zd = differentiate(z, x)
    @show zd
    zd
  end
end

cas = fold_t ∪ mult_t ∪ plus_t ∪ minus_t ∪ mulplus_t ∪ pow_t ∪ div_t ∪ trig_t ∪ diff_t

function customlt(x, y)
  if typeof(x) == Expr && typeof(y) == Expr
    false
  elseif typeof(x) == typeof(y)
    isless(x, y)
  elseif x isa Symbol && y isa Number
    false
  elseif x isa Expr && y isa Number
    false
  elseif x isa Expr && y isa Symbol
    false
  else
    true
  end
end

canonical_t = @theory x y n xs ys begin
  # restore n-arity
  (x * x) --> x^2
  (x^n::Number * x) --> x^(n + 1)
  (x * x^n::Number) --> x^(n + 1)
  (x + (+)(ys...)) --> +(x, ys...)
  ((+)(xs...) + y) --> +(xs..., y)
  (x * (*)(ys...)) --> *(x, ys...)
  ((*)(xs...) * y) --> *(xs..., y)

  (*)(xs...) => Expr(:call, :*, sort!(xs; lt = customlt)...)
  (+)(xs...) => Expr(:call, :+, sort!(xs; lt = customlt)...)
end

function simplcost(n::ENodeTerm, g::EGraph)
  cost = 0 + arity(n)
  if operation(n) == :∂
    cost += 20
  end
  for id in arguments(n)
    eclass = g[id]
    !hasdata(eclass, simplcost) && (cost += Inf; break)
    cost += last(getdata(eclass, simplcost))
  end
  return cost
end

simplcost(n::ENodeLiteral, g::EGraph) = 0

function simplify(ex; steps = 4)
  params = SaturationParams(
    scheduler = ScoredScheduler,
    eclasslimit = 5000,
    timeout = 7,
    schedulerparams = (1000, 5, Schedulers.exprsize),
    #stopwhen=stopwhen,
  )
  hist = UInt64[]
  push!(hist, hash(ex))
  for i in 1:steps
    g = EGraph(ex)
    @profview saturate!(g, cas, params)
    ex = extract!(g, simplcost)
    ex = rewrite(ex, canonical_t)
    if !TermInterface.istree(ex)
      return ex
    end
    if hash(ex) ∈ hist
      println("loop detected $ex")
      return ex
    end
    println(ex)
    push!(hist, hash(ex))
  end

end

@test :(4a) == simplify(:(2a + a + a))
@test :(a * b * c) == simplify(:(a * c * b))
@test :(2x) == simplify(:(1 * x * 2))
@test :((a * b)^2) == simplify(:((a * b)^2))
@test :((a * b)^6) == simplify(:((a^2 * b^2)^3))
@test :(a + b + d) == simplify(:(a + b + (0 * c) + d))
@test :(a + b) == simplify(:(a + b + (c * 0) + d - d))
@test :(a) == simplify(:((a + d) - d))
@test :(a + b + d) == simplify(:(a + b * c^0 + d))
@test :(a * b * x^(d + y)) == simplify(:(a * x^y * b * x^d))
@test :(a * b * x^74103) == simplify(:(a * x^(12 + 3) * b * x^(42^3)))

@test 1 == simplify(:((x + y)^(a * 0) / (y + x)^0))
@test 2 == simplify(:(cos(x)^2 + 1 + sin(x)^2))
@test 2 == simplify(:(cos(y)^2 + 1 + sin(y)^2))
@test 2 == simplify(:(sin(y)^2 + cos(y)^2 + 1))

@test :(y + sec(x)^2) == simplify(:(1 + y + tan(x)^2))
@test :(y + csc(x)^2) == simplify(:(1 + y + cot(x)^2))

# simplify(:( ∂(x^2, x)))

@time simplify(:(∂(x^(cos(x)), x)))

@test :(2x^3) == simplify(:(x * ∂(x^2, x) * x))

# @simplify ∂(y^3, y) * ∂(x^2 + 2, x) / y * x

# @simplify (6 * x * x * y)

# @simplify ∂(y^3, y) / y

# # ex = :( ∂(x^(cos(x)), x) )
# ex = :( (6 * x * x * y) )
# g = EGraph(ex)
# saturate!(g, cas)
# g.classes
# extract!(g, simplcost; root=g.root)

# params = SaturationParams(
#     scheduler=BackoffScheduler,
#     eclasslimit=5000,
#     timeout=7,
#     schedulerparams=(1000,5),
#     #stopwhen=stopwhen,
# )

# ex = :((x+y)^(a*0) / (y+x)^0)
# g = EGraph(ex)
# @profview println(saturate!(g, cas, params))

# ex = extract!(g, simplcost)
# ex = rewrite(ex, canonical_t; clean=false)

# FIXME this is a hack to get the test to work.
if VERSION < v"1.9.0-DEV"
  function EGraphs.make(::Val{:type_analysis}, g::EGraph, n::ENodeTerm)
    Any
  end

  function EGraphs.make(::Val{:type_analysis}, g::EGraph, n::ENodeLiteral)
    v = n.value
    if v == :im
      typeof(im)
    else
      typeof(v)
    end
  end

  function EGraphs.make(::Val{:type_analysis}, g::EGraph, n::ENodeTerm{Expr})
    if exprhead(n) != :call
      # println("$n is not a call")
      t = Any
      # println("analyzed type of $n is $t")
      return t
    end
    sym = operation(n)
    if !(sym isa Symbol)
      # println("head $sym is not a symbol")
      t = Any
      # println("analyzed type of $n is $t")
      return t
    end

    symval = getfield(@__MODULE__, sym)
    child_classes = map(x -> g[x], arguments(n))
    child_types = Tuple(map(x -> getdata(x, :type_analysis, Any), child_classes))

    # t = t_arr[1]
    t = Core.Compiler.return_type(symval, child_types)

    if t == Union{}
      throw(MethodError(symval, child_types))
    end
    # println("analyzed type of $n is $t")
    return t
  end

  EGraphs.join(::Val{:type_analysis}, from, to) = typejoin(from, to)

  EGraphs.islazy(::Val{:type_analysis}) = true

  function infer(e)
    g = EGraph(e)
    analyze!(g, :type_analysis)
    getdata(g[g.root], :type_analysis)
  end

  ex1 = :(cos(1 + 3.0) + 4 + (4 - 4im))
  ex2 = :("ciao" * 2)
  ex3 = :("ciao" * " mondo")

  @test ComplexF64 == infer(ex1)
  @test_throws MethodError infer(ex2)
  @test String == infer(ex3)
end