text:
function get_rve_dim_tag(bbox::NTuple{6,Float64})
    x_min, y_min, z_min, x_max, y_max, z_max = bbox
    if (x_min != x_max) && (y_min != y_max) && (z_min != z_max)
        return "XYZ"
    elseif (x_min == x_max && y_min != y_max && z_min != z_max)
        return "YZ"
    elseif (x_min != x_max && y_min == y_max && z_min != z_max)
        return "ZX"
    elseif (x_min != x_max && y_min != y_max && z_min == z_max)
        return "XY"
    end
end

"""
Check periodicity and get the node pairs to be contrained
Tasks:
    1. Get the groups of nodes on faces, edges and vertices
    2. Verify, if the number of nodes on opposite faces, edges match
    3. Verify, if the vertices groups has single node
    4. Verify, if a node on one face is pairing with morethan one node on opposite faces
    5.

Returns:
--------

3D Node groups as a dictionary with the following keys\n
`outer_nodes`,\n
faces: `xp`, `yp`, `zp`, `xn`, `yn`, `zn` \n
edges: `zpxp`, `zpxn`, `znxp`, `znxn`, `xpyp`, `xpyn`, `xnyp`, `xnyn`, `ypzp`, `ypzn`, `ynzp`, `ynzn`,\n
vertices: `xpypzp`, `xpypzn`, `xpynzp`, `xnypzp`, `xnynzp`, `xnypzn`, `xnypzp`, `xnynzn`.

2D Node groups for zx-plane as a dictionary with the following keys \n
`outer_nodes`,\n
edges: `xp`, `xn`, `yp`, `yn`, \n
vertices: `xpyp`, `xpyn`, `xnyp`, `xnyn`.

"""
function get_rve_node_groups(
    node_tags::Vector{Int},
    node_coordinates::Matrix{Float64},
    rve_bbox::NTuple{6,Float64};
    small_par::Float64=1e-06
)::Dict{String,Matrix{IF64}}
    x_min, y_min, z_min, x_max, y_max, z_max = rve_bbox
    xlb_node::Bool, ylb_node::Bool, zlb_node::Bool = false, false, false
    xub_node::Bool, yub_node::Bool, zub_node::Bool = false, false, false
    #
    a_node_coordinates::NTuple{3,Float64} = (0.0, 0.0, 0.0)
    xn_ntags::Vector{Int} = Int[]
    xn_ncoor::Vector{Float64} = Float64[]
    yn_ntags::Vector{Int} = Int[]
    yn_ncoor::Vector{Float64} = Float64[]
    zn_ntags::Vector{Int} = Int[]
    zn_ncoor::Vector{Float64} = Float64[]
    xp_ntags::Vector{Int} = Int[]
    xp_ncoor::Vector{Float64} = Float64[]
    yp_ntags::Vector{Int} = Int[]
    yp_ncoor::Vector{Float64} = Float64[]
    zp_ntags::Vector{Int} = Int[]
    zp_ncoor::Vector{Float64} = Float64[]
    #
    # edges about x-axis
    ynzn_ntags::Vector{Int} = Int[]
    ynzn_ncoor::Vector{Float64} = Float64[]
    ypzn_ntags::Vector{Int} = Int[]
    ypzn_ncoor::Vector{Float64} = Float64[]
    ypzp_ntags::Vector{Int} = Int[]
    ypzp_ncoor::Vector{Float64} = Float64[]
    ynzp_ntags::Vector{Int} = Int[]
    ynzp_ncoor::Vector{Float64} = Float64[]
    # edges about y-axis
    znxn_ntags::Vector{Int} = Int[]
    znxn_ncoor::Vector{Float64} = Float64[]
    zpxn_ntags::Vector{Int} = Int[]
    zpxn_ncoor::Vector{Float64} = Float64[]
    zpxp_ntags::Vector{Int} = Int[]
    zpxp_ncoor::Vector{Float64} = Float64[]
    znxp_ntags::Vector{Int} = Int[]
    znxp_ncoor::Vector{Float64} = Float64[]
    # edges about z-axis
    xnyn_ntags::Vector{Int} = Int[]
    xnyn_ncoor::Vector{Float64} = Float64[]
    xpyn_ntags::Vector{Int} = Int[]
    xpyn_ncoor::Vector{Float64} = Float64[]
    xpyp_ntags::Vector{Int} = Int[]
    xpyp_ncoor::Vector{Float64} = Float64[]
    xnyp_ntags::Vector{Int} = Int[]
    xnyp_ncoor::Vector{Float64} = Float64[]
    #
    xpypzp_ntags::Vector{Int} = Int[]
    xpypzp_ncoor::Vector{Float64} = Float64[]
    xpypzn_ntags::Vector{Int} = Int[]
    xpypzn_ncoor::Vector{Float64} = Float64[]
    xpynzn_ntags::Vector{Int} = Int[]
    xpynzn_ncoor::Vector{Float64} = Float64[]
    xpynzp_ntags::Vector{Int} = Int[]
    xpynzp_ncoor::Vector{Float64} = Float64[]
    xnypzp_ntags::Vector{Int} = Int[]
    xnypzp_ncoor::Vector{Float64} = Float64[]
    xnypzn_ntags::Vector{Int} = Int[]
    xnypzn_ncoor::Vector{Float64} = Float64[]
    xnynzn_ntags::Vector{Int} = Int[]
    xnynzn_ncoor::Vector{Float64} = Float64[]
    xnynzp_ntags::Vector{Int} = Int[]
    xnynzp_ncoor::Vector{Float64} = Float64[]
    #
    outer_ntags::Vector{Int} = Int[]
    outer_ncoor::Vector{Float64} = Float64[]
    #
    # Making the groups of nodes
    for (k, anode_tag) in enumerate(node_tags)
        nx, ny, nz = node_coordinates[:, k]
        xlb_node = x_min - small_par < nx < x_min + small_par
        ylb_node = y_min - small_par < ny < y_min + small_par
        zlb_node = z_min - small_par < nz < z_min + small_par
        xub_node = x_max - small_par < nx < x_max + small_par
        yub_node = y_max - small_par < ny < y_max + small_par
        zub_node = z_max - small_par < nz < z_max + small_par
        # a_node_coordinates = (nx, ny, nz)
        if (xlb_node || xub_node || ylb_node || yub_node || zlb_node || zub_node)
            #
            push!(outer_ntags, anode_tag)
            push!(outer_ncoor, nx, ny, nz)
            # only faces
            if (xlb_node && !ylb_node && !zlb_node && !yub_node && !zub_node)
                push!(xn_ntags, anode_tag)
                push!(xn_ncoor, nx, ny, nz)  # exclusively xn
            elseif (xub_node && !ylb_node && !zlb_node && !yub_node && !zub_node)
                push!(xp_ntags, anode_tag)
                push!(xp_ncoor, nx, ny, nz)  # exclusively xp
            elseif (ylb_node && !xlb_node && !zlb_node && !xub_node && !zub_node)
                push!(yn_ntags, anode_tag)
                push!(yn_ncoor, nx, ny, nz)  # exclusively yn
            elseif (yub_node && !xlb_node && !zlb_node && !xub_node && !zub_node)
                push!(yp_ntags, anode_tag)
                push!(yp_ncoor, nx, ny, nz)  # exclusively yp
            elseif (zlb_node && !xlb_node && !ylb_node && !xub_node && !yub_node)
                push!(zn_ntags, anode_tag)
                push!(zn_ncoor, nx, ny, nz)  # exclusively zn
            elseif (zub_node && !xlb_node && !ylb_node && !xub_node && !yub_node)
                push!(zp_ntags, anode_tag)
                push!(zp_ncoor, nx, ny, nz)  # exclusively zp
            #
            # only edges
            elseif (ylb_node && zlb_node && !xlb_node && !xub_node)
                push!(ynzn_ntags, anode_tag)
                push!(ynzn_ncoor, nx, ny, nz)  # exclusively ynzn
            elseif (yub_node && zlb_node && !xlb_node && !xub_node)
                push!(ypzn_ntags, anode_tag)
                push!(ypzn_ncoor, nx, ny, nz)  # exclusively ypzn
            elseif (yub_node && zub_node && !xlb_node && !xub_node)
                push!(ypzp_ntags, anode_tag)
                push!(ypzp_ncoor, nx, ny, nz)  # exclusively ypzp
            elseif (zub_node && ylb_node && !xlb_node && !xub_node)
                push!(ynzp_ntags, anode_tag)
                push!(ynzp_ncoor, nx, ny, nz)  # exclusively ynzp
            #
            elseif (xlb_node && zlb_node && !ylb_node && !yub_node)
                push!(znxn_ntags, anode_tag)
                push!(znxn_ncoor, nx, ny, nz)  # exclusively znxn
            elseif (xlb_node && zub_node && !ylb_node && !yub_node)
                push!(zpxn_ntags, anode_tag)
                push!(zpxn_ncoor, nx, ny, nz)  # exclusively zpxn
            elseif (xub_node && zub_node && !ylb_node && !yub_node)
                push!(zpxp_ntags, anode_tag)
                push!(zpxp_ncoor, nx, ny, nz)  # exclusively zpxp
            elseif (xub_node && zlb_node && !ylb_node && !yub_node)
                push!(znxp_ntags, anode_tag)
                push!(znxp_ncoor, nx, ny, nz)  # exclusively znxp
            #
            elseif (xlb_node && ylb_node && !zlb_node && !zub_node)
                push!(xnyn_ntags, anode_tag)
                push!(xnyn_ncoor, nx, ny, nz)  # exclusively xnyn
            elseif (xub_node && ylb_node && !zlb_node && !zub_node)
                push!(xpyn_ntags, anode_tag)
                push!(xpyn_ncoor, nx, ny, nz)  # exclusively xpyn
            elseif (xub_node && yub_node && !zlb_node && !zub_node)
                push!(xpyp_ntags, anode_tag)
                push!(xpyp_ncoor, nx, ny, nz)  # exclusively xpyp
            elseif (xlb_node && yub_node && !zlb_node && !zub_node)
                push!(xnyp_ntags, anode_tag)
                push!(xnyp_ncoor, nx, ny, nz)  # exclusively xnyp
            #
            # only vertices
            elseif (xub_node && yub_node && zub_node)
                push!(xpypzp_ntags, anode_tag)
                push!(xpypzp_ncoor, nx, ny, nz)  # exclusively xpypzp
            elseif (xub_node && yub_node && zlb_node)
                push!(xpypzn_ntags, anode_tag)
                push!(xpypzn_ncoor, nx, ny, nz)  # exclusively xpypzn
            elseif (xub_node && ylb_node && zlb_node)
                push!(xpynzn_ntags, anode_tag)
                push!(xpynzn_ncoor, nx, ny, nz)  # exclusively xpynzn
            elseif (xub_node && ylb_node && zub_node)
                push!(xpynzp_ntags, anode_tag)
                push!(xpynzp_ncoor, nx, ny, nz)  # exclusively xpynzp
            elseif (xlb_node && yub_node && zub_node)
                push!(xnypzp_ntags, anode_tag)
                push!(xnypzp_ncoor, nx, ny, nz)  # exclusively xnypzp
            elseif (xlb_node && yub_node && zlb_node)
                push!(xnypzn_ntags, anode_tag)
                push!(xnypzn_ncoor, nx, ny, nz)  # exclusively xnypzn
            elseif (xlb_node && ylb_node && zlb_node)
                push!(xnynzn_ntags, anode_tag)
                push!(xnynzn_ncoor, nx, ny, nz)  # exclusively xnynzn
            elseif (xlb_node && ylb_node && zub_node)
                push!(xnynzp_ntags, anode_tag)
                push!(xnynzp_ncoor, nx, ny, nz)  # exclusively xnynzp
            end
        end
    end
    #
    return Dict(
        "xp" => IF64[reshape(xp_ntags, 1, :); reshape(xp_ncoor, 3, :)],
        "xn" => IF64[reshape(xn_ntags, 1, :); reshape(xn_ncoor, 3, :)],
        "yp" => IF64[reshape(yp_ntags, 1, :); reshape(yp_ncoor, 3, :)],
        "yn" => IF64[reshape(yn_ntags, 1, :); reshape(yn_ncoor, 3, :)],
        "zp" => IF64[reshape(zp_ntags, 1, :); reshape(zp_ncoor, 3, :)],
        "zn" => IF64[reshape(zn_ntags, 1, :); reshape(zn_ncoor, 3, :)],
        "ypzn" => IF64[reshape(ypzn_ntags, 1, :); reshape(ypzn_ncoor, 3, :)],
        "ynzp" => IF64[reshape(ynzp_ntags, 1, :); reshape(ynzp_ncoor, 3, :)],
        "ypzp" => IF64[reshape(ypzp_ntags, 1, :); reshape(ypzp_ncoor, 3, :)],
        "ynzn" => IF64[reshape(ynzn_ntags, 1, :); reshape(ynzn_ncoor, 3, :)],
        "znxp" => IF64[reshape(znxp_ntags, 1, :); reshape(znxp_ncoor, 3, :)],
        "zpxn" => IF64[reshape(zpxn_ntags, 1, :); reshape(zpxn_ncoor, 3, :)],
        "zpxp" => IF64[reshape(zpxp_ntags, 1, :); reshape(zpxp_ncoor, 3, :)],
        "znxn" => IF64[reshape(znxn_ntags, 1, :); reshape(znxn_ncoor, 3, :)],
        "xnyp" => IF64[reshape(xnyp_ntags, 1, :); reshape(xnyp_ncoor, 3, :)],
        "xpyn" => IF64[reshape(xpyn_ntags, 1, :); reshape(xpyn_ncoor, 3, :)],
        "xpyp" => IF64[reshape(xpyp_ntags, 1, :); reshape(xpyp_ncoor, 3, :)],
        "xnyn" => IF64[reshape(xnyn_ntags, 1, :); reshape(xnyn_ncoor, 3, :)],
        "xpynzn" => IF64[reshape(xpynzn_ntags, 1, :); reshape(xpynzn_ncoor, 3, :)],
        "xnypzn" => IF64[reshape(xnypzn_ntags, 1, :); reshape(xnypzn_ncoor, 3, :)],
        "xnynzp" => IF64[reshape(xnynzp_ntags, 1, :); reshape(xnynzp_ncoor, 3, :)],
        "xnypzp" => IF64[reshape(xnypzp_ntags, 1, :); reshape(xnypzp_ncoor, 3, :)],
        "xpynzp" => IF64[reshape(xpynzp_ntags, 1, :); reshape(xpynzp_ncoor, 3, :)],
        "xpypzn" => IF64[reshape(xpypzn_ntags, 1, :); reshape(xpypzn_ncoor, 3, :)],
        "xpypzp" => IF64[reshape(xpypzp_ntags, 1, :); reshape(xpypzp_ncoor, 3, :)],
        "xnynzn" => IF64[reshape(xnynzn_ntags, 1, :); reshape(xnynzn_ncoor, 3, :)],
        "outer_nodes" => IF64[reshape(outer_ntags, 1, :); reshape(outer_ncoor, 3, :)],
    )
end

#
function get_rve_node_groups(
    node_tags::Vector{Int},
    node_coordinates::Matrix{Float64},
    rve_bbox::NTuple{4,Float64};
    small_par::Float64=1e-06
)::Dict{String,Matrix{IF64}}
    x_min, y_min, x_max, y_max = rve_bbox
    xlb_node::Bool, ylb_node::Bool = false, false
    xub_node::Bool, yub_node::Bool = false, false
    #
    a_node_coordinates::NTuple{3,Float64} = (0.0, 0.0, 0.0)
    xn_ntags::Vector{Int} = Int[]
    xn_ncoor::Vector{Float64} = Float64[]
    yn_ntags::Vector{Int} = Int[]
    yn_ncoor::Vector{Float64} = Float64[]
    xp_ntags::Vector{Int} = Int[]
    xp_ncoor::Vector{Float64} = Float64[]
    yp_ntags::Vector{Int} = Int[]
    yp_ncoor::Vector{Float64} = Float64[]
    #
    # Vertices
    xnyn_ntags::Vector{Int} = Int[]
    xnyn_ncoor::Vector{Float64} = Float64[]
    xpyn_ntags::Vector{Int} = Int[]
    xpyn_ncoor::Vector{Float64} = Float64[]
    xpyp_ntags::Vector{Int} = Int[]
    xpyp_ncoor::Vector{Float64} = Float64[]
    xnyp_ntags::Vector{Int} = Int[]
    xnyp_ncoor::Vector{Float64} = Float64[]
    #
    outer_ntags::Vector{Int} = Int[]
    outer_ncoor::Vector{Float64} = Float64[]
    #
    # Making the groups of nodes
    for (k, anode_tag) in enumerate(node_tags)
        nx, ny = node_coordinates[:, k]
        xlb_node = x_min - small_par < nx < x_min + small_par
        ylb_node = y_min - small_par < ny < y_min + small_par
        xub_node = x_max - small_par < nx < x_max + small_par
        yub_node = y_max - small_par < ny < y_max + small_par
        # a_node_coordinates = (nx, ny, nz)
        if (xlb_node || xub_node || ylb_node || yub_node)
            #
            push!(outer_ntags, anode_tag)
            push!(outer_ncoor, nx, ny,)
            # only edges
            if (xlb_node && !ylb_node && !yub_node)   # exclusively xn
                push!(xn_ntags, anode_tag)
                push!(xn_ncoor, nx, ny,)
            elseif (xub_node && !ylb_node && !yub_node)  # exclusively xp
                push!(xp_ntags, anode_tag)
                push!(xp_ncoor, nx, ny,)
            elseif (ylb_node && !xlb_node && !xub_node)  # exclusively yn
                push!(yn_ntags, anode_tag)
                push!(yn_ncoor, nx, ny,)
            elseif (yub_node && !xlb_node && !xub_node)  # exclusively yp
                push!(yp_ntags, anode_tag)
                push!(yp_ncoor, nx, ny,)
            #
            # only vertices
            elseif (xlb_node && ylb_node)  # exclusively xnyn
                push!(xnyn_ntags, anode_tag)
                push!(xnyn_ncoor, nx, ny,)
            elseif (xub_node && ylb_node)  # exclusively xpyn
                push!(xpyn_ntags, anode_tag)
                push!(xpyn_ncoor, nx, ny,)
            elseif (xub_node && yub_node)  # exclusively xpyp
                push!(xpyp_ntags, anode_tag)
                push!(xpyp_ncoor, nx, ny,)
            elseif (xlb_node && yub_node)  # exclusively xnyp
                push!(xnyp_ntags, anode_tag)
                push!(xnyp_ncoor, nx, ny,)
            end
        end
    end
    #
    return Dict(
        "xp" => IF64[reshape(xp_ntags, 1, :); reshape(xp_ncoor, 2, :)],
        "xn" => IF64[reshape(xn_ntags, 1, :); reshape(xn_ncoor, 2, :)],
        "yp" => IF64[reshape(yp_ntags, 1, :); reshape(yp_ncoor, 2, :)],
        "yn" => IF64[reshape(yn_ntags, 1, :); reshape(yn_ncoor, 2, :)],
        "xnyp" => IF64[reshape(xnyp_ntags, 1, :); reshape(xnyp_ncoor, 2, :)],
        "xpyn" => IF64[reshape(xpyn_ntags, 1, :); reshape(xpyn_ncoor, 2, :)],
        "xpyp" => IF64[reshape(xpyp_ntags, 1, :); reshape(xpyp_ncoor, 2, :)],
        "xnyn" => IF64[reshape(xnyn_ntags, 1, :); reshape(xnyn_ncoor, 2, :)],
        "outer_nodes" => IF64[reshape(outer_ntags, 1, :); reshape(outer_ncoor, 2, :)],
    )
end

""" Created Vector of Node Pairs for node groups in 3D """
function node_pairs(
    png::Matrix{IF64},  # png: parent node group
    cng::Matrix{IF64},  # cng: child node group
    shift::NTuple{3,Float64},
    pair_tag::String;
    ϵ::Float64=1e-06
)::Vector{FENodePair}
    δx, δy, δz = shift
    node_pairs::Vector{FENodePair} = []
    #
    parent_id, child_id = split(pair_tag, "_")
    # rounding to a place higher than tolerance order. This is to avoid the influence of machine
    # precision on the sorting process.
    png[2:end, :] = round.(png[2:end, :]; digits=(1 + ceil(Int, log10(1 / ϵ))))  # rounding to a place higher than tolerance order
    cng[2:end, :] = round.(cng[2:end, :]; digits=(1 + ceil(Int, log10(1 / ϵ))))
    #
    sorted_png = sortslices(png, dims=2, by=x -> (x[2], x[3], x[4]), rev=false)
    sorted_cng = sortslices(cng, dims=2, by=x -> (x[2], x[3], x[4]), rev=false)
    if size(sorted_png, 2) != size(sorted_cng, 2)
        throw(FEPrepError("Unequal number of nodes on opposite groups"))
    end
    num_nodes = size(png, 2)
    for i = 1:num_nodes
        apg_nt, apg_nx, apg_ny, apg_nz = sorted_png[:, i]
        acg_nt, acg_nx, acg_ny, acg_nz = sorted_cng[:, i]  # TODO check if vectorization increases performance
        if (
            abs(abs(apg_nx - acg_nx) - δx) < ϵ &&
            abs(abs(apg_ny - acg_ny) - δy) < ϵ &&
            abs(abs(apg_nz - acg_nz) - δz) < ϵ
        )
            push!(node_pairs, FENodePair("$parent_id-$apg_nt-$child_id-$acg_nt", apg_nt, acg_nt))
        else
            throw(MeshPeriodicityError(ang1_nt, pair_tag))
        end
    end
    return node_pairs
end

""" Created Vector of Node Pairs for node groups in 2D """
function node_pairs(
    png::Matrix{IF64},  # png: parent node group
    cng::Matrix{IF64},  # cng: child node group
    shift::NTuple{2,Float64},
    pair_tag::String;
    ϵ::Float64=1e-06
)::Vector{FENodePair}
    δx, δy = shift
    node_pairs::Vector{FENodePair} = []
    #
    parent_id, child_id = split(pair_tag, "_")
    # rounding to a place higher than tolerance order. This is to avoid the influence of machine
    # precision on the sorting process.
    png[2:end, :] = round.(png[2:end, :]; digits=(1 + ceil(Int, log10(1 / ϵ))))  # rounding to a place higher than tolerance order
    cng[2:end, :] = round.(cng[2:end, :]; digits=(1 + ceil(Int, log10(1 / ϵ))))
    #
    sorted_png = sortslices(png, dims=2, by=x -> (x[2], x[3]), rev=false)
    sorted_cng = sortslices(cng, dims=2, by=x -> (x[2], x[3]), rev=false)
    if size(sorted_png, 2) != size(sorted_cng, 2)
        throw(FEPrepError("Unequal number of nodes on opposite groups"))
    end
    num_nodes = size(png, 2)
    for i = 1:num_nodes
        apg_nt, apg_nx, apg_ny = sorted_png[:, i]
        acg_nt, acg_nx, acg_ny = sorted_cng[:, i]  # TODO check if vectorization increases performance
        if (
            abs(abs(apg_nx - acg_nx) - δx) < ϵ &&
            abs(abs(apg_ny - acg_ny) - δy) < ϵ
        )
            push!(node_pairs, FENodePair("$parent_id-$apg_nt-$child_id-$acg_nt", apg_nt, acg_nt))
        else
            throw(MeshPeriodicityError(ang1_nt, pair_tag))
        end
    end
    return node_pairs
end

function node_pairs(
    node_groups::Dict{String,Matrix{IF64}},
    bbox::NTuple{6,Float64},
)::Vector{FENodePair}
    x_min, y_min, z_min, x_max, y_max, z_max = bbox
    lx, ly, lz = x_max - x_min, y_max - y_min, z_max - z_min
    #
    parent_child_map = Dict(
            "xn" => [("xp", (lx, 0.0, 0.0))],
            "yn" => [("yp", (0.0, ly, 0.0))],
            "zn" => [("zp", (0.0, 0.0, lz))],
            "xnyn" => [
                ("xnyp", (0.0, ly, 0.0)),
                ("xpyn", (lx, 0.0, 0.0)),
                ("xpyp", (lx, ly, 0.0)),
            ],
            "ynzn" => [
                ("ynzp", (0.0, 0.0, lz)),
                ("ypzn", (0.0, ly, 0.0)),
                ("ypzp", (0.0, ly, lz)),
            ],
            "znxn" => [
                ("znxp", (lx, 0.0, 0.0)),
                ("zpxn", (0.0, 0.0, lz)),
                ("zpxp", (lx, 0.0, lz)),
            ],
            "xnynzn" => [
                ("xnynzp", (0.0, 0.0, lz)),
                ("xnypzn", (0.0, ly, 0.0)),
                ("xnypzp", (0.0, ly, lz)),
                ("xpynzn", (lx, 0.0, 0.0)),
                ("xpynzp", (lx, 0.0, lz)),
                ("xpypzn", (lx, ly, 0.0)),
                ("xpypzp", (lx, ly, lz)),
            ],
        )
    fe_node_pairs::Vector{FENodePair} = FENodePair[]
    for (a_parent, children) in parent_child_map
        for (a_child, distance) in children
            if (a_parent in keys(node_groups) &&
                a_child in keys(node_groups))
                append!(
                    fe_node_pairs,
                    node_pairs(
                        node_groups[a_parent],
                        node_groups[a_child],
                        distance,
                        a_parent * "_" * a_child,
                    ),
                )
            end
        end
    end
    return fe_node_pairs
end

function node_pairs(
    node_groups::Dict{String,Matrix{IF64}},
    bbox::NTuple{4,Float64},
)::Vector{FENodePair}
    x_min, y_min, x_max, y_max= bbox
    lx, ly= x_max - x_min, y_max - y_min
    #
    parent_child_map = Dict(
            "xn" => [("xp", (lx, 0.0))],
            "yn" => [("yp", (0.0, ly))],
            "xnyn" => [
                ("xnyp", (0.0, ly)),
                ("xpyn", (lx, 0.0)),
                ("xpyp", (lx, ly)),
            ],
        )
    fe_node_pairs::Vector{FENodePair} = FENodePair[]
    for (a_parent, children) in parent_child_map
        for (a_child, distance) in children
            if (a_parent in keys(node_groups) &&
                a_child in keys(node_groups))
                append!(
                    fe_node_pairs,
                    node_pairs(
                        node_groups[a_parent],
                        node_groups[a_child],
                        distance,
                        a_parent * "_" * a_child,
                    ),
                )
            end
        end
    end
    return fe_node_pairs
end

"""
    make_rve_node_pairs(ntags, ncoor)

 Returns a Vector of finite element node pairs (FENodePair).

"""
function make_rve_node_pairs(
    node_tags::Vector{Int},
    node_coordinates::Matrix{Float64};
    small_para::Float64=1e-06
)::Vector{FENodePair}
    if length(node_tags) != size(node_coordinates, 2)
        throw(
            DimensionMismatch(
                "Number of node tags and number of columns in coordinate materix are not equal. Please ensure that, size of the coordinate matrix is (3, num_tags).
                Ensure that each column of the coordinate matrix represents a node coordinates.
                This is to take the advantage of column major representation of Julia.",
            ),
        )
    end
    ((x_min, x_max), (y_min, y_max), (z_min, z_max)) = extrema(node_coordinates, dims=2)
    rve_bbox::NTuple{6,Float64} = (x_min, y_min, z_min, x_max, y_max, z_max)
    rve_dim_tag = get_rve_dim_tag(rve_bbox)
    if rve_dim_tag == "XYZ"
        rve_node_groups = get_rve_node_groups(node_tags, node_coordinates, rve_bbox; small_par=small_para)
        return node_pairs(rve_node_groups, rve_bbox)
    elseif rve_dim_tag == "XY"
        bbox = (x_min, y_min, x_max, y_max)
        rve_node_groups = get_rve_node_groups(node_tags, node_coordinates, bbox; small_par=small_para)
        return node_pairs(rve_node_groups, bbox)
    end
end

# ==================================================================================#
#
#                       Prepare Finite Element Sets
#
# ==================================================================================#

"""
    get_felement_sets() -> Vector{AbstractFElementSet}

 Args:
 -----
 - **element_connectivity**::Dict{Int64, Matrix{Int64}},
     Element conne gmsh ids are used to identify element sets which will eventually be converted to ABAQUS convention of element naming.
 - **nodal_data**::Dict{Int64, Tuple{Vararg{Float64}}},
 - **elset_tag**::Union{Int64, String};
 - **num_ip**::Int64=-1,
 - **gq_order**::Int64=-1,
 - **print_debug_info**::Bool=false,

"""

function make_finite_element_set(
    element_type::Int,
    element_connectivity::Matrix{Int},
    nodal_data::Dict{Int,Vector{Float64}};
    elset_tag::Union{String,Int}="Element-Set",
    material::Union{Nothing,Material}=nothing,
    num_ip::Int64=-1,
    gq_order::Int64=-1
)::FiniteElementSet
    #
    elements::Vector{AbstractFElement} = AbstractFElement[]
    ivols::Vector{Float64} = Float64[]
    shfg_xyz::Vector{SMatrix} = SMatrix[]
    num_node_per_fele::Int = size(element_connectivity, 1) - 1  # first row is for ele_tag => size-1
    # num_fele::Int = size(element_connectivity, 2)
    num_shape_func::Int = num_node_per_fele
    dim, fele_data_type, integration_points = get_felement_properties(
        element_type, gq_order, num_ip
    )
    parent_csys_shf_data = get_shape_functions_info(
        fele_data_type,
        integration_points
    )
    #
    # Running over each element of the connectivity table.
    for a_ele_info in eachcol(element_connectivity)
        empty!(ivols)
        empty!(shfg_xyz)
        ele_tag = a_ele_info[1]
        ele_nodes = a_ele_info[2:end]
        nodal_coordinates = [nodal_data[i][1:dim] for i in ele_nodes]
        for (a_ip, shfg_parent) in parent_csys_shf_data
            Jacobian = zeros(Float64, dim, dim)
            # run over pairs of coordinates and shape function data
            for (k2, a_node_coor) in enumerate(nodal_coordinates)
                Jacobian +=
                    (shfg_parent[2:end, k2:k2] * reshape(collect(a_node_coor), 1, :))
            end
            Jacobian_det = det(Jacobian)
            @assert Jacobian_det > 0.0 "Determinant of a Jacobian must be >= 0.0 but the element $(ele_tag) has $Jacobian_det."
            push!(ivols, Jacobian_det * a_ip.wt)
            aip_shf_grad_xyz =
                [shfg_parent[1:1, :]; inv(Jacobian) * shfg_parent[2:end, :]]
            push!(shfg_xyz, SMatrix{1 + dim,num_shape_func,Float64}(aip_shf_grad_xyz))
        end
        push!(
            elements,
            fele_data_type(
                ele_tag,
                SVector{num_node_per_fele,Int}(ele_nodes),
                tuple(ivols...),
                tuple(shfg_xyz...),
            )
        )  # end of updating elements vector.
    end
    return FiniteElementSet(
        elset_tag,
        elements,
        material,
    )
end

function make_finite_element_sets(
    element_connectivities::Dict{Int64,Matrix{Int64}},
    nodal_data::Dict{Int64,Vector{Float64}},
    ele_material::Material,
    elset_tag::Union{Int64,String};
    num_ip::Int64=-1,
    gq_order::Int64=-1
)::Vector{FiniteElementSet}
    #
    return [make_finite_element_set(
        a_ele_type,
        a_ele_type_ele_connectivity,
        nodal_data;
        elset_tag=elset_tag,
        material=ele_material,
        num_ip=num_ip,
        gq_order=gq_order
    ) for (a_ele_type, a_ele_type_ele_connectivity) in element_connectivities]
end

# =================================================================================
#                       Prepare Node Sets
# =================================================================================

"""
It returns a dictionary with the following key-value pairs

key:
----
    Node tag
Value:
-----
    A tuple of nodal coordinates. If not specified, it will have x, y, z coordinates but one can choose to have either (x, y) or (y, z) or (z, x) or (x, ) or (y, ) or (z, ) using ** option

"""
function make_fe_node_set(
    node_tags::Union{Vector{Int64},Vector{Int32}},
    nodal_coor::Matrix{Float64};
    nodal_dof::Tuple{Vararg{String}}=("x", "y", "z")
)::Dict{Int64,Vector{Float64}}  # TODO make node sets object oriented
    #
    @assert size(nodal_coor, 1) >= 2 "Number of rows must be 2 or 3 being in the order of x, y, z"
    @assert length(node_tags) == size(nodal_coor, 2) "Number of node tags and number of columns in nodal coordinates matrix doesn't match; $(length(node_tags)) ≠ $(size(nodal_coor)[1]))"
    #
    req_dofs::Vector{Int64} = []
    for a_dof in nodal_dof
        if uppercase(a_dof) == "X"
            push!(req_dofs, 1)
        elseif uppercase(a_dof) == "Y"
            push!(req_dofs, 2)
        elseif uppercase(a_dof) == "Z"
            push!(req_dofs, 3)
        end
    end
    #
    nodal_coor = nodal_coor[req_dofs, :]  # extracting dof, as required.
    #
    node_tags = convert(Vector{Int64}, node_tags)
    #
    return Dict((i, j) for (i, j) in zip(node_tags, eachcol(nodal_coor)))
end

# =================================================================================
#                       Misc.
# =================================================================================

function get_total_volume_of_elements(
    finite_element_sets::Vector{FiniteElementSet};
    print_info::Bool=false
)
    ele_set_volumes::Vector{Float64} = Float64[]
    element_counter::Int = 0
    a_ele_set_vol::Float64 = 0.0
    if print_info
        println("\n EVALUATING VOLUMES OF ELEMENT SETS...!")
    end
    for a_finite_element_set in finite_element_sets
        a_ele_set_vol = 0.0
        for a_finte_element in a_finite_element_set.elements
            a_ele_set_vol += sum(a_finte_element.ip_volumes)
            element_counter += 1
        end
        if print_info
            println(
                "Volume of element set with ID: ", a_finite_element_set.tag, " is ", a_ele_set_vol
            )
        end
        push!(ele_set_volumes, a_ele_set_vol)
    end
    total_volume::Float64 = sum(ele_set_volumes)
    if print_info
        println("Sum of volumes of ", element_counter, " elements is  ", total_volume)
    end
    return total_volume
end

# ================================================
#                   Archives
# ================================================

# function get_felement_sets_a(
#     element_connectivity::Dict{Int64,Matrix{Int64}},
#     nodal_data::Dict{Int64,Vector{Float64}},
#     elset_tag::Union{Int64,String};
#     num_ip::Int64 = -1,
#     gq_order::Int64 = -1,
#     print_debug_info::Bool = false
# )::Vector{AbstractFElementSet}
#     #
#     println("Preparing $elset_tag...!")
#     element_sets::Vector{AbstractFElementSet} = []
#     elements::Vector{AbstractFElement} = []
#     num_nodes::Int64 = 0
#     num_shf::Int64 = 0
#     num_felements::Int64 = 0
#     ivols::Vector{Float64} = Float64[]
#     shfg_xyz::Vector{SMatrix} = SMatrix[]
#     ele_tag::Int = 0
#     ele_nodes::Vector{Int} = Int[]
#     #
#     # TODO add element type assertion
#     tot_vol::Float64 = 0.0   # *****************
#     #
#     # Running over each element type, with their element_connectivity information.
#     for (ele_type, ele_con) in element_connectivity
#         empty!(elements)
#         num_nodes = size(ele_con)[1] - 1  # as the first row corresponds to element tag
#         num_felements = size(ele_con)[2]
#         num_shf = num_nodes
#         dim, fele_data_type, integration_points = get_felement_properties(
#             ele_type,
#             gq_order,
#             num_ip
#         )
#         #
#         # For the given element type, get shape functions and its gradients at
#         # gauss points in parent coordinate system
#         shf_data_parent_csys = get_shape_functions_info(
#             fele_data_type,
#             integration_points
#         )
#         #
#         if print_debug_info
#             println(
#                 "Element Data Type: $fele_data_type,
#                 Number of elements :: $(num_felements),
#                 Number of nodes/element :: $num_nodes,
#                 Integration points :: $integration_points,",
#             )
#         end
#         #
#         # Running over each element of the given type.
#         for a_col in eachcol(ele_con)
#             empty!(ivols)
#             empty!(shfg_xyz)
#             #
#             ele_tag, ele_nodes = a_col[1], a_col[2:end]
#             nodal_coordinates = (nodal_data[i][1:dim] for i in ele_nodes)  # NTuple{N, Vector{Float64}}
#             #
#             # run over each gauss point and get Jacobian and its determinant
#             for (a_ip, shfg_parent) in shf_data_parent_csys
#                 Jacobian = zeros(Float64, dim, dim)
#                 # run over pairs of coordinates and shape function data
#                 for (k2, a_node_coor) in enumerate(nodal_coordinates)
#                     Jacobian +=
#                         (shfg_parent[2:end, k2:k2] * reshape(collect(a_node_coor), 1, :))
#                 end
#                 Jacobian_det = det(Jacobian)

#                 @assert Jacobian_det > 0.0 "Determinant of a Jacobian must be >= 0.0 but the element $(ele_tag) has $Jacobian_det."

#                 push!(ivols, Jacobian_det * a_ip.wt)
#                 #
#                 aip_shf_grad_xyz =
#                     [shfg_parent[1:1, :]; inv(Jacobian) * shfg_parent[2:end, :]]
#                 push!(shfg_xyz, SMatrix{1 + dim,num_shf,Float64}(aip_shf_grad_xyz))
#             end
#             elements = [
#                 elements
#                 fele_data_type(
#                     tag = ele_tag,
#                     node_tags = SVector{num_nodes}(ele_nodes),
#                     ip_volumes = ivols,
#                     sfgrad_xyz = shfg_xyz,
#                 )
#             ]
#             tot_vol += sum(ivols)  # *****************
#         end
#         #
#         push!(element_sets, FiniteElementSet(elset_tag, elements))
#     end
#     println("Total Volume by brute sum ", tot_vol)   # *****************
#     return element_sets
# end

# function get_felement_sets(
#     element_connectivity::Dict{Int64,Matrix{Int64}},
#     nodal_data::Dict{Int64, Vector{Float64}},
#     elset_tag::Union{Int64,String};
#     num_ip::Int64 = -1,
#     gq_order::Int64 = -1,
#     print_debug_info::Bool = false,
# )::Vector{AbstractFElementSet}
#     #
#     println("Preparing $elset_tag...!")
#     element_sets::Vector{AbstractFElementSet} = AbstractFElementSet[]
#     elements::Vector{AbstractFElement} = AbstractFElement[]
#     num_nodes::Int64 = 0
#     num_shf::Int64 = 0
#     num_felements::Int64 = 0
#     ivols::Vector{Float64} = []
#     shfg_xyz::Vector{SMatrix} = []
#     #
#     #
#     for (ele_type, ele_con) in element_connectivity
#         # @assert ele_type in fe_libray_g "Element type $ele_type (GMSH id) is not implemented as of now..!"
#         empty!(elements)
#         num_nodes = size(ele_con)[1] - 1  # as the first row corresponds to element tag
#         num_felements = size(ele_con)[2]
#         dim, fele_data_type, integration_points =
#             get_felement_properties(ele_type, gq_order, num_ip)
#         num_ip = length(integration_points)
#         num_shf = num_nodes
#         if print_debug_info
#             println(
#                 "For $fele_data_type, \nNumber of elements :: $(num_felements), \nNumber of nodes/element :: $num_nodes, \nIntegration points :: $integration_points,",
#             )
#         end
#         #
#         # get shape functions and its gradients at gauss points in parent coordinate system
#         shf_data_parent_csys = get_shape_functions_info(fele_data_type, integration_points)

#         # get Jacobian and its determinant
#         for (ele_tag, ele_nodes...) in eachcol(ele_con)
#             nodal_coordinates = (nodal_data[i][1:dim] for i in ele_nodes)
#             empty!(ivols)
#             empty!(shfg_xyz)
#             # run over each gauss point
#             for (a_ip, shfg_parent) in shf_data_parent_csys
#                 Jacobian = zeros(Float64, dim, dim)
#                 # run over pairs of coordinates and shape function data

#                 for (k2, a_node_coor) in enumerate(nodal_coordinates)
#                     Jacobian +=
#                         (shfg_parent[2:end, k2:k2] * reshape(collect(a_node_coor), 1, :))
#                 end

#                 Jacobian_det = det(Jacobian)

#                 @assert Jacobian_det > 0.0 "Determinant of a Jacobian must be >= 0.0 but $Jacobian_det observed on element $(ele_tag)"

#                 push!(ivols, Jacobian_det * a_ip.wt)
#                 #
#                 aip_shf_grad_xyz =
#                     [shfg_parent[1:1, :]; inv(Jacobian) * shfg_parent[2:end, :]]
#                 push!(shfg_xyz, SMatrix{1 + dim,num_shf,Float64}(aip_shf_grad_xyz))
#             end
#             elements = [
#                 elements
#                 fele_data_type(
#                     tag = ele_tag,
#                     node_tags = SVector(ele_nodes...),
#                     ip_volumes = ivols,
#                     sfgrad_xyz = shfg_xyz,
#                 )
#             ]
#         end
#         #
#         push!(element_sets, FiniteElementSet(elset_tag, elements))
#     end
#     return element_sets
# end


text:
using LaTeXStrings
import JLD
import PyPlot

include("utilities.jl")

times = Dict()
chs = Dict()

#do the ugta plot
network_dir = "TSA250_50"
fig, ax = PyPlot.subplots()
resultsfilename = "results/" * replace(network_dir, "/"=>"_") * "_results.jld"
times, chs = JLD.load(resultsfilename, "times", "chs")
names = [map(x->x[1][1] * "+" * x[2][1], Base.product(preconditioned_solvers, preconditioners))[:]; map(x->x[1], unpreconditioned_solvers)]
names = sort(names; by=n->times[n])
ax.barh(names, map(x->times[x], names), color=map(n->ifelse(chs[n][end] < 1e-10, "C0", "C1"), names))
c0_patch = PyPlot.matplotlib.patches.Patch(color="C0", label="converged")
c1_patch = PyPlot.matplotlib.patches.Patch(color="C1", label="did not converge")
ax.legend(handles=[c0_patch, c1_patch], fontsize=14)
ax.set_xscale("log")
ax.tick_params(axis="x", labelsize=14)
ax.tick_params(axis="y", labelsize=14)
ax.set_xlim(1e2, 3e4)
ax.set_xlabel("Run Time (seconds)", fontsize=16)
fig.tight_layout()
fig.savefig("ugta.pdf")
display(fig)
println()
PyPlot.close(fig)

#do the plot on the grid scaling for the thomas dfns
numnodes = Dict("thomas_tpl_p2_x01"=>62070, "thomas_tpl_p3_x01"=>497635, "thomas_tpl_p5_x01"=>2309464, "thomas_tpl_p10_x01"=>7123318, "thomas_tpl_p20_x01"=>13038898)
network_dirs = ["thomas_tpl_p2_x01", "thomas_tpl_p3_x01", "thomas_tpl_p5_x01", "thomas_tpl_p10_x01", "thomas_tpl_p20_x01"]
fig, ax = PyPlot.subplots()
for (i, network_dir) in enumerate(network_dirs)
	resultsfilename = "results/" * replace(network_dir, "/"=>"_") * "_results.jld"
	thesetimes, thesechs = JLD.load(resultsfilename, "times", "chs")
	times[network_dir] = thesetimes
	chs[network_dir] = thesechs
end
for name in names
	if sum(map(network_dir->chs[network_dir][name][end] < 1e-10, network_dirs)) > 2
		ax.loglog(map(network_dir->numnodes[network_dir], network_dirs), map(network_dir->ifelse(chs[network_dir][name][end] < 1e-10, times[network_dir][name], NaN), network_dirs), ".-", label=name, alpha=0.75, linewidth=3, ms=10)
	end
end
ax.legend(fontsize=14)
ax.tick_params(axis="x", labelsize=14)
ax.tick_params(axis="y", labelsize=14)
ax.set_xlabel("Number of grid nodes", fontsize=16)
ax.set_ylabel("Run time (seconds)", fontsize=16)
fig.tight_layout()
fig.savefig("node_scaling.pdf")
display(fig)
println()
PyPlot.close(fig)

#do the plot on the b scaling for the homogenous dfns
network_dirs = ["var_b_0.0", "var_b_0.5", "var_b_1.0"]
bs = Dict(zip(network_dirs, ["0", "0.5", "1.0"]))
fig, ax = PyPlot.subplots()
for (i, network_dir) in enumerate(network_dirs)
	resultsfilename = "results/" * replace(network_dir, "/"=>"_") * "_results.jld"
	thesetimes, thesechs = JLD.load(resultsfilename, "times", "chs")
	times[network_dir] = thesetimes
	chs[network_dir] = thesechs
end
for name in names
	if sum(map(network_dir->chs[network_dir][name][end] < 1e-10, network_dirs)) > 2
		ax.semilogy(map(network_dir->bs[network_dir], network_dirs), map(network_dir->ifelse(chs[network_dir][name][end] < 1e-10, times[network_dir][name], NaN), network_dirs), ".-", label=name, alpha=0.75, linewidth=3, ms=10)
	end
end
ax.legend(fontsize=14)
ax.tick_params(axis="x", labelsize=14)
ax.tick_params(axis="y", labelsize=14)
ax.set_xlabel("Variance of the log-aperture", fontsize=16)
ax.set_ylabel("Run time (seconds)", fontsize=16)
fig.tight_layout()
fig.savefig("b_scaling.pdf")
display(fig)
println()
PyPlot.close(fig)

ext: jl


text:
function romberg(f, a, b, n)
	r = zeros(2, n)
	h = b - a
	r[1, 1] = (h/2) * (f(a) + f(b))
	for i ∈ 2:n
		tsum = 0
		for k ∈ 1:(2^(i - 2))
			tsum += f(a + (k - 0.5) * h)
		end
		r[2, 1] = (1/2)*(h*tsum + r[1, 1])
		for j ∈ 2:i
			r[2, j] = r[2, j - 1] + (r[2, j - 1] - r[1, j - 1])/(4^(j - 1) - 1)
		end
		h = h/2
		for j ∈ 1:i
			r[1, j] = r[2, j]
		end
	end
	return r[2, n]
end

f(x) = sin(x)
println(romberg(f, 0, 3, 20))


text:
function classify(number)
    factors = []
    if number == 1
        return "deficient"
    elseif number <= 0
        return throw(DomainError(number,"Classification is only possible for natural numbers."))
    end

    for i in 1:(number-1)
        if number%i == 0
            push!(factors, i)
        end
    end
    if sum(factors) > number
        return "abundant"
    elseif sum(factors) < number
        return "deficient"
    else
        return "perfect"
    end
end

max_stars_repo_name: Akshat-mehrotra/julia

max_issues_repo_name: Akshat-mehrotra/julia

max_forks_repo_name: Akshat-mehrotra/julia


text:
module BlossomV

if isfile(joinpath(Pkg.dir("BlossomV"),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("BlossomV not properly installed. Please run Pkg.build(\"BlossomV\")")
end

export
    PerfectMatchingCtx,
    Matching, solve, add_edge, get_match, get_all_matches

type PerfectMatchingCtx
    ptr::Ptr{Void}

    function PerfectMatchingCtx(ptr::Ptr{Void})
        self = new(ptr)
        finalizer(self, destroy)
        self
    end
end

function destroy(matching::PerfectMatchingCtx)
    if matching.ptr == C_NULL
        return
    end
    ccall((:matching_destruct, _jl_blossom5), Void, (Ptr{Void},), matching.ptr)
    matching.ptr = C_NULL
    nothing
end

dense_num_edges(node_num) =  node_num*(node_num-1)÷2

Matching(node_num::Integer) = Matching(node_num, dense_num_edges(node_num))
Matching(node_num::Integer, edge_num_max::Integer) = Matching(Int32(node_num), Int32(edge_num_max))
function Matching(node_num::Int32, edge_num_max::Int32)
    assert(node_num % 2 == 0)
    ptr = ccall((:matching_construct, _jl_blossom5), Ptr{Void}, (Int32, Int32), node_num, edge_num_max)
    PerfectMatchingCtx(ptr)
end

function add_edge(matching::PerfectMatchingCtx, first_node::Integer, second_node::Integer, cost::Integer)
	add_edge(matching, Int32(first_node), Int32(second_node), Int32(cost))
end

function add_edge(matching::PerfectMatchingCtx, first_node::Int32, second_node::Int32, cost::Int32)
    first_node != second_node || error("Can not have an edge between $(first_node) and itself")
    first_node >= 0  || error("first_node less than zero (value: $(first_node)). Indexes are zero-based")
    second_node >= 0  || error("second_node less than zero (value: $(second_node)). Indexes are zero-based")
    cost >= 0  || error("Cost must be positive. edge between $(first_node) and $(second_node) is $cost")

	ccall((:matching_add_edge, _jl_blossom5), Int32, (Ptr{Void}, Int32, Int32, Int32), matching.ptr, first_node, second_node, cost)
end

function solve(matching::PerfectMatchingCtx)
    ccall((:matching_solve, _jl_blossom5), Void, (Ptr{Void},), matching.ptr)
    nothing
end

get_match(matching::PerfectMatchingCtx, node::Integer) = get_match(matching, Int32(node))
function get_match(matching::PerfectMatchingCtx, node::Int32)
    ccall((:matching_get_match, _jl_blossom5), Int32, (Ptr{Void}, Int32), matching.ptr, node)
end

function get_all_matches(m::PerfectMatchingCtx, n_nodes::Integer)
    ret = Matrix{Int32}(2, n_nodes÷2)
    assigned = falses(n_nodes)
    kk = 1
    for ii in 0:n_nodes-1
        assigned[ii+1] && continue
        jj = get_match(m, ii)
        @assert(!assigned[jj+1])
        assigned[ii+1] = true
        assigned[jj+1] = true

        ret[1,kk] = ii
        ret[2,kk] = jj
        kk+=1
    end
    @assert(kk-1 == n_nodes÷2)
    return ret
end
end # module


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

__precompile__()

module MetropolisWithinGibbs
# ----------------------------------------------------------------------------------------------------------------------
# This module define the functions used to estimate state-space models with a Metropolis-Within-Gibbs algorithm.
#
# The Metropolis step of the Metropolis-Within-Gibbs follows closely the algorithm used in An and Schorfheide (2007).
# The major difference between the two of them is that An and Schorfheide (2007) draw the parameters directly with
# the Metropolis algorithm. Here, the supports of most parameters are constrained in some dimension.
# Therefore, the parameters are transformed, so that their support is unconstrained, and then the Metropolis algorithm
# is executed. The sum of the natural logarithm of the Jacobian is used to evaluate correctly the prior densities.
# DSGE software (e.g. YADA) uses similar modifications to deal with constrained parameters.
# ----------------------------------------------------------------------------------------------------------------------

	using DataFrames, Distributions, Distributed, LinearAlgebra;
	local_path = dirname(@__FILE__);

	# -----------------------------------------------------------------------------------------------------------------
	# Types
	# -----------------------------------------------------------------------------------------------------------------

	mutable struct ParSsm{X <: Float64}
		y::Array{Union{X, Missing}, 2}
		d::Union{Array{X, 1}}
		Z::Union{Array{X, 1}, Array{X, 2}}
		R::Union{Array{X, 1}, Array{X, 2}}
		c::Union{Array{X, 1}}
		T::Union{Array{X, 1}, Array{X, 2}}
		Q::Union{Array{X, 1}, Array{X, 2}}
		α¹::Union{Array{X, 1}}
		P¹::Union{Array{X, 1}, Array{X, 2}}
		P̄¹::Union{Array{X, 1}, Array{X, 2}}
		λ::Array{X, 1}
		ρ::Array{X, 1}
		logprior::X
		loglik::X
		logposterior::X
	end

	struct SizeParSsm{X <: Int64}
		d::X
		Z::X
		R::X
		c::X
		T::X
		Q::X
		λ::X
		ρ::X
		θ::X
	end

	struct BoolParSsm{X1 <: BitArray{1}, X2 <: BitArray{2}}
		d::X1
		Z::X2
		R::X2
		c::X1
		T::X2
		Q::X2
		λ::X1
		ρ::X1
	end

	struct PriorOpt{X <: Float64} # the logpdf for λ and ρ are constants
		N::Distributions.Normal{X}
		IG::Distributions.InverseGamma{X} #TruncatedNormal{X}
		λ::X
		ρ::X
	end

	import Base.copy
	Base.copy(x::T) where T = T([getfield(x, k) for k ∈ fieldnames(T)]...)

	#=
	Base.copy(m::ParSsm) 	= ParSsm([ copy(getfield(m, k)) for k = 1:length(fieldnames(m)) ]...);
	Base.copy(m::SizeParSsm) = SizeParSsm([ copy(getfield(m, k)) for k = 1:length(fieldnames(m)) ]...);
	Base.copy(m::BoolParSsm) = BoolParSsm([ copy(getfield(m, k)) for k = 1:length(fieldnames(m)) ]...);
	Base.copy(m::PriorOpt) 	= PriorOpt([ copy(getfield(m, k)) for k = 1:length(fieldnames(m)) ]...);
	=#

	# -----------------------------------------------------------------------------------------------------------------
	# Subroutines
	# -----------------------------------------------------------------------------------------------------------------

	# Estimation
	include("$local_path/subroutines/estimation/kalman_diffuse!.jl");
	include("$local_path/subroutines/estimation/mwg_main.jl");
	include("$local_path/subroutines/estimation/mwg_run.jl");

	# Julia's subroutines
	include("$local_path/subroutines/extra/ex_blkdiag.jl");
	include("$local_path/subroutines/extra/ex_inv.jl");
	include("$local_path/subroutines/extra/ex_ismember.jl");

	# Get
	include("$local_path/subroutines/get/get_logjacobian.jl");
	include("$local_path/subroutines/get/get_par_bound.jl");
	include("$local_path/subroutines/get/get_par_unb.jl");
	include("$local_path/subroutines/get/get_progress.jl");
	include("$local_path/subroutines/get/get_random_disturbance.jl");
	include("$local_path/subroutines/get/get_mwg_jump.jl");

	# Set
	include("$local_path/subroutines/set/set_par.jl");

	# -----------------------------------------------------------------------------------------------------------------
	# Export functions
	# -----------------------------------------------------------------------------------------------------------------

	export kalman_diffuse!, mwg_main, ex_blkdiag, ex_inv, ex_ismember, set_par_fast!, ParSsm, SizeParSsm, BoolParSsm, PriorOpt;
end


text:
function hclust2newick(hc_struct::Hclust{Float64}, flat::Bool=true)

#=
  transliterated from  https://rdrr.io/bioc/ctc/src/R/hc2Newick.R

  propertynames(hc_struct)

  merges  :: Array{Int64,2}
  heights :: Array{Float64,1}
  order   :: Array{Int64,1}
  linkage :: Symbol

  merge   :: Array{Int64,2}   (not plural guess R compatibility v.s. Julia convention)
  height  :: Array{Float64,1} (not plural  also identical with heights)
  method  :: Symbol
  labels  :: UnitRange{Int64}

=#

  dist = 0
  if  isnothing(hc_struct.labels)
    labels = 1:length(hc_struct.order)
  else
    labels = hc_struct.labels
  end

  function putparenthesis(i)
    ## recursive function
    ## i: index of row in hc_struct.merge
    j = hc_struct.merge[i, 1]
    k = hc_struct.merge[i, 2]

    if j < 0
      left = labels[-j]
      if k > 0
        dist = hc_struct.height[i] - hc_struct.height[k]
      else
        dist = hc_struct.height[i]
      end
    else
      left = putparenthesis(j)
    end

    if k < 0
      right = labels[-k]
      if j > 0
        dist = hc_struct.height[i] - hc_struct.height[j]
      else
        dist = hc_struct.height[i]
      end
    else
      right = putparenthesis(k)
    end
    if flat
      return join(["(", left, ":", dist/2, ",", right, ":", dist/2, ")"], "")  # R paste()
    else
      return (left=left, right=right, dist=dist)  # R list() ... mutability wories
      # named tuple? let?  StructArrays?
    end
  end

  n = putparenthesis(first(size(hc_struct.merge)))  # R nrows
  if flat
    n = join([n, ";"], "")  # R paste()
  end

  return n

end;

#= write_newick(io:IO, hc_struct::Hclust{Float64}, flat::Bool=true) =#

function write_newick(io::IO, hc_struct::Hclust{Float64}, flat=true)
	#open('newick.tree', w) do io
    	write(io, hclust2newick(hc_struct, flat))
	#end
end


text:
@doc """
    vardecomp2level(data::Array{Float64, 2}, factor::Array{Float64}, betas::Array{Float64}, factorassign::Array{Float64})

Description:
Compute the portion of the variation of each observable series that may be attributed to their corresponding/assigned latent factors across all levels.

Inputs:
- data = Matrix with each column representing a data series.
- factor = Matrix containing latent factor estimates.
- betas = Matrix containing observation equation coefficient parameter estimates.
- factorassign = Matrix containing the indeces of factors across all levels (columns) assigned to each observable series (rows).

Output:
- vardecomps = Matrix containing the variance contributions of factors across all levels (columns) corresponding to each observable series (rows).
"""
function vardecomp2level(data, factor, betas, factorassign)

    T, N = size(data)
    T, Nf = size(factor)

    # Record data variances
    varvars = zeros(N)
    for i in 1:N
        errors = data[:,i] - betas[i, 2] * factor[:, 1] + betas[i, 3] * factor[:, 1+factorassign[i, 2]]
        varvars[i] = betas[i, 2]^2 * var(factor[:, 1]) + betas[i, 3]^2 * var(factor[:, 1+factorassign[i, 2]]) + var(errors)
    end

    # Record factor estimate variances
    varfacs = zeros(Nf)
    for j in 1:Nf
        varfacs[j] = var(factor[:, j])
    end

    # Compute variance decompositions
    vardecomps = zeros(N, 2)
    for i in 1:N
        vardecomps[i, 1] = (betas[i, 2]^2 * varfacs[1]) / varvars[i]
        vardecomps[i, 2] = (betas[i, 3]^2 * varfacs[1+factorassign[i, 2]]) / varvars[i]
    end

    return vardecomps
end;
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################
######################


text:
#!/usr/bin/env julia

using Distributed

push!(LOAD_PATH, pwd() )

global iterations = 100000
global top_length = 32768

# Create a random chromosome
function random_chromosome(length)

    chromosome = Bool[]

    for i in 1:1:length
        if ( rand(0:1) > 0.5 )
      	    push!(chromosome,true)
        else
	    push!(chromosome,false)
        end
    end
    chromosome
end

# Mutate and compute time
function time_mutations(number,length)
  inicioTiempo = time()

  for i in 0:1:number
      indi = random_chromosome(length)
      ones = @distributed (+) for i in indi
          i
      end
  end
  time()-inicioTiempo

end

function run_onemax()
    length = 16
    while length <= top_length
        println("julia-Onemax-distributed, " * string(length) * ", " * string(time_mutations( iterations,length )))
        length *= 2
    end
end

run_onemax()


text:
function number_of_steps(permanent_function, n)

    """gives the number of steps in an exact permanent function"""

    if permanent_function == ryser
        n*2^n
    elseif permanent_function == naive
        factorial(big(n))
    elseif permanent_function == naive_tensor
        factorial(big(n))^2
    else
        error("not implemented")
    end
end


text:
using Genie.Router
using Reactive
using DataStructures

route("/") do

end

route("/test") do
    "Hello World!"
end

function service_create(str, n, image, internalport)
    p = port_base[str]
    s = String.(split(str, "_"))
    function instance(i)
        op = p + i
        if image == "mdpetters/apn"
            read(
                `docker service create -d --name=$(str)-$(i) --publish $(op):$(internalport) --network=$(s[1]) $(image) jupyter notebook --NotebookApp.token= --NotebookApp.password=`,
            )
        else
            read(
                `docker service create -d --name=$(str)-$(i) --publish $(op):$(internalport) --network=$(s[1]) $(image)`,
            )
        end
    end
    map(instance, 0:n-1)
end

IP = "127.0.0.1"

# Service virtualDMA: "swarm_vdma"
lab1 = "swarm_vdma"
n_clients1 = 1
vdma_ports = Stack{Int}()
image1 = "mdpetters/virtualdma:server"

lab2 = "swarm_htestbed"
testbed_ports = Stack{Int}()
n_clients2 = 2
image2 = "mdpetters/testbed:server"

lab3 = "swarm_apn"
apn_ports = Stack{Int}()
n_clients3 = 2
image3 = "mdpetters/apn"

lab4 = "swarm_invert"
invert_ports = Stack{Int}()
n_clients4 = 1
image4 = "mdpetters/inverttdma:server"

resolve_ports = Dict(lab1 => vdma_ports, lab2 => testbed_ports, lab3 => apn_ports, lab4 => invert_ports)
port_base = Dict{String,Int}(lab1 => 1000, lab2 => 1020, lab3 => 1040, lab4 => 1260)

map(i -> push!(resolve_ports[lab1], i), port_base[lab1]:port_base[lab1]+n_clients1-1)
map(i -> push!(resolve_ports[lab2], i), port_base[lab2]:port_base[lab2]+n_clients2-1)
map(i -> push!(resolve_ports[lab3], i), port_base[lab3]:port_base[lab3]+n_clients3-1)
map(i -> push!(resolve_ports[lab4], i), port_base[lab4]:port_base[lab4]+n_clients4-1)

run(`docker network create -d overlay swarm --attachable`)
service_create(lab1, n_clients1, image1, 1234)
service_create(lab2, n_clients2, image2, 1234)
service_create(lab3, n_clients3, image3, 8888)
service_create(lab4, n_clients4, image4, 1234)

sleep(60 * 15)

# Service
route("virtualTDMA") do
    if isempty(vdma_ports)
        "Sorry, all containers are checked out. This resource is currently unavailable. Please check back again later. If this issue persists, please contact mdpetter@ncsu.edu"
    else
        p1 = pop!(resolve_ports[lab1])
        println("Checking out port $(p1)")
        println("Available ports")
        println(resolve_ports[lab1])
        Genie.Renderer.redirect("http://$(IP):$(p1)/open?path=webapp.jl")
    end
end

route("hygroscopicityTestbed") do
    if isempty(testbed_ports)
        "Sorry, all containers are checked out. This resource is currently unavailable. Please check back again later. If this issue persists, please contact mdpetter@ncsu.edu"
    else
        p2 = pop!(resolve_ports[lab2])
        println("Checking out port $(p2)")
        println("Available ports")
        println(resolve_ports[lab2])
        Genie.Renderer.redirect("http://$(IP):$(p2)/open?path=webapp.jl")
    end
end

route("apn") do
    if isempty(apn_ports)
        "Sorry, all containers are checked out. This resource is currently unavailable. Please check back again later. If this issue persists, please contact mdpetter@ncsu.edu"
    else
        p3 = pop!(resolve_ports[lab3])
        println("Checking out port $(p3)")
        println("Available ports")
        println(resolve_ports[lab3])
        Genie.Renderer.redirect("http://$(IP):$(p3)/tree?")
    end
end

route("invertHTDMA") do
	if isempty(testbed_ports)
		"Sorry, all containers are checked out. This resource is currently unavailable. Please check back again later. If this issue persists, please contact mdpetter@ncsu.edu"
	else
		p4 = pop!(resolve_ports[lab4])
		println("Checking out port $(p4)")
		println("Available ports")
		println(resolve_ports[lab4])
		Genie.Renderer.redirect("http://$(IP):$(p4)/open?path=webapp.jl")
	end
end

include("monitor.jl")
containerActivity = Signal(Dict{String,Number}())
timer = fps(1.0 / (15.0 * 60.0))
mylog = map(_ -> containerManager(), timer)


text:
"""
    function discrete_gamma_rates(α::T, β::S,k::Int64, method::Symbol=:mean; sig::Int64=1)::Array{Float64} where {T<:Real, S<:Real}

This function calculates the rates of the discretized gamma distribution.
It returns a vector of either the mean or the median weights in each category,
such that each category is of equal proportion.

See: Yang, 1994, Maximum likelihood phylogenetic estimation from DNA sequences
with variable rates over sites: Approximate methods. (https://doi.org/10.1007/BF00160154)
"""
function discrete_gamma_rates(α::T, β::S,k::Int64, method::Symbol=:mean; sig::Int64=1)::Array{Float64} where {T<:Real, S<:Real}
    meanvals = Array{Float64, 1}(undef, k)
    factor::Float64 = α/β*k

    if method==:median
        meanvals .= median_boundaries(α, β, k)
        t = sum(meanvals)
        meanvals .*= (factor / t)
    else
        meanvals .= mean_boundaries(α, β, k)
        @inbounds @simd for i in 1:k-1
            meanvals[i] = gamma_inc(α+1,meanvals[i]*β, sig)[1]
        end
        meanvals[k] = 1.0
        @inbounds for i in k:-1:2
            meanvals[i] -= meanvals[i-1]
            meanvals[i] *= factor
        end
        meanvals[1] *= factor
    end
    meanvals
end

#### helper functions for boundaries of rate categories ####

function mean_boundaries(α::T, β::S, k::Int64)::Array{Float64} where {T<:Real, S<:Real}
    ch = Chisq((2*α)/(2*β))
    boundaries = Array{Float64, 1}(undef, k)
    for i in 1:k-1
        boundaries[i] = quantile(ch, i/k)/(2*β)
    end
    boundaries
end

function median_boundaries(α::T, β::S, k::Int64)::Array{Float64} where {T<:Real, S<:Real}
    ch = Chisq((2*α)/(2*β))
    boundaries = Array{Float64, 1}(undef, k)
    for i in 1:k
        boundaries[i] = quantile(ch, ((i-1)*2+1)/(2*k))
    end
    return boundaries
end


text:
using Curvelet
using Base.Test

N           = 512
sample_data = randn(N,N)
plan        = Curvelet.planCurveletTransform(N, Curvelet.C1)
tmp         = Curvelet.curveletTransform(sample_data, plan)
sample_hat  = Curvelet.inverseCurveletTransform(tmp, plan)
@test norm(sample_hat.-sample_hat) < 0.01


text:
using ConfParser
using Base.Test

conf = ConfParse("confs/config.http")
parse_conf!(conf)

@test retrieve(conf, "email")    == "juliarocks@socks.com"
@test retrieve(conf, "password") == "qwerty"
@test commit!(conf, "email", "newemail@test.com") == true

save!(conf, "confs/out.conf")


text:
type VineGraphChart <: JFinM_Charts.AbstractChart
    chartType::String
    cmd::String
    shape::String
    maxLayer::Int
    styles::Array{String, 1}
    colors::Array{String, 1}
    lens::Array{String, 1}
    weights::Array{String, 1}
end

## define default values
VineGraphChart() = VineGraphChart("VineGraphChart",
                                  "neato",
                                  "circle",
                                  2,
                                  ["bold", "dotted"],
                                  ["black", "black"],
                                  ["1.0", "1.8"],
                                  ["1", "1"]
                                  )

function GViz(vn::Copulas.Vine,
              chrt::VineGraphChart)
    nVars = size(vn.trees, 1)
    layers = Copulas.linkLayers(vn)
    dotCode = ""
    dotCode = string(dotCode, "graph G {\n")
    dotCode = string(dotCode, "node [shape=$(chrt.shape)];\n")
    for ii=1:nVars
        for jj=ii+1:nVars
            currLayer = layers[ii, jj]
            if currLayer <= chrt.maxLayer
                tmpCode = "$ii -- $jj [style=$(chrt.styles[currLayer]); color=$(chrt.colors[currLayer]); len=$(chrt.lens[currLayer]); weight=$(chrt.weights[currLayer])];\n"
                    dotCode = string(dotCode, tmpCode)
            end
        end
    end
    dotCode = string(dotCode, "}")
    return GViz("VineGraphChart", "neato", dotCode)
end


text:
using Documenter, Hospitalplanning

makedocs(;
    modules=[Hospitalplanning],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/HBreddam/Hospitalplanning.jl/blob/{commit}{path}#L{line}",
    sitename="Hospitalplanning.jl",
    authors="Henrik Bøgedal Breddam",
    assets=String[],
)

deploydocs(;
    repo="github.com/HBreddam/Hospitalplanning.jl",
)


text:
import DataStructures: OrderedDict

"""
```
HetDSGELag{T} <: AbstractModel{T}
```

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"an_schorfheide\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.

"""
mutable struct HetDSGELag{T} <: AbstractModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values

    # Temporary to get it to work. Need to
    # figure out a more flexible way to define
    # "grids" that are not necessarily quadrature
    # grids within the model
    grids::OrderedDict{Symbol,Union{Grid, Array}}

    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    state_variables::Vector{Symbol}                        # Vector of symbols of the state variables
    jump_variables::Vector{Symbol}                         # Vector of symbols of the jump variables
    normalized_model_states::Vector{Symbol}                # All of the distributional model
                                                           # state variables that need to be normalized
    endogenous_states_unnormalized::OrderedDict{Symbol,UnitRange} # Vector of unnormalized
                                                           # ranges of indices
    endogenous_states::OrderedDict{Symbol,UnitRange}       # Vector of ranges corresponding
                                                           # to normalized (post Klein solution) indices
    exogenous_shocks::OrderedDict{Symbol,Int}              #
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,UnitRange}  #
    observables::OrderedDict{Symbol,Int}                   #

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
end

description(m::HetDSGELag) = "HetDSGELag, $(m.subspec)"

"""
`init_model_indices!(m::HetDSGELag)`

Arguments:
`m:: HetDSGELag`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::HetDSGELag)
    # Endogenous states
    endogenous_states = collect([:l′_t1, :μ′_t1, :k′_t, :R′_t1, :i′_t1,
                                 :π′_t1,:π′_t2,:π′_t3,:y′_t1,:y′_t2,:y′_t3, :y′_t4,
                                 :z′_t1,:z′_t2,:z′_t3,:w′_t1,:I′_t1,:B′,:G′,:z′_t,:MU′,:LAMW′,
                                 :LAMF′,:MON′,:l′_t,:μ′_t,:R′_t,:i′_t,:t′_t,:w′_t,:L′_t,:π′_t,
                                 :wageinflation′_t,:mu′_t,:y′_t,:I′_t,:mc′_t,:Q′_t,:capreturn′_t])
                                 #=:l_t1,:μ_t1,:k_t,:R_t1,:i_t1,:π_t1,:π_t2,:π_t3,:y_t1,:y_t2,
                                 :y_t3,:y_t4,:z_t1,:z_t2,:z_t3,:w_t1,:I_t1,:B,:G,:Z,:MU,:LAMW,
                                 :LAMF,:MON,:ELL,:μ_t,:R_t,:i_t,
                                 :t_t,:w_t,:L_t,:π_t,:wageinflation_t,:mu_t,:y_t,:I_t,:mc_t,:Q_t,:capreturn_t]) =#

    # endogenous scalar states: k_t, R_t1, i_t1, π_t1, π_t2, π_t3, y_t1, y_t2, y_t3, y_t4, z_1, z_2, z_3, w_t1, x_t1

    # Exogenous shocks
    exogenous_shocks = collect([:b_sh,:g_sh,:z_sh,:μ_sh,:λ_w_sh, :λ_f_sh,:MON_sh])

    # Equilibrium conditions
    equilibrium_conditions = collect([:eq_euler,:eq_kolmogorov_fwd,:eq_lag_ell,:eq_lag_wealth,
                                      :eq_market_clearing,:eq_lambda, :eq_transfers,
                                      :eq_investment,:eq_tobin_q,:eq_capital_accumulation,
                                      :eq_wage_phillips,:eq_price_phillips,:eq_marginal_cost,
                                      :eq_gdp,:eq_optimal_kl,:eq_taylor, :eq_fisher,
                                      :eq_nominal_wage_inflation,:LR,:LI,:LPI,:L2PI,:L3PI,:LY,:L2Y,
                                      :L3Y,:L4Y,:LZ,:L2Z,:L3Z,:LW,:LX,:F33,:F34,:F35,:F36,:F37,
                                      :F38, :F39])

    # Observables
    observables = keys(m.observable_mappings)

    ########################################################################################
    # Setting indices of endogenous_states and equilibrium conditions manually for now

    setup_indices!(m)
    endo = m.endogenous_states_unnormalized
    eqcond = equilibrium_conditions

    ########################################################################################

    m.normalized_model_states = [:μ′_t1, :l′_t1]

    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
end

function HetDSGELag(subspec::String="ss0";
                   custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                   testing = false)

    # Model-specific specifications
    spec               = "het_dsge_lag"
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = HetDSGELag{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(),
            # grids and keys
            OrderedDict{Symbol,Union{Grid, Array}}(), OrderedDict{Symbol,Int}(),

            # normalized_model_states, state_inds, jump_inds
            Vector{Symbol}(), Vector{Symbol}(), Vector{Symbol}(),

            # model indices
            # endogenous states unnormalized, endogenous states normalized
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,UnitRange}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,UnitRange}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}())

    m <= Setting(:nx, 300)
    m <= Setting(:ns, 2)

    default_settings!(m)

    # # Set observable transformations
    init_observable_mappings!(m)

    # Initialize model indices
    init_model_indices!(m)

    # Set settings
    model_settings!(m)
    # default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Initialize parameters
    init_parameters!(m)

    # Initialize aggregate steady state parameters (necessary for grid construction)
    aggregate_steadystate!(m)

    # Initialize grids
    init_grids!(m)

    # # Solve for the steady state
    steadystate!(m)

    # # So that the indices of m.endogenous_states reflect the normalization
    normalize_model_state_indices!(m)

    init_subspec!(m)

    return m
end

"""
```
init_parameters!(m::HetDSGELag)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::HetDSGELag)

    # Initialize parameters
    m <= parameter(:r, 0.01, fixed = true,
                   description = "r: Steady-state real interest rate.", tex_label = "r")
    m <= parameter(:α, 0.3, fixed = true,
                   description = "α: Capital share", tex_label = "\\alpha")
    m <= parameter(:H, 1.0, fixed = true,
                   description = "H: Aggregate hours worked", tex_label = "H")
    m <= parameter(:δ, 0.03, fixed = true,
                   description = "δ: Depreciation rate", tex_label = "\\delta")
    m <= parameter(:μ_sp, 0.0, fixed = true, description = "μ_sp: The trend in the skill process",
                   tex_label = "\\mu_{sp}")
    m <= parameter(:ρ_sp, 0.7, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(), BetaAlt(0.5, 0.2), fixed=true,
                   description="ρ_sp: AR(1) coefficient in the skill process.",
                   tex_label="\\rho_{sp}")
    m <= parameter(:σ_sp, 0.01, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed=true,
                   description="σ_sp: The standard deviation of the skill process.",
                   tex_label="\\sigma_{sp}")
    m <= parameter(:γ, 0.0, description = "γ: TFP growth")
    #GoverY = 0.01
    m <= parameter(:g, 1/(1-0.01), description = "g: Steady-state government spending/gdp")
    m <= parameter(:η, 0.1, description = "η: Borrowing constraint (normalized by TFP)")

    m <= parameter(:ρ_B, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = "ρ_B: Persistence of discount factor shock", tex_label = "\\rho_B")
    m <= parameter(:ρ_G, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = "ρ_G: Persistence of govt spending shock", tex_label = "\\rho_G")
    m <= parameter(:ρ_z, 0.5,(1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = "ρ_z: Persistence of tfp growth shock", tex_label = "\\rho_z")
    m <= parameter(:ρ_μ, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = "ρ_μ: Persistence of investment shock", tex_label = "\\rho_\\mu")
    m <= parameter(:ρ_lamw, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = "ρ_λ_w: Persistence of wage mkup shock", tex_label = "\\rho_{\\lambda_w}")
    m <= parameter(:ρ_lamf, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = " # persistence of price mkup shock", tex_label = "\\rho_{\\lambda_f}")
    m <= parameter(:ρ_mon, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                     fixed=true, description = " # persistence of mon policy shock", tex_label = "\\rho_{mon}")

    m <= parameter(:σ_b, 0.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=true, description = "# standard deviation of discount factor shock", tex_label = "\\sigma_B")
    m <= parameter(:σ_g, 0.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=true, description = "# standard deviation of govt spending shock", tex_label = "\\sigma_G")
    m <= parameter(:σ_z, 0.5,(1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=false, description = "# standard deviation of tfp growth shock", tex_label = "\\sigma_z")
    m <= parameter(:σ_μ, 0.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=true, description = " # standard deviation of investment shock", tex_label = "\\sigma_\\mu")
    m <= parameter(:σ_λ_w, 0.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=true, description = "# standard deviation of wage mkup shock", tex_label = "\\sigma_{\\lambda_w}")
    m <= parameter(:σ_λ_f, 0.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=true, description = " # standard deviation of price mkup shock", tex_label = "\\sigma_{\\lambda_f}")
    m <= parameter(:σ_MON, 0.5, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.1),
                   fixed=true, description = " # standard deviation of mon policy shock", tex_label = "\\sigma_{mon}")

 m <= parameter(:spp, 4., description = "# second derivative of investment adjustment cost")
 m <= parameter(:lamw, 1.5, description = "# wage markup")
 m <= parameter(:ϕh , 2., description = "# inverse frisch elasticity")
 m <= parameter(:Φw , 10., description = "# rotemberg cost for wages")
 m <= parameter(:lamf, 1.5 , description = "# price markup")
 m <= parameter(:Φp , 1., description = "# rotemberg cost for prices")
 m <= parameter(:ρR , 0.75, description = "# persistence in taylor rule")
 m <= parameter(:ψπ , 1.5, description = "# weight on inflation in taylor rule")
 m <= parameter(:ψy , 0.5, description = "# weight on output growth in taylor rule")

    # Setting steady-state parameters
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)

    # Steady state grids for functional/distributional variables and market-clearing discount rate
    m <= SteadyStateParameterGrid(:lstar, fill(NaN, nx*ns), description = "Steady-state expected discounted
                                  marginal utility of consumption", tex_label = "l_*")
    m <= SteadyStateParameterGrid(:cstar, fill(NaN, nx*ns), description = "Steady-state consumption",
                                  tex_label = "c_*")
    m <= SteadyStateParameterGrid(:μstar, fill(NaN, nx*ns), description = "Steady-state cross-sectional
                                  density of cash on hand", tex_label = "\\mu_*")
    m <= SteadyStateParameter(:βstar, NaN, description = "Steady-state discount factor",
                              tex_label = "\\beta_*")
end

"""
```
aggregate_steadystate!(m::HetDSGELag)
```
Steady state of aggregate scalar variables (Break out the "analytic" steady-state solution here instead of in steadystate!() since some of these parameters are required to construct the grids that are used in the functional/distributional steady state variable calculation (where init_grids! is called before steadystate! in model construction).
"""
function aggregate_steadystate!(m::HetDSGELag)
    m <= SteadyStateParameter(:Rkstar, m[:r] + m[:δ], description = "Rental rate on capital", tex_label = "Rk_*")
    m <= SteadyStateParameter(:ωstar, (m[:α]^(m[:α]/(1-m[:α])))*(1-m[:α])*m[:Rkstar]^(-m[:α]/(1-m[:α])), description = "Real wage", tex_label = "\\omega_*")
    m <= SteadyStateParameter(:klstar, (m[:α]/(1-m[:α]))*(m[:ωstar]/m[:Rkstar])*exp(m[:γ]), description = "Capital/Labor ratio", tex_label = "kl_*")
    m <= SteadyStateParameter(:kstar, m[:klstar]*m[:H], description = "Capital", tex_label = "k_*")
    m <= SteadyStateParameter(:xstar, (1-(1-m[:δ])*exp(-m[:γ]))*m[:kstar], description = "Investment", tex_label = "x_*")
    m <= SteadyStateParameter(:ystar, (exp(-m[:α]*m[:γ])*(m[:kstar]^m[:α])*m[:H]^(1-m[:α])), description = "GDP", tex_label = "y_*")
    m <= SteadyStateParameter(:Tstar, m[:Rkstar]*m[:kstar]*exp(-m[:γ]) - m[:xstar] - (1-(1/m[:g]))*m[:ystar], description = "Net transfer to households", tex_label = "T_*")

    return m
end

"""
```
init_grids!(m::HetDSGELag)
```
"""
function init_grids!(m::HetDSGELag)

    nx      = get_setting(m, :nx)
    ns      = get_setting(m, :ns)
    λ       = get_setting(m, :λ)

    grids = OrderedDict()

    # Skill grid
    #lsgrid, sprob, sscale = tauchen86(m[:μ_sp].value, m[:ρ_sp].value, m[:σ_sp].value, ns, λ)
    sprob = [[0.9 0.1];[0.1 0.9]]
    sgrid = [0.8;1.2]
    (λs, vs) = eigen(Matrix{Float64}(sprob'))
    order_λs = sortperm(λs, rev = true)
    vs = vs[:,order_λs]
    ss_skill_distr = vs[:,1]/sum(vs[:,1])
    means = ss_skill_distr'*sgrid
    meanz = (get_setting(m, :zhi)+get_setting(m, :zlo))/2.
    sgrid = sgrid/(meanz*means) # so that skills integrate to 1
    sscale = sgrid[2] - sgrid[1]
    swts = (sscale/ns)*ones(ns)
    #sgrid = exp.(lsgrid)
    grids[:sgrid] = Grid(sgrid, swts, sscale)

 # Calculate endogenous grid bounds (bounds are analytic functions of model parameters,
    # which ensure there is no mass in the distributions across x where there shouldn't be)
    # Calculate the lowest possible cash on hand in steady state
    zlo = get_setting(m, :zlo)
    smin = minimum(sgrid)*zlo #exp(m[:μ_sp]/(1-m[:ρ_sp]) - get_setting(m, :λ)*sqrt(m[:σ_sp]^2/(1-m[:ρ_sp])^2))*get_setting(m, :zlo)
    m <= Setting(:xlo, m[:ωstar]*smin*m[:H] - (1+m[:r])*m[:η]*exp(-m[:γ]) + m[:Tstar] + sgrid[1]*m[:ωstar]*m[:H]*0.05, "Lower bound on cash on hand")
    m <= Setting(:xhi, max(get_setting(m, :xlo)*2, get_setting(m, :xlo)+5.), "Upper Bound on cash on hand")
    m <= Setting(:xscale, get_setting(m, :xhi) - get_setting(m, :xlo), "Size of the xgrid")
    xlo     = get_setting(m, :xlo)
    xhi     = get_setting(m, :xhi)
    xscale  = get_setting(m, :xscale)

    # Cash on hand grid
    grids[:xgrid] = Grid(uniform_quadrature(xscale), xlo, xhi, nx, scale = xscale)

    # Markov transition matrix for skill
    grids[:fgrid] = sprob./swts

    # Total grid vectorized across both dimensions
    grids[:sgrid_total] = kron(sgrid, ones(nx))
    grids[:xgrid_total] = kron(ones(ns), grids[:xgrid].points)
    grids[:weights_total] = kron(swts, grids[:xgrid].weights)

    m.grids = grids
end

function model_settings!(m::HetDSGELag)
    default_settings!(m)

    # Defaults
    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), DSGE_DATE_FORMAT)
    m <= Setting(:data_vintage, vint, true, "vint", "Data vintage")

    saveroot = normpath(joinpath(dirname(@__FILE__), "../../../","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "../../../","save","input_data"))

    m <= Setting(:saveroot, saveroot, "Root of data directory structure")
    m <= Setting(:dataroot, datapath, "Input data directory path")

    # Solution method
    m <= Setting(:solution_method, :klein)

    m <= Setting(:krylov_reduce, false)

    # Likelihood method
    m <= Setting(:use_chand_recursion, true)

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0,
                 "Number of anticipated policy shocks")
    # m <= Setting(:n_anticipated_shocks_padding, 20,
                 # "Padding for anticipated policy shocks")

    # Whether one wishes to re-compute βstar
    m <= Setting(:use_last_βstar, false, "Flag to avoid recomputing steady-state.")

    # Number of states and jumps
    m <= Setting(:normalize_distr_variables, true, "Whether or not to perform the
                 normalization of the distributional states in the Klein solution step")
    m <= Setting(:n_predetermined_variables, 0, "Number of predetermined variables after
                 removing the densities. Calculated with the Jacobian normalization.
                 This set to 0 is just the default setting, since it will always be
                 overwritten once the Jacobian is calculated.")

    endo = m.endogenous_states_unnormalized
    state_indices = [endo[:l′_t1]; endo[:μ′_t1];  endo[:k′_t]; endo[:R′_t1];endo[:i′_t1];endo[:π′_t1];
                     endo[:π′_t2];endo[:π′_t3];endo[:y′_t1];endo[:y′_t2];endo[:y′_t3];endo[:y′_t4];
                     endo[:z′_t1];endo[:z′_t2];endo[:z′_t3];endo[:w′_t1];endo[:I′_t1];endo[:B′];
                     endo[:G′];endo[:z′_t];endo[:MU′];endo[:LAMW′]; endo[:LAMF′];endo[:MON′]]
                     #=endo[:l_t1];
                     endo[:μ_t1];endo[:k_t];endo[:R_t1];endo[:i_t1]; endo[:π_t1];endo[:π_t2];
                     endo[:π_t3]; endo[:y_t1];endo[:y_t2];endo[:y_t3]; endo[:y_t4];endo[:z_t1];
                     endo[:z_t2]; endo[:z_t3];endo[:w_t1];endo[:I_t1]; endo[:B];endo[:G];endo[:Z];
                     endo[:MU]; endo[:LAMW];endo[:LAMF];endo[:MON]]=#

    jump_indices = [endo[:l′_t]; endo[:μ′_t];endo[:R′_t];endo[:i′_t];endo[:t′_t];endo[:w′_t];
                    endo[:L′_t];endo[:π′_t];endo[:wageinflation′_t];endo[:mu′_t];endo[:y′_t];
                    endo[:I′_t];endo[:mc′_t]; endo[:Q′_t];endo[:capreturn′_t]]
                    #=endo[:ELL];
                    endo[:μ_t];endo[:R_t];endo[:i_t];endo[:t_t];endo[:w_t]; endo[:L_t];endo[:π_t];
                    endo[:wageinflation_t];endo[:mu_t];endo[:y_t];endo[:I_t];endo[:mc_t];endo[:Q_t];
                    endo[:capreturn_t]]=#

    m <= Setting(:state_indices, state_indices, "Which indices of m.endogenous_states correspond to
                 backward looking state variables")
    m <= Setting(:jump_indices, jump_indices, "Which indices of m.endogenous_states correspond to jump
                 variables")
    #m.state_variables = m.endogenous_states.keys[get_setting(m, :state_indices)]
    #m.jump_variables = m.endogenous_states.keys[get_setting(m, :jump_indices)]

    # Mollifier setting parameters
    m <= Setting(:In, 0.443993816237631, "Normalizing constant for the mollifier")
    m <= Setting(:zlo, 1.0, "Lower bound on second income shock to mollify actual income")
    m <= Setting(:zhi, 5.0, "Upper bound on second income shock to mollify actual income")

    # s: Skill Distribution/ "Units of effective labor" Grid Setup
    m <= Setting(:ns, 2, "Skill distribution grid points")
    m <= Setting(:λ, 2.0, "The λ parameter in the Tauchen distribution calculation")

    # x: Cash on Hand Grid Setup
    m <= Setting(:nx, 300, "Cash on hand distribution grid points")

    m <= Setting(:trunc_distr, true)
    m <= Setting(:rescale_weights, true)
    m <= Setting(:mindens, 1e-8)

    # Function-valued variables include distributional variables
    m <= Setting(:n_function_valued_backward_looking_states, 2, "Number of function-valued
                 backward looking state variables")
    m <= Setting(:n_backward_looking_distributional_vars, 2, "Number of state variables that are
                 distributional variables.")
    m <= Setting(:n_function_valued_jumps, 2, "Number of function-valued jump variables")
    m <= Setting(:n_jump_distributional_vars, 2, "Number of jump variables that are distributional
                 variables.")

    # Note, these settings assume normalization.
    # The n degrees of freedom removed depends on the distributions/dimensions
    # of heterogeneity that we have discretized over, in this case,
    # cash on hand and the skill distribution. In general the rule of
    # thumb is, remove one degree of freedom for the first endogenous distribution (cash on
    # hand), then one additional degree of freedom for each exogenous distribution (skill
    # distribution). Multiple endogenous distributions only permit removing a single degree
    # of freedom since it is then non-trivial to obtain the marginal distributions.
    m <= Setting(:n_degrees_of_freedom_removed, 1, "Number of degrees of freedom from the
                 distributional variables to remove.")
    n_dof_removed = get_setting(m, :n_degrees_of_freedom_removed)

    ####################################################
    # Calculating the number of backward-looking states
    ####################################################
    n_backward_looking_distr_vars = get_setting(m, :n_backward_looking_distributional_vars)
    m <= Setting(:backward_looking_states_normalization_factor,
                 n_dof_removed*n_backward_looking_distr_vars, "The number of dimensions removed from the
                 backward looking state variables for the normalization.")

    n = get_setting(m, :nx)*get_setting(m, :ns)
    n_backward_looking_vars = length(get_setting(m, :state_indices))
    n_backward_looking_function_valued_vars = get_setting(m, :n_function_valued_backward_looking_states)
    n_backward_looking_scalar_vars = n_backward_looking_vars - n_backward_looking_function_valued_vars

    m <= Setting(:n_backward_looking_states, n*n_backward_looking_distr_vars +
                 n_backward_looking_scalar_vars - get_setting(m, :backward_looking_states_normalization_factor),
                 "Number of state variables, in the true sense (fully
                  backward looking) accounting for the discretization across the grid
                  of function-valued variables and the normalization of
                  distributional variables.")

    ##################################
    # Calculating the number of jumps
    ##################################
    n_jump_distr_vars = get_setting(m, :n_jump_distributional_vars)
    m <= Setting(:jumps_normalization_factor,
                 n_dof_removed*n_jump_distr_vars, "The number of dimensions removed from the
                 jump variables for the normalization.")

    n_jump_vars = length(get_setting(m, :jump_indices))
    n_jump_function_valued_vars = get_setting(m, :n_function_valued_jumps)
    n_jump_scalar_vars = n_jump_vars - n_jump_function_valued_vars

    m <= Setting(:n_jumps, n*n_jump_function_valued_vars +
                 n_jump_scalar_vars - get_setting(m, :jumps_normalization_factor),
                 "Number of jump variables (forward looking) accounting for
                  the discretization across the grid of function-valued variables and the
                  normalization of distributional variables.")

    m <= Setting(:n_model_states, get_setting(m, :n_backward_looking_states) + get_setting(m, :n_jumps),
                 "Number of 'states' in the state space model. Because backward and forward
                 looking variables need to be explicitly tracked for the Klein solution
                 method, we have n_states and n_jumps")
end

function setup_indices!(m::HetDSGELag)
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)
    endo = m.endogenous_states_unnormalized
    eqconds = m.equilibrium_conditions
    nxns = nx*ns

    # Endogenous function-valued states
    endo[:l′_t1]    = 1:nxns            # lagged ell function
    endo[:μ′_t1]    = nxns+1:2*nxns     # lagged wealth distribution
    #endogenous scalar-valued states
    endo[:k′_t]  = 2*nxns+1:2*nxns+1             # capital –dont get confused with steadystate object K
    endo[:R′_t1] = 2*nxns+2:2*nxns+2              # lagged real interest rate
    endo[:i′_t1] = 2*nxns+3:2*nxns+3             # lagged nominal interest rate
    endo[:π′_t1] = 2*nxns+4:2*nxns+4             # lagged inflation
    endo[:π′_t2] = 2*nxns+5:2*nxns+5             # double-lagged inflation
    endo[:π′_t3] = 2*nxns+6:2*nxns+6             # triple-lagged inflation
    endo[:y′_t1] = 2*nxns+7:2*nxns+7             # lagged gdp
    endo[:y′_t2] = 2*nxns+8:2*nxns+8             # double-lagged gdp
    endo[:y′_t3] = 2*nxns+9:2*nxns+9             # triple-lagged gdp
    endo[:y′_t4] = 2*nxns+10:2*nxns+10           # quad-lagged gdp
    endo[:z′_t1] = 2*nxns+11:2*nxns+11           # lag tfp growth
    endo[:z′_t2] = 2*nxns+12:2*nxns+12           # double-lag tfp growth
    endo[:z′_t3] = 2*nxns+13:2*nxns+13           # triple-lag tfp growth
    endo[:w′_t1] = 2*nxns+14:2*nxns+14           # lag real wages
    endo[:I′_t1] = 2*nxns+15:2*nxns+15           # lag investment–don't get this confused with i, the nominal interest rate
    # exogenous scalar-valued states:
    endo[:B′]    = 2*nxns+16:2*nxns+16        # discount factor shock
    endo[:G′]    = 2*nxns+17:2*nxns+17        # govt spending
    endo[:z′_t]    = 2*nxns+18:2*nxns+18        # tfp growth
    endo[:MU′]   = 2*nxns+19:2*nxns+19        # investment shock
    endo[:LAMW′] = 2*nxns+20:2*nxns+20        # wage markup
    endo[:LAMF′] = 2*nxns+21:2*nxns+21        # price markup
    endo[:MON′]  = 2*nxns+22:2*nxns+22        # monetary policy shock

 # function-valued jumps
    endo[:l′_t]  = 2*nxns+23:3*nxns+22 # ell function
    endo[:μ′_t]  = 3*nxns+23:4*nxns+22 # wealth distribution
    #scalar-valued jumps
    endo[:R′_t]  = 4*nxns+23:4*nxns+23        # real interest rate
    endo[:i′_t]  = 4*nxns+24:4*nxns+24        # nominal interest rate
    endo[:t′_t]  = 4*nxns+25:4*nxns+25        # transfers + dividends
    endo[:w′_t]  = 4*nxns+26:4*nxns+26        # real wage
    endo[:L′_t]  = 4*nxns+27:4*nxns+27        # hours worked
    endo[:π′_t]  = 4*nxns+28:4*nxns+28        # inflation
    endo[:wageinflation′_t] = 4*nxns+29:4*nxns+29        # nominal wage inflation
    endo[:mu′_t] =  4*nxns+30:4*nxns+30        # average marginal utility
    endo[:y′_t]  = 4*nxns+31:4*nxns+31       # gdp
    endo[:I′_t]  = 4*nxns+32:4*nxns+32        # investment
    endo[:mc′_t] = 4*nxns+33:4*nxns+33        # marginal cost - this is ζ in HetDSGELagₖd.pdf
    endo[:Q′_t]  = 4*nxns+34:4*nxns+34        # Tobin's qfunction
    endo[:capreturn′_t] = 4*nxns+35:4*nxns+35        # return on capital

    nvars = 4*nxns+35
    nscalars = 35 # num eqs which output scalars
    nyscalars = 13 # num scalar jumps
    nxscalars = nscalars - nyscalars # num scalar states

#=endo[:l_t1] = endo[:l′_t1] .+ nvars #LELL
endo[:μ_t1] = endo[:μ′_t1] .+ nvars #LM
endo[:k_t] = endo[:k′_t] .+ nvars #KK
endo[:R_t1] = endo[:R′_t1] .+ nvars #LRR
endo[:i_t1] = endo[:i′_t1] .+ nvars #LII
endo[:π_t1] = endo[:π′_t1]  .+ nvars #LPI
endo[:π_t2] = endo[:π′_t2]  .+ nvars #L2PI
endo[:π_t3] = endo[:π′_t3]  .+ nvars #L3PI
endo[:y_t1] = endo[:y′_t1] .+ nvars #LY
endo[:y_t2] = endo[:y′_t2] .+ nvars #L2Y
endo[:y_t3] = endo[:y′_t3] .+ nvars #l3Y
endo[:y_t4] = endo[:y′_t4] .+ nvars #L4Y
endo[:z_t1] = endo[:z′_t1] .+ nvars #LZ
endo[:z_t2] = endo[:z′_t2] .+ nvars #L2Z
endo[:z_t3] = endo[:z′_t3] .+ nvars #L3Z
endo[:w_t1] = endo[:w′_t1] .+ nvars #LW
endo[:I_t1] = endo[:I′_t1] .+ nvars #LX
endo[:B] = endo[:BP] .+ nvars #B
endo[:G] = endo[:GP] .+ nvars #G
endo[:Z] = endo[:ZP] .+ nvars #Z
endo[:MU] = endo[:MUP] .+ nvars #MU
endo[:LAMW] = endo[:LAMWP] .+ nvars #LAMW
endo[:LAMF] = endo[:LAMFP] .+ nvars #LAMF
endo[:MON] = endo[:MONP] .+ nvars #MON
endo[:ELL] = endo[:l′_t] .+ nvars #ELL
endo[:μ_t] = endo[:μ′_t] .+ nvars #M
endo[:R_t] = endo[:R′_t] .+ nvars #RR
endo[:i_t] = endo[:i′_t] .+ nvars #II
endo[:t_t] = endo[:t′_t]  .+ nvars #TT
endo[:w_t] = endo[:w′_t] .+ nvars #W
endo[:L_t] = endo[:L′_t] .+ nvars #HH
endo[:π_t] = endo[:π′_t] .+ nvars #PI
endo[:wageinflation_t] = endo[:wageinflation′_t] .+ nvars #PIW
endo[:mu_t] = endo[:mu′_t] .+ nvars #LAM
endo[:y_t] = endo[:y′_t] .+ nvars #Y
endo[:I_t] = endo[:I′_t]  .+ nvars #X
endo[:mc_t] = endo[:mc′_t]  .+ nvars #MC
endo[:Q_t] = endo[:Q′_t]  .+ nvars #Q
endo[:capreturn_t] = endo[:capreturn′_t] .+ nvars #RK =#

# create objects needed for solve.jl
# we will order function blocks as follows:
# 1. all function blocks which output a function (first real eqs, then lags)
# 2. all function blocks which map functions to scalars
# 3. all scalar blocks involving endogenous vbls (first real eqs, then lags)
# 4. shock processes
funops = 1:4 # which operators output a function

    # function blocks which output a function
    eqconds[:eq_euler]              = 1:nxns
    eqconds[:eq_kolmogorov_fwd]     = nx*ns+1:2*nx*ns
    eqconds[:eq_lag_ell]            = 2*nxns+1:3*nxns
    eqconds[:eq_lag_wealth]         = 3*nxns+1:4*nxns
    # function blocks which map functions to scalars
    eqconds[:eq_market_clearing]    = 4*nxns+1:4*nxns+1
    eqconds[:eq_lambda]             = 4*nxns+2:4*nxns+2
    #scalar blocks involving endogenous variables
    eqconds[:eq_transfers]            = 4*nxns+3:4*nxns+3 # transfers
    eqconds[:eq_investment]          = 4*nxns+4:4*nxns+4 # investment
    eqconds[:eq_tobin_q]              = 4*nxns+5:4*nxns+5 # tobin's q
    eqconds[:eq_capital_accumulation] = 4*nxns+6:4*nxns+6 # capital accumulation
    eqconds[:eq_wage_phillips] = 4*nxns+7:4*nxns+7 # wage phillips curve
    eqconds[:eq_price_phillips] = 4*nxns+8:4*nxns+8 # price phillips curve
    eqconds[:eq_marginal_cost]  = 4*nxns+9:4*nxns+9 # marginal cost
    eqconds[:eq_gdp]  = 4*nxns+10:4*nxns+10 # gdp
    eqconds[:eq_optimal_kl] = 4*nxns+11:4*nxns+11 # optimal K/L ratio
    eqconds[:eq_taylor] = 4*nxns+12:4*nxns+12 # taylor rule
    eqconds[:eq_fisher] = 4*nxns+13:4*nxns+13 # fisher eqn
    eqconds[:eq_nominal_wage_inflation] = 4*nxns+14:4*nxns+14 # nominal wage inflation
    # lagged variables
    eqconds[:LR] = 4*nxns+15:4*nxns+15 # LR
    eqconds[:LI] = 4*nxns+16:4*nxns+16 # LI
    eqconds[:LPI] = 4*nxns+17:4*nxns+17 # LPI
    eqconds[:L2PI] = 4*nxns+18:4*nxns+18 # L2PI
    eqconds[:L3PI] = 4*nxns+19:4*nxns+19 # L3PI
    eqconds[:LY] = 4*nxns+20:4*nxns+20 # LY
    eqconds[:L2Y] = 4*nxns+21:4*nxns+21 # L2Y
    eqconds[:L3Y] = 4*nxns+22:4*nxns+22 # L3Y
    eqconds[:L4Y] = 4*nxns+23:4*nxns+23 # L4Y
    eqconds[:LZ]  = 4*nxns+24:4*nxns+24 # LZ
    eqconds[:L2Z] = 4*nxns+25:4*nxns+25 # L2Z
    eqconds[:L3Z] = 4*nxns+26:4*nxns+26 # L3Z
    eqconds[:LW]  = 4*nxns+27:4*nxns+27 # LW
    eqconds[:LX] = 4*nxns+28:4*nxns+28 # LX
    # shocks
    eqconds[:F33] = 4*nxns+29:4*nxns+29 # discount factor B
    eqconds[:F34] = 4*nxns+30:4*nxns+30 # govt spending G
    eqconds[:F35] = 4*nxns+31:4*nxns+31 # tfp growth Z
    eqconds[:F36] = 4*nxns+32:4*nxns+32 # investment MU
    eqconds[:F37] = 4*nxns+33:4*nxns+33 # wage mkup LAMW
    eqconds[:F38] = 4*nxns+34:4*nxns+34 # price mkup LAMF
    eqconds[:F39] = 4*nxns+35:4*nxns+35 # monetary policy MON

# Total grid x*s
m <= Setting(:n, get_setting(m, :nx) * get_setting(m, :ns), "Total grid size, multiplying
                     across grid dimensions.")

m <= Setting(:nvars, 4*get_setting(m, :n) +35, "num variables")
m <= Setting(:nscalars, 35, " # num eqs which output scalars")
m <= Setting(:nyscalars, 13, "num scalar jumps")
m <= Setting(:nxscalars, 35 - 13, "num scalar states")

m.endogenous_states = deepcopy(endo)

end