text:
export eval_f_on_dof_pts, get_bdedge_integration_pts,
    gauss_weights, compute_fem_boundary_mass_matrix1, compute_fem_boundary_mass_term1,
    eval_scalar_on_boundary_edge, eval_strain_on_boundary_edge, eval_normal_and_shear_stress_on_boundary_edge

"""
    eval_f_on_gauss_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
"""
function eval_f_on_gauss_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
    xy = gauss_nodes(mesh)
    if tensor_input
        return f(constant(xy[:,1]), constant(xy[:,2]))
    end
    f.(xy[:,1], xy[:,2])
end

"""
    eval_f_on_fem_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
"""
function eval_f_on_fem_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
    xy = mesh.nodes
    if tensor_input
        return f(constant(xy[:,1]), constant(xy[:,2]))
    end
    f.(xy[:,1], xy[:,2])
end

"""
    eval_f_on_fvm_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
"""
function eval_f_on_fvm_pts(f::Function, mesh::Mesh; tensor_input::Bool = false)
    xy = fvm_nodes(mesh)
    if tensor_input
        return f(constant(xy[:,1]), constant(xy[:,2]))
    end
    f.(xy[:,1], xy[:,2])
end

"""
    eval_f_on_dof_pts(f::Function, mesh::Mesh)

Evaluates `f` on the DOF points.

- For P1 element, the DOF points are FEM points and therefore `eval_f_on_dof_pts` is equivalent to `eval_on_on_fem_pts`.
- For P2 element, the DOF points are FEM points plus the middle point for each edge.

Returns a vector of length `mesh.ndof`.
"""
function eval_f_on_dof_pts(f::Function, mesh::Mesh)
    if size(mesh.conn, 2)==3
        return eval_f_on_fem_pts(f, mesh)
    end
    xy = mesh.nodes
    xy2 = zeros(mesh.nedge, 2)
    for i = 1:mesh.nedge
        xy2[i,:] = (mesh.nodes[mesh.edges[i,1], :] + mesh.nodes[mesh.edges[i,2], :])/2
    end
    xy = [xy;xy2]
    f.(xy[:,1], xy[:,2])
end

"""
    compute_fem_source_term1(f::PyObject, mesh::Mesh)
"""
function compute_fem_source_term1(f::PyObject, mesh::Mesh)
    @assert length(f)==get_ngauss(mesh)
    fem_source_scalar_ = load_op_and_grad(libmfem,"fem_source_scalar")
    f = convert_to_tensor(Any[f], [Float64]); f = f[1]
    out = fem_source_scalar_(f)
    n = mesh.ndof
    set_shape(out, (n,))
end

"""
    compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh)
"""
function compute_fem_source_term1(f::Array{Float64,1}, mesh::Mesh)
    @assert length(f)==get_ngauss(mesh)
    out = zeros(mesh.ndof)
    @eval ccall((:FemSourceScalar_forward_Julia, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}), $out, $f)
    out
end

"""
    compute_fem_source_term(f1::Union{PyObject,Array{Float64,2}}, f2::Union{PyObject,Array{Float64,2}}, mesh::Mesh)
"""
function compute_fem_source_term(f1::Union{PyObject,Array{Float64,1}}, f2::Union{PyObject,Array{Float64,1}}, mesh::Mesh)
    [compute_fem_source_term1(f1, mesh); compute_fem_source_term1(f2, mesh)]
end

function compute_fem_source_term(f::Union{PyObject,Array{Float64,1}},  mmesh::Mesh)
    @assert length(f)==2*get_ngauss(mmesh)
    compute_fem_source_term(f[1:get_ngauss(mmesh)], f[get_ngauss(mmesh)+1:end], mmesh)
end

"""
    compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh)
"""
function compute_fem_laplace_matrix1(kappa::PyObject, mesh::Mesh)
    @assert length(kappa) == get_ngauss(mesh)
    fem_laplace_scalar_ = load_op_and_grad(AdFem.libmfem,"fem_laplace_scalar", multiple=true)
    kappa = convert_to_tensor(Any[kappa], [Float64]); kappa = kappa[1]
    indices, vv = fem_laplace_scalar_(kappa)
    n = mesh.ndof
    RawSparseTensor(indices, vv, n, n)
end

"""
    compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh)
"""
function compute_fem_laplace_matrix1(kappa::Array{Float64,1}, mesh::Mesh)
    @assert length(kappa) == get_ngauss(mesh)
    N = get_ngauss(mesh) * size(mesh.conn, 2)^2
    indices = zeros(Int64, 2N)
    vv = zeros(N)
    @eval ccall((:FemLaplaceScalar_forward_Julia, $LIBMFEM), Cvoid, (Ptr{Int64}, Ptr{Cdouble}, Ptr{Cdouble}), $indices, $vv, $kappa)
    indices = reshape(indices, 2, N)'|>Array
    sparse(indices[:,1] .+ 1, indices[:,2] .+ 1, vv, mesh.ndof, mesh.ndof)
end

"""
    compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)
"""
function compute_fem_laplace_matrix(kappa::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)
    Z = compute_fem_laplace_matrix1(kappa, mesh)
    if isa(Z, SparseMatrixCSC)
        [Z spzeros(mesh.ndof, mesh.ndof)
        spzeros(mesh.ndof, mesh.ndof) Z]
    else
        [Z spzero(mesh.ndof)
        spzero(mesh.ndof) Z]
    end
end

compute_fem_laplace_matrix1(mesh::Mesh) = compute_fem_laplace_matrix1(ones(get_ngauss(mesh)), mesh)
compute_fem_laplace_matrix(mesh::Mesh) = compute_fem_laplace_matrix(ones(get_ngauss(mesh)), mesh)

"""
    fem_impose_Dirichlet_boundary_condition1(L::SparseTensor, bdnode::Array{Int64}, mesh::Mesh)

A differentiable kernel for imposing the Dirichlet boundary of a scalar-valued function.
"""
function fem_impose_Dirichlet_boundary_condition1(L::SparseTensor, bdnode::Array{Int64}, mesh::Mesh)
    @warn "Consider imposing boundary conditions using an algebraic approach: impose_Dirichlet_boundary_conditions" maxlog=1
    idx = bdnode
    Lbd = L[:, idx]
    L = scatter_update(L, :, idx, spzero(size(mesh.nodes, 1), length(idx)))
    L = scatter_update(L, idx, :,  spzero(length(idx), size(mesh.nodes, 1)))
    L = scatter_update(L, idx, idx, spdiag(length(idx)))
    L, Lbd
end

"""
    compute_interaction_term(p::Union{PyObject,Array{Float64, 1}}, mesh::Mesh)
"""
function compute_interaction_term(p::Union{PyObject,Array{Float64, 1}}, mesh::Mesh)
    compute_interaction_term_mfem_ = load_op_and_grad(AdFem.libmfem,"compute_interaction_term_mfem")
    p = convert_to_tensor(Any[p], [Float64]); p = p[1]
    out = compute_interaction_term_mfem_(p)
    set_shape(out, (2mesh.ndof, ))
end

"""
    compute_fem_mass_matrix1(rho::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)
"""
function compute_fem_mass_matrix1(rho::Union{PyObject, Array{Float64, 1}}, mesh::Mesh)
    compute_fem_mass_matrix_mfem_ = load_op_and_grad(AdFem.libmfem,"compute_fem_mass_matrix_mfem", multiple=true)
    rho = convert_to_tensor(Any[rho], [Float64]); rho = rho[1]
    indices, vals = compute_fem_mass_matrix_mfem_(rho)
    n = mesh.ndof
    A = RawSparseTensor(indices, vals, n, n)
    A
end

function compute_fem_mass_matrix1(mmesh::Mesh)
    compute_fem_mass_matrix1(ones(get_ngauss(mmesh)), mmesh)
end

function compute_fem_mass_matrix(mmesh::Mesh)
    C = compute_fem_mass_matrix1(mmesh)
    Z = spzero(size(C, 1))
    [C Z;Z C]
end

function compute_fem_mass_matrix(ρ::Union{PyObject, Array{Float64, 1}}, mmesh::Mesh)
    C = compute_fem_mass_matrix1(ρ, mmesh)
    Z = spzero(size(C, 1))
    [C Z;Z C]
end

"""
    compute_fem_advection_matrix1(u::Union{Array{Float64,1}, PyObject},v::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
    compute_fem_advection_matrix1(u::Array{Float64,1}, v::Array{Float64,1}, mesh::Mesh)
"""
function compute_fem_advection_matrix1(u::Union{Array{Float64,1}, PyObject},v::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
    @assert length(u)==length(v)==get_ngauss(mesh)
    compute_fem_advection_matrix_mfem_ = load_op_and_grad(AdFem.libmfem,"compute_fem_advection_matrix_mfem", multiple=true)
    u,v = convert_to_tensor(Any[u,v], [Float64,Float64])
    indices, vals = compute_fem_advection_matrix_mfem_(u,v)
    n = mesh.ndof
    RawSparseTensor(indices, vals, n, n)
end

function compute_fem_advection_matrix1(u::Array{Float64,1}, v::Array{Float64,1}, mesh::Mesh)
    @assert length(u)==length(v)==get_ngauss(mesh)
    N = get_ngauss(mesh) * mesh.elem_ndof^2
    indices = zeros(Int64, 2N)
    vv = zeros(N)
    @eval ccall((:ComputeFemAdvectionMatrixMfem_forward_Julia, $LIBMFEM),
        Cvoid, (Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), $indices, $vv, $u, $v)
    sparse(indices[1:2:end] .+ 1, indices[2:2:end] .+ 1, vv, mesh.ndof, mesh.ndof)
end

"""
    compute_interaction_matrix(mesh::Mesh)
"""
function compute_interaction_matrix(mesh::Mesh)
    elem_dof = size(mesh.conn, 2)
    N = get_ngauss(mesh) * 2 * elem_dof
    ii = zeros(Int64, N)
    jj = zeros(Int64, N)
    vv = zeros(Float64, N)
    @eval ccall((:ComputeInteractionMatrixMfem, $LIBMFEM), Cvoid, (Ptr{Int64}, Ptr{Int64}, Ptr{Cdouble}), $ii, $jj, $vv)
    m = size(mesh.elems, 1)
    sparse(ii, jj, vv, m, 2mesh.ndof)
end

"""
    eval_grad_on_gauss_pts1(u::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
"""
function eval_grad_on_gauss_pts1(u::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
    @assert length(u)==mesh.ndof
    fem_grad_mfem_ = load_op_and_grad(AdFem.libmfem,"fem_grad_mfem")
    u = convert_to_tensor(Any[u], [Float64]); u = u[1]
    out = fem_grad_mfem_(u)
    m = size(gauss_nodes(mesh), 1)
    set_shape(out, (m, 2))
end

"""
    eval_grad_on_gauss_pts(u::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
"""
# function eval_grad_on_gauss_pts(u::Union{Array{Float64,1}, PyObject}, mesh::Mesh)
#     n = size(mesh.nodes, 1)
#     m = size(gauss_nodes(mesh), 1)
#     r1 = eval_grad_on_gauss_pts1(u[1:n], mesh)
#     r2 = eval_grad_on_gauss_pts1(u[n+1:end], mesh)
#     out = zeros(m, 2, 2)
#     for i = 1:m
#         out[i, 1, 1] = r1[i, 1]  # out is numerical array but r1 is tensor, cannot assign value
#         out[i, 1, 2] = r1[i, 2]
#         out[i, 2, 1] = r2[i, 1]
#         out[i, 2, 2] = r2[i, 2]
#     end
#     return out
# end

@doc raw"""
    compute_fem_stiffness_matrix(kappa::PyObject, mesh::Mesh)
    compute_fem_stiffness_matrix(kappa::Array{Float64, 3}, mesh::Mesh)
    compute_fem_stiffness_matrix(kappa::Array{Float64, 2}, mesh::Mesh)

Computes the stiffness matrix. Here, the acceptable sizes of $\kappa$
- a $3\times 3$ matrix, which is the pointwise uniform stiffness matrix
- a $N_g \times 3 \times 3$ tensor, which includes a specific stiffness matrix at each Gauss node
"""
function compute_fem_stiffness_matrix(kappa::PyObject, mesh::Mesh)
    if length(size(kappa))==2
        kappa = reshape(repeat(reshape(kappa, (-1,)), get_ngauss(mesh)), (get_ngauss(mesh), 3, 3))
    end
    @assert size(kappa) == (get_ngauss(mesh), 3, 3)
    compute_fem_stiffness_matrix_mfem_ = load_op_and_grad(AdFem.libmfem,"compute_fem_stiffness_matrix_mfem", multiple=true)
    kappa = convert_to_tensor(Any[kappa], [Float64]); kappa = kappa[1]
    kappa = reshape(kappa, (-1,))
    indices, vv = compute_fem_stiffness_matrix_mfem_(kappa)
    RawSparseTensor(indices, vv, 2mesh.ndof, 2mesh.ndof)
end

function compute_fem_stiffness_matrix(kappa::Array{Float64, 2}, mesh::Mesh)
    @assert size(kappa) == (3, 3)
    K = zeros(get_ngauss(mesh), 3, 3)
    for i = 1:get_ngauss(mesh)
        K[i, :, :] = kappa
    end
    compute_fem_stiffness_matrix(K, mesh)
end

function compute_fem_stiffness_matrix(kappa::Array{Float64, 3}, mesh::Mesh)
    @assert size(kappa) == (get_ngauss(mesh), 3, 3)
    K = zeros(length(kappa))
    s = 1
    for i = 1:size(kappa, 1)
        for p = 1:3
            for q = 1:3
                K[s] = kappa[i, p, q]
                s += 1
            end
        end
    end
    N = (size(mesh.conn, 2) * 2)^2 * get_ngauss(mesh);
    indices = zeros(Int64, 2N)
    vv = zeros(N)
    @eval ccall((:ComputeFemStiffnessMatrixMfem_forward_Julia, $LIBMFEM), Cvoid,
            (Ptr{Clonglong}, Ptr{Cdouble}, Ptr{Cdouble}), $indices, $vv, $K)
    indices = reshape(indices, 2, N)'|>Array
    sparse(indices[:,1] .+ 1, indices[:, 2] .+ 1, vv, 2mesh.ndof, 2mesh.ndof)
end

@doc raw"""
    get_bdedge_integration_pts(bdedge::Array{Int64, 2}, mesh::Mesh)

Returns Gauss quadrature points on the boundary edge as a $n\times 2$ matrix. Here $n$ is the number of rows for `bdedge`.
"""
function get_bdedge_integration_pts(bdedge::Array{Int64, 2}, mmesh::Mesh)
    N = @eval ccall((:get_LineIntegralN, $(AdFem.LIBMFEM)), Cint, ())
    xs = zeros(N)
    ws = zeros(N)
    @eval ccall((:get_LineIntegralPnW, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}), $xs, $ws)

    bdN = size(bdedge, 1)
    x = zeros(N * bdN)
    y = zeros(N * bdN)
    for i = 1:bdN
        l = mmesh.nodes[bdedge[i,1],:]
        r = mmesh.nodes[bdedge[i,2],:]
        xx = @. (1-xs) * l[1] + xs * r[1]
        yy = @. (1-xs) * l[2] + xs * r[2]
        x[(i-1)*N+1:i*N] = xx
        y[(i-1)*N+1:i*N] = yy
    end
    [x y]
end

@doc raw"""
    _traction_get_nodes(bdedge::Array{Int64, 2}, mesh::Mesh)

Returns $x$ and $y$ coordinates for each element that contains `bdedge`.

For example, if `bdedge = [2,4]`, this function finds the corresponding element, e.g., one that contains vertices
2, 4, 5. Then two arrays are returned:
- `node_x`: $x$ coordinates of node 2, 4, and 5
- `node_y`: $y$ coordinates of node 2, 4, and 5
"""
function _traction_get_nodes(bdedge::Array{Int64, 2}, mesh::Mesh)
    nodes_x = zeros(size(bdedge, 1), 3)
    nodes_y = zeros(size(bdedge, 1), 3)
    _edge_to_elem_map = Dict{Tuple{Int64, Int64}, Int64}()
    dic = (x,y)->(minimum([x,y]), maximum([x,y]))
    add_dic = (e1, e2, e3)->begin
        if haskey(_edge_to_elem_map, dic(e1, e2))
            delete!(_edge_to_elem_map, dic(e1, e2))
        else
            _edge_to_elem_map[dic(e1, e2)] = e3
        end
    end
    for i = 1:mesh.nelem
        e1, e2, e3 = mesh.elems[i,:]
        add_dic(e1, e2, e3)
        add_dic(e3, e2, e1)
        add_dic(e1, e3, e2)
    end
    for i = 1:size(bdedge, 1)
        idx = _edge_to_elem_map[dic(bdedge[i,1], bdedge[i,2])]
        nodes_x[i, :] = [mesh.nodes[bdedge[i,1], 1]; mesh.nodes[bdedge[i,2], 1]; mesh.nodes[idx, 1]]
        nodes_y[i, :] = [mesh.nodes[bdedge[i,1], 2]; mesh.nodes[bdedge[i,2], 2]; mesh.nodes[idx, 2]]
    end
    nodes_x, nodes_y
end

@doc raw"""
    bcedge(mesh::Mesh)

Returns all boundary edges as a set of integer pairs (edge vertices).
"""
function bcedge(mesh::Mesh)
    _edge_to_elem_map = Set{Tuple{Int64, Int64}}()
    dic = (x,y)->(minimum([x,y]), maximum([x,y]))
    add_dic = (e1, e2)->begin
        if dic(e1, e2) in _edge_to_elem_map
            delete!(_edge_to_elem_map, dic(e1, e2))
        else
            push!(_edge_to_elem_map, dic(e1, e2))
        end
    end
    for i = 1:mesh.nelem
        e1, e2, e3 = mesh.elems[i,:]
        add_dic(e1, e2)
        add_dic(e3, e2)
        add_dic(e1, e3)
    end
    out = zeros(Int64, length(_edge_to_elem_map), 2)
    for (k,s) in enumerate(_edge_to_elem_map)
        out[k,:] = [s[1] s[2]]
    end
    out
end

@doc raw"""
    bcedge(f::Function, mesh::Mesh)

Returns all edge indices that satisfies `f(x1, y1, x2, y2) = true`
Here the edge endpoints are given by $(x_1, y_1)$ and $(x_2, y_2)$.
"""
function bcedge(f::Function, mmesh::Mesh)
    out = bcedge(mmesh)
    edges = []
    for i = 1:size(out, 1)
        e1, e2 = out[i,:]
        x1, y1 = mmesh.nodes[e1, :]
        x2, y2 = mmesh.nodes[e2, :]
        if f(x1, y1, x2, y2)
            push!(edges, [e1 e2])
        end
    end
    vcat(edges...)
end

@doc raw"""
    bcnode(mesh::Mesh; by_dof::Bool = true)

Returns all boundary node indices.

If `by_dof = true`, `bcnode` returns the global indices for boundary DOFs.

- For `P2` elements, the returned values are boundary node DOFs + boundary edge DOFs (offseted by `mesh.nnode`)
- For `BDM1` elements, the returned values are boundary edge DOFs + boundary edge DOFs offseted by `mesh.nedge`
"""
function bcnode(mmesh::Mesh; by_dof::Bool = true)
    bdedge = bcedge(mmesh)
    if by_dof && mmesh.elem_type == P2
        edgedof = get_edge_dof(bdedge, mmesh) .+ mmesh.nnode
        [collect(Set(bdedge[:])); edgedof]
    elseif by_dof && mmesh.elem_type == BDM1
        bd = get_edge_dof(bdedge, mmesh)
        [bd; bd .+ mmesh.nedge]
    else
        collect(Set(bdedge[:]))
    end
end

"""
    bcnode(f::Function, mesh::Mesh; by_dof::Bool = true)

Returns the boundary node DOFs that satisfies `f(x,y) = true`.

!!! note

    For BDM1 element and `by_dof = true`, because the degrees of freedoms are associated with edges, `f` has the signature

    ```julia
    f(x1::Float64, y1::Float64, x2::Float64, y2::Float64)::Bool
    ```

    `bcnode` only returns DOFs on edges such that `f(x1, y1, x2, y2)=true`.
"""
function bcnode(f::Function, mesh::Mesh; by_dof::Bool = true)
    if mesh.elem_type==BDM1
        return _bcnode_bdm1(f, mesh; by_dof=by_dof)
    end
    @assert mesh.elem_type in [P1, P2]
    nd = bcnode(mesh, by_dof=by_dof)
    out = Int64[]
    for i = 1:length(nd)
        if nd[i]<=mesh.nnode
            if f(mesh.nodes[nd[i], 1], mesh.nodes[nd[i],2])
                push!(out, nd[i])
            end
        else
            a, b = mesh.edges[nd[i]-mesh.nnode, :]
            xy = (mesh.nodes[a, :] + mesh.nodes[b, :])/2
            if f(xy[1], xy[2])
                push!(out, nd[i])
            end
        end
    end
    out
end

function _bcnode_bdm1(f::Function, mmesh::Mesh; by_dof::Bool = true)
    if !by_dof
        return bcnode(mmesh, by_dof = false)
    end
    bdedge = bcedge((x1, y1, x2, y2)->f(x1, y1, x2, y2), mmesh)
    if length(bdedge)==0
        return Int64[]
    end
    e = get_edge_dof(bdedge, mmesh)
    [e; e.+mmesh.nedge]
end

@doc raw"""
    compute_fem_traction_term1(t::Array{Float64, 1},
    bdedge::Array{Int64,2}, mesh::Mesh)

Computes the boundary integral

$$\int_{\Gamma} t(x, y) \delta u dx$$

Returns a vector of size `dof`.
"""
function compute_fem_traction_term1(t::Array{Float64, 1},
            bdedge::Array{Int64,2}, mesh::Mesh)
    # sort bdedge so that bdedge[i,1] < bdedge[i,2]
    out = zeros(mesh.ndof)
    bdN = size(bdedge, 1)
    if mesh.elem_type == P1
        @eval ccall((:ComputeFemTractionV2_Julia, $LIBMFEM), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Clonglong}, Cint), $out, $t, $(bdedge'[:]), $(bdN))
        return out
    end
    for i = 1:size(bdedge, 1)
        if bdedge[i,1]>bdedge[i,2]
            bdedge[i,:] = [bdedge[i,2]; bdedge[i,1]]
        end
    end
    order = mesh.lorder
    D = _edge_dict(mesh)
    node_x = zeros(size(bdedge, 1))
    node_y = zeros(size(bdedge, 1))
    ngauss = Int64(@eval ccall((:ComputeFemTractionTermMfem_forward_getNGauss, $LIBMFEM), Cint,
            (Cint,), Int32($order)))

    @assert length(t) == ngauss * bdN
    if mesh.elem_type==P1
        dofs = zeros(Int64, 2bdN)
        for i = 1:bdN
            dofs[2*i-1] = bdedge[i,1] - 1
            dofs[2*i] = bdedge[i,2] - 1
        end
    elseif mesh.elem_type==P2
        dofs = zeros(Int64, 3bdN)
        for i = 1:bdN
            dofs[3*i-2] = bdedge[i,1] - 1
            dofs[3*i-1] = bdedge[i,2] - 1
            dofs[3*i] = D[(bdedge[i,1], bdedge[i,2])] + mesh.nnode - 1
        end
    elseif mesh.elem_type==BDM1
        dofs = zeros(Int64, 2bdN)
        for i = 1:bdN
            e = D[(bdedge[i,1], bdedge[i,2])]
            dofs[2*i-1] = e - 1
            dofs[2*i] = e - 1 + mesh.nedge
        end
    end

    bdnode_x, bdnode_y = _traction_get_nodes(bdedge, mesh)
    bdnode_x = bdnode_x'[:]
    bdnode_y = bdnode_y'[:]

    if mesh.elem_type == P2
        @eval ccall((:ComputeFemTractionTermMfem_forward_Julia, $LIBMFEM), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint),
                $out, $t, Int32.($dofs), $bdnode_x, $bdnode_y, Int32($bdN), Int32($order))
    elseif mesh.elem_type == BDM1
        sn = get_boundary_edge_orientation(bdedge, mesh)
        @eval ccall((:ComputeBDMTractionTermMfem_forward_Julia, $LIBMFEM), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint),
            $out, $sn, $t, Int32.($dofs), $bdnode_x, $bdnode_y, Int32($bdN), Int32($order))
    end
    return out
end

@doc raw"""
    compute_fem_traction_term1(t::PyObject,bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the boundary integral

$$\int_{\Gamma} t(x, y) \delta u dx$$

Returns a vector of size `dof`.

Here $t$ is defined on boundary Gauss points, e.g., the quantity computed from [`eval_f_on_boundary_edge`](@ref).

!!! info
    Currently, only P1 element is supported.

# Example

Here is an example to evaluate the boundary integral on a domain $[0,1]^2$

$$\int_{0}^1 t(x, 1) \delta u dx$$

```julia
mmesh = Mesh(10,10,0.1)
a = constant(1.0)
f = (x,y)->x+y*a
top = bcedge((x1,y1,x2,y2)->(y1>0.99) && (y2>0.99), mmesh)
t = eval_f_on_boundary_edge(f, top, mmesh; tensor_input = true)
T = compute_fem_traction_term1(t, top, mmesh)
```
"""
function compute_fem_traction_term1(t::PyObject,bdedge::Array{Int64, 2}, mmesh::Mesh)
    N = @eval ccall((:get_LineIntegralN, $(AdFem.LIBMFEM)), Cint, ())
    if size(bdedge, 1)*N != length(t)
        @error "Size of `t` (=$(length(t))) does not match #edge x #gauss points per edge (=$(size(bdedge, 1)) × $N)"
    end
    @assert mmesh.elem_type == P1
    compute_fem_traction_term_v_ = load_op_and_grad(AdFem.libmfem,"compute_fem_traction_term_v")
    t,edgeid = convert_to_tensor(Any[t,bdedge], [Float64,Int64])
    out = compute_fem_traction_term_v_(t,edgeid)
    set_shape(out, (mmesh.ndof,))
end

@doc raw"""
    compute_fem_traction_term(t1::PyObject, t2::PyObject, bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the boundary integral

$$\int_\Gamma \mathbf{t}(x,y)\cdot \delta \mathbf{u} d\mathbf{x}$$

Here $\mathbf{t}(x,y)$ is the boundary traction term

$$\mathbf{t}(x,y) = \begin{bmatrix} t_1(x,y) \\ t_2(x,y) \end{bmatrix}$$

and

$$\mathbf{u} = \begin{bmatrix} u_1(x,y) \\ u_2(x,y) \end{bmatrix}$$

$t_1$ and $t_2$ are defined on boundary Gauss points. See [`eval_f_on_boundary_edge`](@ref).
"""
function compute_fem_traction_term(t1::PyObject, t2::PyObject, bdedge::Array{Int64, 2}, mmesh::Mesh)
    T1 = compute_fem_traction_term1(t1, bdedge, mmesh)
    T2 = compute_fem_traction_term1(t2, bdedge, mmesh)
    [T1;T2]
end

"""
    compute_fem_traction_term(t::Array{Float64, 2},
    bdedge::Array{Int64,2}, mesh::Mesh)
"""
function compute_fem_traction_term(t::Array{Float64, 2},
    bdedge::Array{Int64,2}, mesh::Mesh)
    @assert size(t,2)==2
    [compute_fem_traction_term1(t[:,1], bdedge, mesh);
    compute_fem_traction_term1(t[:,2], bdedge, mesh)]
end

"""
    compute_fem_traction_term(t1::Array{Float64, 1}, t2::Array{Float64, 1},
    bdedge::Array{Int64,2}, mesh::Mesh)
"""
function compute_fem_traction_term(t1::Array{Float64, 1}, t2::Array{Float64, 1},
    bdedge::Array{Int64,2}, mesh::Mesh)
    [compute_fem_traction_term1(t1, bdedge, mesh);
    compute_fem_traction_term1(t2, bdedge, mesh)]
end

"""
    eval_f_on_boundary_edge(f::Function, bdedge::Array{Int64, 2}, mesh::Mesh; tensor_input::Bool = false)

Evaluates `f` on the boundary **Gauss points**. Here `f` has the signature

```f(Float64, Float64)::Float64```

or

```f(PyObject, PyObject)::PyObject```
"""
function eval_f_on_boundary_edge(f::Function, bdedge::Array{Int64, 2},
        mesh::Mesh; tensor_input::Bool = false)
    pts = get_bdedge_integration_pts(bdedge, mesh)
    if tensor_input
        f(constant(pts[:,1]), constant(pts[:,2]))
    else
        f.(pts[:,1], pts[:,2])
    end
end

"""
    eval_f_on_boundary_node(f::Function, bdnode::Array{Int64}, mesh::Mesh)
"""
function eval_f_on_boundary_node(f::Function, bdnode::Array{Int64}, mesh::Mesh)
    out = zeros(length(bdnode))
    for i = 1:length(bdnode)
        if bdnode[i]>mesh.nnode
            a, b = mesh.edges[bdnode[i]-mesh.nnode, :]
            xy = (mesh.nodes[a, :] + mesh.nodes[b, :])/2
            out[i] = f(xy[1], xy[2])
        else
            a = bdnode[i]
            out[i] = f(mesh.nodes[a, 1], mesh.nodes[a, 2])
        end
    end
    out
end

@doc raw"""
    compute_von_mises_stress_term(K::Array{Float64, 3}, u::Array{Float64, 1}, mesh::Mesh)
    compute_von_mises_stress_term(K::Array{Float64, 2}, u::Array{Float64, 1}, mesh::Mesh)

Computes the von Mises stress term

$$\sigma_v = \sqrt{\sigma_1^2 - \sigma_1 \sigma_2 + \sigma_2^2 + 3\sigma_{12}^2}$$
"""
function compute_von_mises_stress_term(K::Array{Float64, 3}, u::Array{Float64, 1}, mesh::Mesh)
    @assert length(u) == 2mesh.ndof
    @assert size(K) == (get_ngauss(mesh), 3, 3)
    hmat = zeros(9*get_ngauss(mesh))
    for i = 1:get_ngauss(mesh)
        for j = 1:3
            for k = 1:3
                hmat[(i-1)*9+(j-1)*3+k] = K[i, j, k]
            end
        end
    end
    out = zeros(get_ngauss(mesh))
    @eval ccall((:mfem_compute_von_mises_stress_term, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                $out, $hmat, $u)
    out
end

function compute_von_mises_stress_term(K::Array{Float64, 2}, u::Array{Float64, 1}, mesh::Mesh)
    Ks = zeros(get_ngauss(mesh), 3, 3)
    for i = 1:get_ngauss(mesh)
        Ks[i,:,:] = K
    end
    compute_von_mises_stress_term(Ks, u, mesh)
end

"""
    compute_fem_laplace_term1(u::Array{Float64, 1},nu::Array{Float64, 1}, mesh::Mesh)
"""
function compute_fem_laplace_term1(u::Array{Float64, 1},nu::Array{Float64, 1}, mesh::Mesh)
    @assert length(u) == mesh.ndof
    @assert length(nu) == get_ngauss(mesh)
    out = zeros(mesh.ndof)
    @eval ccall((:ComputeLaplaceTermMfem_forward_Julia, $LIBMFEM), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), $out, $nu, $u)
    out
end

"""
    compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
                                nu::Union{PyObject, Array{Float64, 1}},
                                mesh::Mesh)
"""
function compute_fem_laplace_term1(u::Union{PyObject, Array{Float64, 1}},
                                   nu::Union{PyObject, Array{Float64, 1}},
                                   mesh::Mesh)
    @assert length(u) == mesh.ndof
    @assert length(nu) == get_ngauss(mesh)
    compute_laplace_term_mfem_ = load_op_and_grad(AdFem.libmfem,"compute_laplace_term_mfem")
    u,nu = convert_to_tensor(Any[u,nu], [Float64,Float64])
    out = compute_laplace_term_mfem_(u,nu)
    set_shape(out, (mesh.ndof,))
end

"""
    gauss_weights(mmesh::Mesh)

Returns the weights for each Gauss points.
"""
function gauss_weights(mmesh::Mesh)
    w = zeros(get_ngauss(mmesh))
    @eval ccall((:mfem_get_gauss_weights, $LIBMFEM), Cvoid, (Ptr{Cdouble},), $w)
    w
end

"""
    compute_fvm_source_term(f::Array{Float64, 1}, mmesh::Mesh)
"""
function compute_fvm_source_term(f::Array{Float64, 1}, mmesh::Mesh)
    w = gauss_weights(mmesh)
    src = zeros(mmesh.nelem)
    ngauss_per_elem = get_ngauss(mmesh)÷mmesh.nelem
    for i = 1:mmesh.nelem
        idx = (i-1)*ngauss_per_elem+1:i*ngauss_per_elem
        src[i] = sum(f[idx].*w[idx])
    end
    src
end

@doc raw"""
    compute_strain_energy_term(Sigma::Array{Float64, 2}, mmesh::Mesh)

Computes the strain energy term

$$\int_A \sigma : \varepsilon (\delta u) dx$$

Here $\sigma$ is a fourth-order tensor. `Sigma` is a `ngauss × 3` matrix, each row represents  $[\sigma_{11}, \sigma_{22}, \sigma_{12}]$ at
each Gauss point.

The output is a length `2mmesh.ndof` vector.
"""
function compute_strain_energy_term(Sigma::Array{Float64, 2}, mmesh::Mesh)
    @assert size(Sigma,1)==get_ngauss(mmesh)
    @assert size(Sigma,2)==3
    out = zeros(2mmesh.ndof)
    @eval ccall((:ComputeStrainEnergyTermMfem_forward_Julia, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}), $out, $(Array(Sigma')))
    out
end

"""
    compute_strain_energy_term(Sigma::PyObject, mmesh::Mesh)
"""
function compute_strain_energy_term(Sigma::PyObject, mmesh::Mesh)
    @assert size(Sigma,1)==get_ngauss(mmesh)
    @assert size(Sigma,2)==3
    compute_strain_energy_term_mfem_ =  @eval load_op_and_grad($libmfem,"compute_strain_energy_term_mfem")
    sigma = convert_to_tensor(Any[Sigma], [Float64]); sigma = sigma[1]
    se = compute_strain_energy_term_mfem_(sigma)
    set_shape(se, (2mmesh.ndof,))
end

"""
    eval_strain_on_gauss_pts(u::Array{Float64}, mmesh::Mesh)

Evaluates the strain on Gauss points. `u` is a vector of size `2mmesh.ndof`.

The output is a `ngauss × 3` vector.
"""
function eval_strain_on_gauss_pts(u::Array{Float64}, mmesh::Mesh)
    @assert length(u)==2mmesh.ndof
    ε = zeros(3get_ngauss(mmesh))
    @eval ccall((:EvalStrainOnGaussPts_forward_Julia, $LIBMFEM), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}), $ε, $u)
    Array(reshape(ε, 3, get_ngauss(mmesh))')
end

"""
    eval_strain_on_gauss_pts(u::PyObject, mmesh::Mesh)
"""
function eval_strain_on_gauss_pts(u::PyObject, mmesh::Mesh)
    @assert length(u)==2mmesh.ndof
    eval_strain_on_gauss_pts_ = @eval load_op_and_grad($libmfem,"eval_strain_on_gauss_pts")
    u = convert_to_tensor(Any[u], [Float64]); u = u[1]
    ε = eval_strain_on_gauss_pts_(u)
    set_shape(ε, (get_ngauss(mmesh), 3))
end

@doc raw"""
    compute_fem_boundary_mass_matrix1(c::Union{Array{Float64}, PyObject}, bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the matrix

$$\int_\Gamma cu \delta u ds$$

The parameters are
- `bdedge`: a $N_e \times 2$ integer array, the boundary edge to integrate on
- `c`: given by a vector of length $4N_e$; currently, each edge has 4 quadrature points;

The output is a $N_v\times N_v$ sparse matrix.
"""
function compute_fem_boundary_mass_matrix1(c::Union{Array{Float64}, PyObject}, bdedge::Array{Int64, 2}, mmesh::Mesh)
    @assert length(c)==4*size(bdedge, 1)
    compute_boundary_mass_matrix_one_ = @eval load_op_and_grad($libmfem,"compute_boundary_mass_matrix_one", multiple=true)
    c,idx = convert_to_tensor(Any[c,bdedge], [Float64,Int64])
    ij, vv = compute_boundary_mass_matrix_one_(c,idx)
    RawSparseTensor(ij, vv, mmesh.ndof, mmesh.ndof)
end

@doc raw"""
    compute_fem_boundary_mass_term1(u::Union{Array{Float64}, PyObject},
        c::Union{Array{Float64}, PyObject}, bdedge::Array{Int64, 2}, mmesh::Mesh)

Computes the term

$$\int_\Gamma cu \delta u ds$$

The parameters are
- `u` : a vector of length $N_v$
- `bdedge` : a $N_e \times 2$ integer array, the boundary edge to integrate on
- `c`: given by a vector of length $4N_e$; currently, each edge has 4 quadrature points;

The output is a $N_v\times N_v$ sparse matrix.
"""
function compute_fem_boundary_mass_term1(u::Union{Array{Float64}, PyObject},
        c::Union{Array{Float64}, PyObject}, bdedge::Array{Int64, 2}, mmesh::Mesh)
    @assert length(u)==mmesh.nnode
    compute_fem_boundary_mass_matrix1(c, bdedge, mmesh) * u
end

@doc raw"""
    eval_scalar_on_boundary_edge(u::Union{PyObject, Array{Float64, 1}},
            edge::Array{Int64, 1}, mmesh::Mesh)

Returns an array of values on the Gauss quadrature nodes for each edge.

- `u`: nodal values
- `edge`: a $N\times 2$ integer array; each row represents an edge $(x_1, x_2)$

The returned array consists of $(y_1, y_2, \ldots)$

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/AdFem/docs/eval_scalar_on_boundary_edge.png?raw=true)
"""
function eval_scalar_on_boundary_edge(u::Union{PyObject, Array{Float64, 1}},
        edge::Array{Int64, 2}, mmesh::Mesh)
    eval_scalar_on_boundary_edge_ = @eval load_op_and_grad($libmfem,"eval_scalar_on_boundary_edge")
    u,edge = convert_to_tensor(Any[u,edge], [Float64,Int64])
    out = eval_scalar_on_boundary_edge_(u,edge)
    N = @eval ccall((:get_LineIntegralN, $(AdFem.LIBMFEM)), Cint, ())
    reshape(out, (N*size(edge, 1), ))
end

@doc raw"""
    eval_strain_on_boundary_edge(u::Union{PyObject, Array{Float64, 1}},
    edge::Array{Int64, 2}, mmesh::Mesh)

Returns an array of strain tensors on the Gauss quadrature nodes for each edge.

The returned value has size $N_e N_g\times 3$. Here $N_e$ is the number of edges, and $N_g$ is the number of
Gauss points on the edge. Each row of `edge`, $(l,r)$, has the following property: $l < r$
"""
function eval_strain_on_boundary_edge(u::Union{PyObject, Array{Float64, 1}},
    edge::Array{Int64, 2}, mmesh::Mesh)
    @assert length(u) == 2mmesh.ndof
    eval_strain_on_boundary_edge_ = @eval load_op_and_grad($libmfem,"eval_strain_on_boundary_edge")
    u,edge = convert_to_tensor(Any[u,edge], [Float64,Int64])
    out = eval_strain_on_boundary_edge_(u,edge)
    N = @eval ccall((:get_LineIntegralN, $(AdFem.LIBMFEM)), Cint, ())
    reshape(out, (N*size(edge, 1), 3))
end

@doc raw"""
    eval_normal_and_shear_stress_on_boundary_edge(sigma::Union{PyObject, Array{Float64, 2}},
        edge::Array{Int64, 2}, mmesh::Mesh)

Calculates the normal and shear stress on boundary edges.

$$\sigma_n = (\sigma n)\cdot n, \tau = (\sigma n) \cdot m$$

Here
$$m = \begin{bmatrix}n_2\\-n_1\end{bmatrix}$$

![](https://github.com/ADCMEMarket/ADCMEImages/blob/master/AdFem/docs/shear_normal.png?raw=true)
"""
function eval_normal_and_shear_stress_on_boundary_edge(sigma::Union{PyObject, Array{Float64, 2}},
    edge::Array{Int64, 2}, mmesh::Mesh)
    nls = get_edge_normal(edge, mmesh)
    N = @eval ccall((:get_LineIntegralN, $(AdFem.LIBMFEM)), Cint, ())
    normals = zeros(N*size(nls, 1), 2)
    for i = 1:size(nls, 1)
        normals[(i-1)*N+1:i*N, :] = repeat(nls[i,:]', N, 1)
    end
    eval_normal_and_shear_stress_on_boundary_edge_ = @eval load_op_and_grad($libmfem,"eval_normal_and_shear_stress_on_boundary_edge", multiple=true)
    sigma,normals = convert_to_tensor(Any[sigma,normals], [Float64,Float64])
    sn, st = eval_normal_and_shear_stress_on_boundary_edge_(sigma,normals)
    set_shape(sn, (N*size(nls, 1),)), set_shape(st, (N*size(nls, 1), ))
end


text:
function DiffusionDefinition.diffusion_message(::Val{:SIR})
    println("* * * * * * * * * * * * * * * * * * * * * * * * ")
    println("Susceptible-Infectious-Recovered (SIR) diffusion")
    println("* * * * * * * * * * * * * * * * * * * * * * * * ")
    println("Description")
    println("------------")
    println("A two dimensional diffusion (Iₜ,Rₜ):")
    println("dIₜ=[α(1-Iₜ-Rₜ)Iₜ - βIₜ]dt -σ₁√[(1-Iₜ-Rₜ)Iₜ]dWₜ,")
    println("dIₜ=βIₜdt -σ₂√IₜdWₜ,")
    println("corresponding to proportion of Infectious (Iₜ) and Recovered (Rₜ).")
    println("a proportion of susceptible is computed via: Sₜ=1-Iₜ-Rₜ.")
    println("Iₜ,Rₜ,Sₜ∈[0,1].")
    println("The model has four parameters:")
    println("\tα, β, σ1, σ2")
    println()
    println("Constructors")
    println("------------")
    println("You may initialize the process either by passing the parameters")
    println("in the order above, as in:")
    println("α, β, σ₁, σ₂ = ...")
    println("P = SIR(α, β, σ₁, σ₂)")
    println("or by keyword arguments")
    println("p = (α = ..., β = ..., σ1 = ..., σ2 = ...)")
    println("P = SIR(;p...)")
    println()
end

@diffusion_process SIR{T} begin
    :dimensions
    process --> 2
    wiener --> 2

    :parameters
    (α, β, σ1, σ2) --> T

    #=
    :conjugate
    phi(t, u) --> (
        (0.0, 0.0),
        ((1.0 - u[1] - u[2])*u[1], 0.0),
        (-u[1], u[1]),
        (0.0, 0.0),
        (0.0, 0.0),
        (0.0, 0.0),
    )
    nonhypo(x) --> x
    num_non_hypo --> 2
    =#

    :additional
    statespace --> BoundedStateSpace(((1,2), (0.0,0.0)), ((1,2), (1.0,1.0)))
end

function DiffusionDefinition.b(t, x, P::SIR)
    @SVector[
        P.α*(1.0 - x[1] - x[2])*x[1] - P.β*x[1],
        P.β*x[1]
    ]
end

function DiffusionDefinition.σ(t, x, P::SIR)
    @SMatrix [
        -P.σ1*sqrt((1.0 - x[1] - x[2])*x[1])  -P.σ2*sqrt(x[1]);
        0.0   P.σ2*sqrt(x[1])
    ]
end

function DiffusionDefinition.bound_satisfied(P::SIR, x)
    (x[1] + x[2]) < 1.0 && DiffusionDefinition._bound_satisfied(P, x)
end

ext: jl


text:
using Test

# Just test to see that Lyceum can be loaded.
using Lyceum

@testset "Lyceum.jl" begin
    @test true
end


text:
@doc raw"""
    fit!(
        model::StateSpaceModel;
        filter::KalmanFilter=default_filter(model),
        optimizer::Optimizer=Optimizer(Optim.LBFGS())
    )

Estimate the state-space model parameters via maximum likelihood. The resulting optimal
hyperparameters and the corresponding log-likelihood are stored within the model. You can
choose the desired filter method (`UnivariateKalmanFilter`, `ScalarKalmanFilter`, etc.) and
the [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl) optimization algortihm.

# Example
```jldoctest
julia> model = LocalLevel(rand(100))
LocalLevel model

julia> fit!(model)
LocalLevel model

julia> model = LocalLinearTrend(LinRange(1, 100, 100) + rand(100))
LocalLinearTrend model

julia> fit!(model; optimizer = Optimizer(StateSpaceModels.Optim.NelderMead()))
LocalLinearTrend model
```
"""
function fit!(
    model::StateSpaceModel;
    filter::KalmanFilter=default_filter(model),
    optimizer::Optimizer=default_optimizer(model),
)
    isfitted(model) && return model
    @assert has_fit_methods(typeof(model))
    initial_unconstrained_hyperparameter = handle_optim_initial_hyperparameters(model)
    # TODO Should there be a try catch?
    func = TwiceDifferentiable(
        x -> -optim_loglike(model, filter, x), initial_unconstrained_hyperparameter
    )
    opt = optimize(
        func, initial_unconstrained_hyperparameter, optimizer.method, optimizer.options
    )
    # optim_loglike returns the mean loglike for numerical stability purposes
    # the optimum loglike must be rescaled
    opt_loglikelihood = -opt.minimum * size(model.system.y, 1)
    opt_hyperparameters = opt.minimizer
    update_model_hyperparameters!(model, opt_hyperparameters)
    # TODO
    # I leaned that this is not a good way to compute the covariance matrix of paarameters
    # we should investigate other methods
    numerical_hessian = Optim.hessian!(func, opt_hyperparameters)
    std_err = diag(pinv(numerical_hessian))
    fill_results!(model, opt_loglikelihood, std_err)
    return model
end

"""
    has_fit_methods(model_type::Type{<:StateSpaceModel}) -> Bool

Verify if a certain `StateSpaceModel` has the necessary methods to perform `fit!``.
"""
function has_fit_methods(model_type::Type{<:StateSpaceModel})
    tuple_with_model_type = Tuple{model_type}
    m1 = hasmethod(default_filter, tuple_with_model_type)
    m2 = hasmethod(initial_hyperparameters!, tuple_with_model_type)
    m3 = hasmethod(constrain_hyperparameters!, tuple_with_model_type)
    m4 = hasmethod(unconstrain_hyperparameters!, tuple_with_model_type)
    m5 = hasmethod(fill_model_system!, tuple_with_model_type)
    return m1 && m2 && m3 && m4 && m5
end

struct CoefficientTable{Fl<:AbstractFloat}
    names::Vector{String}
    coef::Vector{Fl}
    std_err::Vector{Fl}
    z::Vector{Fl}
    p_value::Vector{Fl}

    function CoefficientTable{Fl}(
        names::Vector{String},
        coef::Vector{Fl},
        std_err::Vector{Fl},
        z::Vector{Fl},
        p_value::Vector{Fl},
    ) where Fl
        @assert length(names) ==
                length(coef) ==
                length(std_err) ==
                length(z) ==
                length(p_value)
        return new{Fl}(names, coef, std_err, z, p_value)
    end
end

function CoefficientTable{Fl}() where Fl
    return CoefficientTable{Fl}(String[], Fl[], Fl[], Fl[], Fl[])
end

Base.length(coef_table::CoefficientTable) = length(coef_table.names)

function Base.isempty(coef_table::CoefficientTable)
    return isempty(coef_table.names) &&
           isempty(coef_table.coef) &&
           isempty(coef_table.std_err) &&
           isempty(coef_table.z) &&
           isempty(coef_table.p_value)
end

mutable struct Results{Fl<:AbstractFloat}
    model_name::String
    coef_table::CoefficientTable{Fl}
    llk::Fl
    aic::Fl
    aicc::Fl
    bic::Fl
    num_observations::Int
    num_hyperparameters::Int
end

function Results{Fl}() where Fl
    return Results{Fl}("", CoefficientTable{Fl}(), Fl(NaN), Fl(NaN), Fl(NaN), Fl(NaN), 0, 0)
end

"""
    results(model::StateSpaceModel)

Query the results of the optimization called by `fit!`.
"""
results(model::StateSpaceModel) = model.results
function Base.isempty(results::Results)
    return isempty(results.coef_table) &&
           isnan(results.llk) &&
           isnan(results.aic) &&
           isnan(results.aicc) &&
           isnan(results.bic) &&
           iszero(results.num_observations) &&
           iszero(results.num_hyperparameters)
end

function fill_results!(model::StateSpaceModel, llk::Fl, std_err::Vector{Fl}) where Fl
    n_obs = length(model.system.y)
    num_hyperparameters = number_free_hyperparameters(model)
    coef_table = build_coef_table(model, std_err)
    # Fill results
    model.results.model_name = model_name(model)
    model.results.coef_table = coef_table
    model.results.llk = llk
    model.results.aic = aic(num_hyperparameters, llk)
    model.results.aicc = aicc(n_obs, num_hyperparameters, llk)
    model.results.bic = bic(n_obs, num_hyperparameters, llk)
    model.results.num_observations = n_obs
    model.results.num_hyperparameters = num_hyperparameters
    return model
end

function aic(k::Int, llk::Fl) where Fl
    return convert(Fl, 2 * k - 2 * llk)
end

function aicc(n::Int, k::Int, llk::Fl) where Fl
    return convert(Fl, aic(k, llk) + (2 * k * (k + 1) / (n - k - 1)))
end

function bic(n::Int, k::Int, llk::Fl) where Fl
    return convert(Fl, log(n) * k - 2 * llk)
end

aic(model::StateSpaceModel) = model.results.aic
aicc(model::StateSpaceModel) = model.results.aicc
bic(model::StateSpaceModel) = model.results.bic

function build_coef_table(model::StateSpaceModel, std_err::Vector{Fl}) where Fl
    all_coef = get_constrained_values(model)
    all_std_err = handle_std_err(model, std_err)
    all_z_stat = handle_z_stat(all_coef, all_std_err)
    all_p_value = handle_p_value(all_coef, all_std_err, all_z_stat)

    return CoefficientTable{Fl}(
        get_names(model), all_coef, all_std_err, all_z_stat, all_p_value
    )
end

function handle_std_err(model::StateSpaceModel, std_err::Vector{Fl}) where Fl
    all_std_err = fill(NaN, number_hyperparameters(model))
    # Put the std_err in the correct position
    for i in 1:length(std_err)
        hyperparamter_on_minimizer = model.hyperparameters.minimizer_hyperparameter_position[i]
        position_on_list_of_hyperparameters = position(
            hyperparamter_on_minimizer, model.hyperparameters
        )
        all_std_err[position_on_list_of_hyperparameters] = std_err[i]
    end
    return Fl.(all_std_err)
end

function handle_z_stat(all_coef::Vector{Fl}, all_std_err::Vector{Fl}) where Fl
    all_z_stat = fill(NaN, length(all_std_err))
    for i in 1:length(all_std_err)
        if !isnan(all_std_err[i])
            all_z_stat[i] = all_coef[i] / all_std_err[i]
        end
    end
    return Fl.(all_z_stat)
end

function handle_p_value(
    all_coef::Vector{Fl}, all_std_err::Vector{Fl}, all_z_stat::Vector{Fl}
) where Fl
    all_p_value = fill(NaN, length(all_std_err))
    for i in 1:length(all_std_err)
        if !isnan(all_std_err[i]) && all_std_err[i] > 0
            dist = Normal(all_coef[i], all_std_err[i])
            all_p_value[i] = 1 - 2 * abs(cdf(dist, all_z_stat[i]) - 0.5)
        end
    end
    return Fl.(all_p_value)
end


text:
# Thrift types are supported
# All Hive types (complex types) are not supported yet

const JTYPES = Dict(
  Int32(0)  => Bool,                # BOOLEAN
  Int32(1)  => Int8,                # TINYINT
  Int32(2)  => Int16,               # SMALLINT
  Int32(3)  => Int32,               # INT
  Int32(4)  => Int64,               # BIGINT
  Int32(5)  => Float32,             # FLOAT
  Int32(6)  => Float64,             # DOUBLE
  Int32(7)  => String,              # STRING
  Int32(8)  => DateTime,            # TIMESTAMP
  Int32(9)  => String,              # BINARY
  Int32(10) => String,              # ARRAY
  Int32(11) => String,              # MAP
  Int32(12) => String,              # STRUCT
  Int32(13) => String,              # UNIONTYPE
  Int32(15) => BigFloat,            # DECIMAL
  Int32(16) => Nothing,             # NULL
  Int32(17) => Date,                # DATE
  Int32(18) => String,              # VARCHAR
  Int32(19) => String,              # CHAR                  (TODO: can be optimized for CHAR(1))
  Int32(20) => String,              # INTERVAL_YEAR_MONTH
  Int32(21) => String               # INTERVAL_DAY_TIME
)

#=
function with_null_check(fn, str::String)
    str = strip(str)
    isempty(str) ? nothing : fn(str)
end
=#

tochar(str::String) = isempty(str) ? Char(0) : first(str)
tobigfloat(str::String) = BigFloat(String(strip(str)))
todate(str::String) = Date(strip(str))
toint8(val::UInt8) = reinterpret(Int8, val)
todatetime(ftime) = Dates.unix2datetime(ftime)
toint32(str::String) = parse(Int32, str)
toint64(str::String) = parse(Int64, str)
tofloat32(str::String) = parse(Float32, str)
tofloat64(str::String) = parse(Float64, str)
function todatetime(str::String)
    str = strip(str)
    occursin('-', str) ? DateTime(replace(str, " "=>"T")) : Dates.unix2datetime(parse(Int, str))
end

const JCONV = Dict(
  Int32(0)  => nothing,             # BOOLEAN
  Int32(1)  => toint8,              # TINYINT
  Int32(2)  => nothing,             # SMALLINT
  Int32(3)  => nothing,             # INT
  Int32(4)  => nothing,             # BIGINT
  Int32(5)  => Float32,             # FLOAT
  Int32(6)  => nothing,             # DOUBLE
  Int32(7)  => nothing,             # STRING
  Int32(8)  => todatetime,          # TIMESTAMP
  Int32(9)  => nothing,             # BINARY
  Int32(10) => nothing,             # ARRAY
  Int32(11) => nothing,             # MAP
  Int32(12) => nothing,             # STRUCT
  Int32(13) => nothing,             # UNIONTYPE
  Int32(15) => tobigfloat,          # DECIMAL
  Int32(16) => nothing,             # NULL
  Int32(17) => todate,              # DATE
  Int32(18) => nothing,             # VARCHAR
  Int32(19) => nothing,             # CHAR
  Int32(20) => nothing,             # INTERVAL_YEAR_MONTH
  Int32(21) => nothing              # INTERVAL_DAY_TIME
)

##
# map column types to Julia types
function julia_type(hs2type::TTypeDesc, typeentry::TPrimitiveTypeEntry)
    if (typeentry._type == 19) && isfilled(typeentry, :typeQualifiers) && isfilled(typeentry.typeQualifiers, :qualifiers)
        qual = typeentry.typeQualifiers.qualifiers
        ("characterMaximumLength" in keys(qual)) && (qual["characterMaximumLength"].i32Value == 1) && (return Char)
    elseif typeentry._type == 15 && isfilled(typeentry, :typeQualifiers) && isfilled(typeentry.typeQualifiers, :qualifiers)
        qual = typeentry.typeQualifiers.qualifiers
        precision = ("precision" in keys(qual)) ? qual["precision"].i32Value : Int32(10)
        scale = ("scale" in keys(qual)) ? qual["scale"].i32Value : Int32(0)
        if scale == 0
            (precision <= 9) && (return Int32)
            (precision <= 18) && (return Int64)
        else
            ((precision+scale) <= 8) && (return Float32)
            ((precision+scale) <= 17) && (return Float64)
        end
    end
    JTYPES[typeentry._type]
end

function julia_conv(hs2type::TTypeDesc, typeentry::TPrimitiveTypeEntry)
    if (typeentry._type == 19) && isfilled(typeentry, :typeQualifiers) && isfilled(typeentry.typeQualifiers, :qualifiers)
        ("characterMaximumLength" in keys(typeentry.typeQualifiers.qualifiers)) && (typeentry.typeQualifiers.qualifiers["characterMaximumLength"].i32Value == 1) && (return tochar)
    elseif typeentry._type == 15 && isfilled(typeentry, :typeQualifiers) && isfilled(typeentry.typeQualifiers, :qualifiers)
        qual = typeentry.typeQualifiers.qualifiers
        precision = ("precision" in keys(qual)) ? qual["precision"].i32Value : Int32(10)
        scale = ("scale" in keys(qual)) ? qual["scale"].i32Value : Int32(0)
        if scale == 0
            (precision <= 9) && (return toint32)
            (precision <= 18) && (return toint64)
        else
            (precision <= 7) && (return tofloat32)
            (precision <= 16) && (return tofloat64)
        end
    end
    JCONV[typeentry._type]
end

function julia_type(hs2type::TTypeDesc, typeentry::TArrayTypeEntry)
    typeptr = typeentry.objectTypePtr
    T = julia_type(hs2type, hs2type.types[typeptr])
    Vector{T}
end

function julia_type(hs2type::TTypeDesc, typeentry::TMapTypeEntry)
    keytypeptr = typeentry.keyTypePtr
    valuetypeptr = typeentry.valueTypePtr
    K = julia_type(hs2type, hs2type.types[keytypeptr])
    V = julia_type(hs2type, hs2type.types[valuetypeptr])
    Dict{K,V}
end

julia_type(hs2type::TTypeDesc, typeentry::TStructTypeEntry) = Any
julia_type(hs2type::TTypeDesc, typeentry::TUnionTypeEntry) = Any
julia_type(hs2type::TTypeDesc, typeentry::TUserDefinedTypeEntry) = Any

function julia_type(hs2type::TTypeDesc, typeentry::TTypeEntry)
    for fld in fieldnames(TTypeEntry)
        isfilled(typeentry, fld) && (return julia_type(hs2type, getfield(typeentry, fld)))
    end
end

function julia_conv(hs2type::TTypeDesc, typeentry::TTypeEntry)
    for fld in fieldnames(TTypeEntry)
        isfilled(typeentry, fld) && (return julia_conv(hs2type, getfield(typeentry, fld)))
    end
end

function julia_type(hs2type::TTypeDesc)
    # The "top" type is always the first element of the list.
    # If the top type is an ARRAY, MAP, STRUCT, or UNIONTYPE type, then subsequent elements represent nested types.
    toptype = hs2type.types[1]
    julia_type(hs2type, toptype)
end

function julia_conv(hs2type::TTypeDesc)
    # The "top" type is always the first element of the list.
    # If the top type is an ARRAY, MAP, STRUCT, or UNIONTYPE type, then subsequent elements represent nested types.
    toptype = hs2type.types[1]
    julia_conv(hs2type, toptype)
end

##
# map column types to Julia type names
julia_type_name(hs2type::TTypeDesc, typeentry::TPrimitiveTypeEntry) = TYPE_NAMES[typeentry._type]

function julia_type_name(hs2type::TTypeDesc, typeentry::TArrayTypeEntry)
    typeptr = typeentry.objectTypePtr
    typename = julia_type_name(hs2type, hs2type.types[typeptr])
    "Vector{$typename}"
end

function julia_type_name(hs2type::TTypeDesc, typeentry::TMapTypeEntry)
    keytypeptr = typeentry.keyTypePtr
    valuetypeptr = typeentry.valueTypePtr
    keytypename = julia_type_name(hs2type, hs2type.types[keytypeptr])
    valuetypename = julia_type_name(hs2type, hs2type.types[valuetypeptr])
    "Dict{$keytypename, $valuetypename}"
end

julia_type_name(hs2type::TTypeDesc, typeentry::TStructTypeEntry) = "Struct"
julia_type_name(hs2type::TTypeDesc, typeentry::TUnionTypeEntry) = "Union"
julia_type_name(hs2type::TTypeDesc, typeentry::TUserDefinedTypeEntry) = "UserDefined"

function julia_type_name(hs2type::TTypeDesc, typeentry::TTypeEntry)
    for fld in fieldnames(TTypeEntry)
        isfilled(typeentry, fld) || continue
        return julia_type_name(hs2type, getfield(typeentry, fld))
    end
end

function julia_type_name(hs2type::TTypeDesc)
    # The "top" type is always the first element of the list.
    # If the top type is an ARRAY, MAP, STRUCT, or UNIONTYPE type, then subsequent elements represent nested types.
    toptype = hs2type.types[1]
    julia_type_name(hs2type, toptype)
end

##
# map column values to Julia types
# map columns into DataVectors
const ColValue = Union{TBoolValue, TByteValue, TI16Value, TI32Value, TI64Value, TDoubleValue, TStringValue}
const Col = Union{TBoolColumn, TByteColumn, TI16Column, TI32Column, TI64Column, TDoubleColumn, TStringColumn, TBinaryColumn}

julia_type(colval::T, convfn::Nothing) where T<:ColValue = isfilled(colval, :value) ? getfield(colval, :value) : nothing
julia_type(colval::T, convfn) where T<:ColValue = isfilled(colval, :value) ? convfn(getfield(colval, :value)) : nothing

function julia_type(colval::TColumnValue, convfn)
    for fld in fieldnames(TColumnValue)
        isfilled(colval, fld) && (return julia_type(getfield(colval, fld), convfn))
    end
    nothing
end

function julia_type(col::T, convfn) where T<:Col
    if (length(col.nulls) < 2) || isempty(col.values)
        map(convfn, col.values)
    else
        nulls = bitset_to_bools(col.nulls, length(col.values))
        data_with_missings(map(convfn, col.values), nulls)
    end
end

function julia_type(col::T, convfn::Nothing) where T<:Col
    if (length(col.nulls) < 2) || isempty(col.values)
        col.values
    else
        nulls = bitset_to_bools(col.nulls, length(col.values))
        data_with_missings(col.values, nulls)
    end
end

function julia_type(col::TColumn, convfn)
    for fld in fieldnames(TColumn)
        isfilled(col, fld) && (return julia_type(getfield(col, fld), convfn))
    end
    nothing
end


text:
function get_link_eq_constraints(graph::ModelGraph)
    link_cons = LinkConstraint[]
    for link_con in getlinkconstraints(graph)
        if isa(link_con.set,MOI.EqualTo)
            push!(link_cons,link_con)
        end
    end
    return link_cons
end

function get_link_ineq_constraints(graph::ModelGraph)
    link_cons = LinkConstraint[]
    for link_con in getlinkconstraints(graph)
        if !(isa(link_con.set,MOI.EqualTo))
            push!(link_cons,link_con)  #It must be an inequality constraint
        end
    end
    return link_cons
end

function prepare_link_matrix(link_constraints::Vector{LinkConstraint},sub_var_map::Dict)
    link_map = Dict()
    link_vars = []   #model node variables
    n_link_vars = 0
    I = Int64[]
    J = Int64[]
    V = Float64[]
    b = Float64[]
    i = 0
    for link in link_constraints  #Just enumerate here
        i += 1
        _unpack_constraint!(link,i,I,J,V,b,link_vars,link_map,sub_var_map)
    end
    A = sparse(I,J,V)

    vars = [sub_var_map[var] for var in link_vars]  #convert to subproblem variables
    return A,b,vars
end

function prepare_link_matrix(link_variables::Vector{LinkVariableRef},master_var_map::Dict)
    link_map = Dict()
    link_vars = []
    n_link_vars = 0
    I = Int64[]
    J = Int64[]
    V = Float64[]
    b = Float64[]
    row = 1
    for link_var in link_variables
        row += _unpack_constraint!(link_var,row,I,J,V,link_vars,link_map,master_var_map)
    end
    H = sparse(I,J,V)
    return H,link_vars #b is always zero for this matrix
end

function _unpack_constraint!(linkvariable::LinkVariableRef,row::Integer,I::Vector,J::Vector,V::Vector,link_vars::Vector,link_map::Dict,master_var_map::Dict)

    child_vars = master_var_map[linkvariable]
    n_sub_problems = length(child_vars)

    j = length(link_vars)
    for var in child_vars
        push!(link_vars,var)    #NOTE: There shouldn't be duplicate child vars from different link variables.      #A linked_variable should only have one parent
        j = length(link_vars)
        link_map[var] = j       #link index of variable var
    end

    for i = 1:n_sub_problems - 1

        var1 = child_vars[i]
        j1 = link_map[var1]             #look up the variable index since it was already added to the matrix

        var2 = child_vars[i+1]
        j2 = link_map[var2]             #look up the variable index since it was already added to the matrix

        coeff1 = -1
        coeff2 = 1

        push!(I,row)
        push!(J,j1)
        push!(V,coeff1)

        push!(I,row)
        push!(J,j2)
        push!(V,coeff2)

        row += 1
    end

    return row - 1
end

function _unpack_constraint!(constraint::LinkConstraint{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},MOI.EqualTo{Float64}},row::Integer,I::Vector,J::Vector,V::Vector,b::Vector,link_vars::Vector,link_map::Dict,var_map::Dict)
    push!(b,constraint.set.value)
    for (var,coeff) in constraint.func.terms
        #check for new link constraint variables
        if !(var in link_vars)
            j = length(link_vars) + 1
            link_map[var] = j   #link index of variable var
            push!(link_vars,var)
        else
            #look up the variable index since it was already added to the matrix
            j = link_map[var]
        end
        push!(I,row)
        push!(J,j)
        push!(V,coeff)
    end
end

function _unpack_constraint!(constraint::LinkConstraint{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},MOI.LessThan{Float64}},row::Integer,I::Vector,J::Vector,V::Vector,b::Vector,link_vars::Vector,link_map::Dict)
    push!(b,constraint.set.upper)
    for (var,coeff) in constraint.func.terms
        #find a new link variable
        if !(var in link_vars)
            j = length(link_vars) + 1
            link_map[var] = j   #link index of variable var
            push!(link_vars,var)
        else
            #get correct variable index
            j = link_map[var]
        end
        push!(I,row)
        push!(J,j)
        push!(V,coeff)
    end
end

function _unpack_constraint!(constraint::LinkConstraint{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},MOI.GreaterThan{Float64}},row::Integer,I::Vector,J::Vector,V::Vector,b::Vector,link_vars::Vector,link_map::Dict)
    push!(b,-1*constraint.set.lower)
    for (var,coeff) in constraint.func.terms
        #find a new link variable
        if !(var in link_vars)
            j = length(link_vars) + 1
            link_map[var] = j   #link index of variable var
            push!(link_vars,var)
        else
            #get correct variable index
            j = link_map[var]
        end
        push!(I,row)
        push!(J,j)
        push!(V,-1*coeff)
    end
end

#TODO: Interval Constraints

#child_vars = linkvariable.graph.linkvariable_map[linkvariable]   #the children variables for this linkvariable

# j = length(link_vars)
# for var in child_vars
#     push!(link_vars,var)    #NOTE: There shouldn't be duplicate child vars from different link variables.      #A linked_variable should only have one parent
#     j = length(link_vars)
#     link_map[var] = j       #link index of variable var
# end

# for i = 1:length(child_vars) - 1    #for each child variable
#     var1 = child_vars[i]
#     j1 = link_map[var1]             #look up the variable index since it was already added to the matrix
#
#     var2 = child_vars[i+1]
#     j2 = link_map[var2]             #look up the variable index since it was already added to the matrix
#
#     coeff1 = -1
#     coeff2 = 1
#
#     push!(I,row)
#     push!(J,j1)
#     push!(V,coeff1)
#
#     push!(I,row)
#     push!(J,j2)
#     push!(V,coeff2)
#
#     row += 1
# end

#Add entry for master variable (link variable)
#master_var = linkvariable.vref                                  #the actual master JuMP variable
#push!(link_vars,master_var)
#coeff_master = 1
#j_master = length(link_vars)
#link_map[master_var] = j_master


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

using Markdown
using InteractiveUtils

# ╔═╡ ec794b0c-67b1-11eb-23e7-459641028308
begin
	using Markdown
	using InteractiveUtils
	using Test
	using LinearAlgebra
end

# ╔═╡ 47f1758c-67b1-11eb-0a89-27d2595e5912
md"""

## Task \# 8

Compute the norms and condition numbers

#### Norm

```math
||A||^2 = max \frac{||Ax||^2}{||x||^2} = max \frac{ x^T A^T A x }{ x^T x } = \lambda_{max}(A^T A) = \sigma^{2}_{max}
```

```math
||A|| = \sqrt{\sigma^{2}_{max}} = \sigma_{max}
```

```math
||A^{-1}|| = \frac{1}{\sigma_{min}}
```

"""

# ╔═╡ 9aa05ff6-67b5-11eb-0b36-95986372b105
function normCalculator(A::Array{Int})
	AᵀA = A' * A

	A_eigenvals, A_eigenvec = eigen(AᵀA)

	A_normFromEigenSquared = √max(A_eigenvals...)

	A_norm = opnorm(A)

	return (A_normFromEigenSquared, A_norm)
end

# ╔═╡ 4533fdcc-67b3-11eb-30f8-851ef7fc43c9
begin
	A₁ = [
		1 7
		1 1
	]

	A₁normFromEigenSquared, A₁norm = normCalculator(A₁)

	@test round.(A₁normFromEigenSquared, digits=5) == round.(A₁norm, digits=5)

	A₁normFromEigenSquared, A₁norm, cond(A₁)
end

# ╔═╡ bf65578a-67b8-11eb-2627-9139089fcf68
md"""
#### Cond is Inf, because of `` || A_2^{-1} || = \sigma_{min} = 0 ``
"""

# ╔═╡ 1c2ed1ac-67b5-11eb-1e77-67e96dbd2802
begin
	A₂ = [
		1 1
		0 0
	]

	A₂normFromEigenSquared, A₂norm = normCalculator(A₂)

	@test round.(A₂normFromEigenSquared, digits=5) == round.(A₂norm, digits=5)

	A₂normFromEigenSquared, A₂norm, cond(A₂)
end

# ╔═╡ 1568a204-67b9-11eb-3d55-d72928803e5e
md"""
#### Cond is 1, because of `` || A_2 || = || A_2^{-1} || ``
"""

# ╔═╡ ed7d86f6-67b3-11eb-273c-43770e6f299c
begin
	A₃ = [
		1 1
		-1 1
	]

	A₃normFromEigenSquared, A₃norm = normCalculator(A₃)

	@test round.(A₃normFromEigenSquared, digits=5) == round.(A₃norm, digits=5)

	A₃normFromEigenSquared, A₃norm, cond(A₃)
end

# ╔═╡ c3dc9efe-67bd-11eb-0224-af9d70653b73
md"""
### Task \# 13
Estimate the condition number
"""

# ╔═╡ e21e095c-67bd-11eb-12e4-e92b18b0b223
begin
	A₄ = [
		1 1
		1 1.0001
	]
	cond(A₄)
end

# ╔═╡ 67060602-67c2-11eb-046e-dfc437a3dc81
md"""
### Task \#17
Fibonacci matrix and the SVD
"""

# ╔═╡ 832167be-67c2-11eb-1693-ab8ea1a74995
begin
	F = [
		1 1
		1 0
	]
	svd(F)
end

# ╔═╡ 8d7f475a-67c3-11eb-243b-53c5d68004f1
md"""
### Task \#18
"""

# ╔═╡ 9a4b6df6-67c3-11eb-27d6-2546dd023fb4
begin
	K = [
		1 2
		2 5
	]
	lu(K)
end

# ╔═╡ Cell order:
# ╠═ec794b0c-67b1-11eb-23e7-459641028308
# ╟─47f1758c-67b1-11eb-0a89-27d2595e5912
# ╠═9aa05ff6-67b5-11eb-0b36-95986372b105
# ╠═4533fdcc-67b3-11eb-30f8-851ef7fc43c9
# ╟─bf65578a-67b8-11eb-2627-9139089fcf68
# ╠═1c2ed1ac-67b5-11eb-1e77-67e96dbd2802
# ╟─1568a204-67b9-11eb-3d55-d72928803e5e
# ╠═ed7d86f6-67b3-11eb-273c-43770e6f299c
# ╟─c3dc9efe-67bd-11eb-0224-af9d70653b73
# ╠═e21e095c-67bd-11eb-12e4-e92b18b0b223
# ╟─67060602-67c2-11eb-046e-dfc437a3dc81
# ╠═832167be-67c2-11eb-1693-ab8ea1a74995
# ╟─8d7f475a-67c3-11eb-243b-53c5d68004f1
# ╠═9a4b6df6-67c3-11eb-27d6-2546dd023fb4


text:
export WeightedMeasure, AbstractWeightedMeasure

"""
    struct WeightedMeasure{R,M} <: AbstractMeasure
        logweight :: R
        base :: M
    end
"""

abstract type AbstractWeightedMeasure <: AbstractMeasure end

logweight(μ::AbstractWeightedMeasure) = μ.logweight
basemeasure(μ::AbstractWeightedMeasure) = μ.base

function logdensity(sm::AbstractWeightedMeasure, x)
    logdensity(sm.base, x) + sm.logweight
end

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

struct WeightedMeasure{R,M} <: AbstractWeightedMeasure
    logweight::R
    base::M
end

function Base.show(io::IO, μ::WeightedMeasure)
    io = IOContext(io, :compact => true)
    print(io, exp(μ.logweight), " * ", μ.base)
end

function Base.show_unquoted(io::IO, μ::WeightedMeasure, indent::Int, prec::Int)
    io = IOContext(io, :compact => true)
    if Base.operator_precedence(:*) ≤ prec
        print(io, "(")
        show(io, μ)
        print(io, ")")
    else
        show(io, μ)
    end
    return nothing
end

function Base.:*(k::T, m::AbstractMeasure) where {T<:Number}
    logk = log(k)
    return weightedmeasure(logk, m)
end

Base.:*(m::AbstractMeasure, k::Real) = k * m

≪(::M, ::WeightedMeasure{R,M}) where {R,M} = true
≪(::WeightedMeasure{R,M}, ::M) where {R,M} = true

sampletype(μ::WeightedMeasure) = sampletype(μ.base)

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

export ParamWeightedMeasure
struct ParamWeightedMeasure{L,N,T,B} <: AbstractWeightedMeasure
    ℓ::L
    par::NamedTuple{N,T}
    base::B
end

function Base.show(io::IO, d::ParamWeightedMeasure)
    io = IOContext(io, :compact => true)
    print(io, "ParamWeighted(", d.ℓ, ", ", d.par, ", ", d.base, ")")
end

basemeasure(d::ParamWeightedMeasure) = d.base

logdensity(d::ParamWeightedMeasure, x) = d.ℓ(d.par)


text:
using FactCheck
using FunctionalCollections

import FunctionalCollections: SparseBitmappedTrie, SparseNode, SparseLeaf,
    bitpos, index, hasindex, arrayof, update

facts("Sparse Bitmapped Vector Tries") do

    context("bitpos and index") do
        l = SparseLeaf{Int}([6, 11, 16, 21, 26],
                            2^5 | 2^10 | 2^15 | 2^20 | 2^25)
        @fact index(l, 6) --> 1
        @fact index(l, 11) --> 2
        @fact index(l, 16) --> 3
        @fact index(l, 21) --> 4
        @fact index(l, 26) --> 5

        @fact hasindex(l, 5)  --> false
        @fact hasindex(l, 6)  --> true
        @fact hasindex(l, 7)  --> false
        @fact hasindex(l, 11) --> true
        @fact hasindex(l, 12) --> false
        @fact hasindex(l, 21) --> true
    end

    context("SparseLeaf update") do
        l = SparseLeaf{Int}()

        @fact update(l, 1, 1)   --> (leaf) -> hasindex(leaf[1], 1)
        @fact update(l, 10, 10) --> (leaf) -> hasindex(leaf[1], 10) &&
                                             !hasindex(leaf[1], 11)

        l = SparseLeaf{Int}([1, 5], 2^0 | 2^4)
        l, _ = update(l, 2, 2)
        @fact arrayof(l) --> [1, 2, 5]
        @fact length(l) --> 3
        @fact hasindex(l, 1) --> true
        @fact hasindex(l, 2) --> true
        @fact hasindex(l, 3) --> false
        @fact hasindex(l, 4) --> false
        @fact hasindex(l, 5) --> true
        @fact hasindex(l, 6) --> false

        @fact arrayof(update(l, 2, 100)[1]) --> [1, 100, 5]
    end

    context("SparseNode update") do
        n, _ = update(SparseNode(AbstractString), 1, "foo")
        @fact length(arrayof(n)) --> 1

        leaf = arrayof(n)[1].arr[1].arr[1].arr[1].arr[1].arr[1].arr[1]
        @fact hasindex(leaf, 1) --> true
        @fact leaf.arr[1] --> "foo"

        n2, _ = update(n, 33, "bar")
        leaf2 = n2.arr[1].arr[1].arr[1].arr[1].arr[1].arr[1].arr[2]
        @fact hasindex(leaf2, 33) --> true
        @fact arrayof(leaf2)[1] --> "bar"
    end

    context("SparseBitmappedTrie get") do
        n, _ = update(SparseNode(Int), 33, 33)
        @fact get(n, 33, "missing") --> 33
        @fact get(update(SparseNode(Int), 12345, 12345)[1], 12345, "missing") --> 12345
        @fact get(n, 12345, "missing") --> "missing"
    end

    context("SparseBitmappedTrie length") do
        n = SparseNode(Int)
        for i=1:1000
            n, _ = update(n, i, i)
        end
        @fact length(n) --> 1000
    end

    context("SparseBitmappedTrie items") do
        n = SparseNode(Int)
        for i=1:1000
            n, _ = update(n, i, i)
        end
        @fact [i for i=n] --> collect(1:1000)
    end

end


text:
using BenchmarkTools

module test
  using Dates
  import Random
  const N = 10^10

  function f0()
    s = 0.
    for i in 1:N
      s += rand()
    end
    return s
  end

  function f1()
    RNG1 = Random.MersenneTwister()
    s = 0.
    for i in 1:N
      s += rand(RNG1)
    end
    return s
  end

  const RNG = Random.MersenneTwister()
  function f2()
    s = 0.
    for i in 1:N
      s += rand(RNG)
    end
    return s
  end

  function t(func)
    nthreads = Threads.nthreads()
    dates = Vector{DateTime}(undef,nthreads)
    task = Vector{Task}(undef,nthreads)
    for i in 1:nthreads
      task[i] = Threads.@spawn func()
      dates[i] = Dates.now()
    end
    s = 0.
    for i in 1:nthreads
      s += fetch(task[i])
    end
    #println(s)
  end
end

println("f1: RNG inside f")
@btime test.t(test.f1)


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

struct PolarisDetector
    hostname::String
end

export PolarisDetector

struct PolarisDetectorInstance
    hostname::String
    data_io::TCPSocket
    control_io::TCPSocket
end

export PolarisDetectorInstance

function Base.open(device_spec::PolarisDetector)
    data_io = connect(device_spec.hostname, 11503)
    control_io = try
        connect(device_spec.hostname, 11502)
    catch err
        close(data_io)
        rethrow()
    end
    PolarisDetectorInstance(device_spec.hostname, data_io, control_io)
end

function Base.close(device::PolarisDetectorInstance)
    close(device.data_io)
    close(device.control_io)
end

function exec_dev_cmd(device::PolarisDetectorInstance, cmd::AbstractString, expected_response::AbstractString)
    @debug "Sending command \"$cmd\"."
    println(device.control_io, cmd)
    resp = readline(device.control_io)
    @debug "Received response \"$resp\"."
    resp != expected_response && throw(ErrorException("Device command \"cmd\" resulted in invalid/unexpected response"))
    nothing
end

@inline Base.setindex!(device::PolarisDetectorInstance, x, property::Symbol) =
    Base.setindex!(device, x, Val(property))

function Base.setindex!(device::PolarisDetectorInstance, ::Tuple{}, ::Val{:reset_index})
    exec_dev_cmd(device, "ResetIndex", "Reset sync pulse timestamp index.")
    ()
end

function Base.setindex!(device::PolarisDetectorInstance, enabled::Bool, ::Val{:sync_pulse_generation})
    if enabled
        exec_dev_cmd(device, "Start", "Start sync pulse.")
    else
        exec_dev_cmd(device, "Stop", "Stop sync pulse.")
    end
    enabled
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9e8c30b8-8aac-4c90-83bb-71578c6c6fe6
begin
	using Plots,PlutoUI
end

# ╔═╡ 11e0f36c-e18f-11ec-15f3-89441c9abc04
md"# 第一个函数"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═11e0f36c-e18f-11ec-15f3-89441c9abc04
# ╠═9e8c30b8-8aac-4c90-83bb-71578c6c6fe6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: phpsmarter/plutonote

max_issues_repo_name: phpsmarter/plutonote

max_forks_repo_name: phpsmarter/plutonote


text:
# ENV["JULIA_DEBUG"] = Metatheory
using Metatheory
using Metatheory.EGraphs
using Metatheory.EGraphs: in_same_set, find_root

@testset "Merging" begin
  testexpr = :((a * 2) / 2)
  testmatch = :(a << 1)
  G = EGraph(testexpr)
  t2, _ = addexpr!(G, testmatch)
  merge!(G, t2.id, EClassId(3))
  @test in_same_set(G.uf, t2.id, EClassId(3)) == true
  # DOES NOT UPWARD MERGE
end

# testexpr = :(42a + b * (foo($(Dict(:x => 2)), 42)))

@testset "Simple congruence - rebuilding" begin
  G = EGraph()
  ec1, _ = addexpr!(G, :(f(a, b)))
  ec2, _ = addexpr!(G, :(f(a, c)))

  testexpr = :(f(a, b) + f(a, c))

  testec, _ = addexpr!(G, testexpr)

  t1, _ = addexpr!(G, :b)
  t2, _ = addexpr!(G, :c)

  c_id = merge!(G, t2.id, t1.id)
  @test in_same_set(G.uf, c_id, t1.id)
  @test in_same_set(G.uf, t2.id, t1.id)
  rebuild!(G)
  @test in_same_set(G.uf, ec1.id, ec2.id)
end

@testset "Simple nested congruence" begin
  apply(n, f, x) = n == 0 ? x : apply(n - 1, f, f(x))
  f(x) = Expr(:call, :f, x)

  G = EGraph(:a)

  t1, _ = addexpr!(G, apply(6, f, :a))
  t2, _ = addexpr!(G, apply(9, f, :a))

  c_id = merge!(G, t1.id, EClassId(1)) # a == apply(6,f,a)
  c2_id = merge!(G, t2.id, EClassId(1)) # a == apply(9,f,a)

  rebuild!(G)

  t3, _ = addexpr!(G, apply(3, f, :a))
  t4, _ = addexpr!(G, apply(7, f, :a))

  # f^m(a) = a = f^n(a) ⟹ f^(gcd(m,n))(a) = a
  @test in_same_set(G.uf, t1.id, EClassId(1)) == true
  @test in_same_set(G.uf, t2.id, EClassId(1)) == true
  @test in_same_set(G.uf, t3.id, EClassId(1)) == true
  @test in_same_set(G.uf, t4.id, EClassId(1)) == false

  # if m or n is prime, f(a) = a
  t5, _ = addexpr!(G, apply(11, f, :a))
  t6, _ = addexpr!(G, apply(1, f, :a))
  c5_id = merge!(G, t5.id, EClassId(1)) # a == apply(11,f,a)

  rebuild!(G)

  @test in_same_set(G.uf, t5.id, EClassId(1)) == true
  @test in_same_set(G.uf, t6.id, EClassId(1)) == true
end


text:
using Downloads, GeoData, NCDatasets, JLD, Printf

"""
    smooth!(A)

Smooth data contained in a matrix with one time step (CFL) of diffusion.
"""
@views function smooth!(A)
    A[2:end-1,2:end-1] .= A[2:end-1,2:end-1] .+ 1.0./4.1.*(diff(diff(A[:,2:end-1], dims=1), dims=1) .+ diff(diff(A[2:end-1,:], dims=2), dims=2))
    A[1,:]=A[2,:]; A[end,:]=A[end-1,:]; A[:,1]=A[:,2]; A[:,end]=A[:,end-1]
    return
end

"""
    mass_balance_constants(xc, yc)

Calculate the distributed mass-balance coefficients for a given spatial grid.
"""
function mass_balance_constants(xc, yc)
    b_max    = 0.15            # max. Mass balance rate
    lat_min, lat_max = 60, 80
    Xc, Yc   = [Float32(x) for x=xc,y=yc], [Float32(y) for x=xc,y=yc]
    Yc2      = Yc .- minimum(Yc); Yc2 .= Yc2/maximum(Yc2)
    grad_b   = (1.3517 .- 0.014158.*(lat_min.+Yc2*(lat_max-lat_min)))./100.0.*0.91 # Mass Bal. gradient, from doi: 10.1017/jog.2016.75
    z_ELA    = 1300.0 .- Yc2*300.0                                 # Educated guess for ELA altitude
    return grad_b, z_ELA, b_max
end

datadir = joinpath(@__DIR__, "../data")
bm_file = joinpath(datadir, "BedMachineGreenland-2017-09-20.nc")

"""
    download_bedmachine_greenland(; force=false, force_earthdata=false, verbose=true)

Try to download the Bedmachine Greenland dataset.  During the short-course
from a server of ours.  Outside the course, it is behind a password-wall and
it is a bit more complicated, try and see the error message.

It will work if:
- there are credentials in the `~/.netrc` file for `urs.earthdata.nasa.gov`, see [1]
- or *during the course* via the copy we'll provide not needing a login

Keyword-arguments:
- set `force==true` to download even if the file is present.

Returns the local path to the file.

Refs:
[1] https://nsidc.org/support/how/how-do-i-programmatically-request-data-services
[2] https://github.com/JuliaLang/Downloads.jl/issues/103#issuecomment-818939509
"""
function download_bedmachine_greenland(; force=false, force_earthdata=false, verbose=true)
    verbose && @warn "This will download the 2GB file BedMachineGreenland-2017-09-20.nc.  Hit Ctrl-C to abort."
    url_course = "https://people.ee.ethz.ch/~werderm/6458-gpu-course/BedMachineGreenland-2017-09-20.nc"
    url = "https://n5eil01u.ecs.nsidc.org/ICEBRIDGE/IDBMG4.003/1993.01.01/BedMachineGreenland-2017-09-20.nc"

    if isfile(bm_file) && filesize(bm_file)==2249644927 && !force
        @warn "File exists and is of right size, not downloading."
        return bm_file
    end

    # try the uploaded file on my web-server
    if !force_earthdata
        try
            println("") # the _download_progress deletes a line, so provide one
            Downloads.download(url_course, bm_file; progress=_download_progress)
            return bm_file
        catch e
            if !(e isa Downloads.RequestError)
                println("Encountered an unexpected error but will try different source.  Error:\n $e")
            end
        end
    end

    # try earthdata/NASA
    # uses https://github.com/JuliaLang/Downloads.jl/issues/103#issuecomment-818939509
    # which may stop working at some point...
    downloader = Downloads.Downloader()
    easy_hook = (easy, info) -> begin
        Downloads.Curl.setopt(easy, Downloads.Curl.CURLOPT_NETRC, Downloads.Curl.CURL_NETRC_OPTIONAL)
        Downloads.Curl.setopt(easy, Downloads.Curl.CURLOPT_COOKIEFILE, "")
    end
    downloader.easy_hook = easy_hook

    try
        Downloads.download(url, out; downloader=downloader, progress=_download_progress)
        return bm_file
    catch e
        if !(e isa Downloads.RequestError)
            println("Encountered an unexpected error:\n $e")
        end
    end

    # no luck, error:
    error("""
        Could not download the Greenland Bedmachine data.  This can have a few reasons:
        - you don't have (correct) credentials to https://urs.earthdata.nasa.gov/home in
          the file `~/.netrc`.
        - something else...

        Either way, it is probably easiest if you download the file by hand (you will still
        need to create a `earthdata.nasa.gov` login).  To do this, go to
        https://nsidc.org/data/idbmg4 and click "Login to Earthdata".  Once logged in visit
        https://n5eil01u.ecs.nsidc.org/ICEBRIDGE/IDBMG4.003/1993.01.01/BedMachineGreenland-2017-09-20.nc
        and save the file.

        Move the file `BedMachineGreenland-2017-09-20.nc` to the folder
        $(joinpath(@__DIR__, "../data")).
        """)
end
function _download_progress(total, now)
    # from https://github.com/timholy/ProgressMeter.jl/blob/45e562e708f70c6e805b7ce94e35e3fe8c889d50/src/ProgressMeter.jl#L590
    print("\r\u1b[K\u1b[A") # clears the line
    @printf("Download complete %1.2f%% \n", now/total *100)
end

"""
    load_data(;downscale=nothing, nx=nothing)

Load the Bedmachine data into memory at a specific resolution.  The resolution can be
chosen by setting either:
- the downscaling factor, between 1 and 329, or
- nx, between 10208 and 32 [default=96]
  - for nx=96 and 160 there are JLD files included in the repo

Further keyword args:
- use_nc: load from the nc-file even for resolutions nx=96 or 160

Note that this will automatically produce a dataset with x-grid
divisible by 32 and y-grid by 8 (as suitable for the GPU), and thus the exact
dimensions/downscaling maybe different.

Return:
- Zbed, Hice, Mask, dx, dy
"""
function load_data(;downscale=nothing, nx=96, use_nc=false)
    if downscale!=nothing && nx!=96
        error("Only choose of the input augments `downscale` and `nx`")
    end
    # GPU is most performant for these divisors of the array-size
    nx_div = 32
    ny_div = 8

    res = get_resolution(;ds=downscale, nx=nx)
    nx = res.nx32

    if !use_nc && nx in [96, 160] # get from jld files in the repo
        data = if nx == 96
            load(joinpath(datadir, "BedMachineGreenland_96_176.jld")) # ultra low res data
        else
            load(joinpath(datadir, "BedMachineGreenland_160_288.jld")) # low res data
        end
        Hice, Mask, Zbed = data["Hice"], data["Mask"], data["Zbed"]
        xc, yc, dx, dy   = data["xc"], data["yc"], data["dx"], data["dy"]
    else  # get from nc file
        if !(isfile(bm_file) && filesize(bm_file)==2249644927)
            println("""Downloading the Bedmachine Greenland dataset.  This may take a while (2GB in size).
                          Hit Ctrl+C to abort.
                    """)
            download_bedmachine_greenland(verbose=false)
        end

        bms = NCDstack(bm_file)
        Zbed = missing2nan(downscale_and_crop(bms[:bed], res))
        Hice = missing2nan(downscale_and_crop(bms[:thickness], res))
        Mask = missing2nan(downscale_and_crop(bms[:mask], res))
        # mask: 0=ocean, 1=ice-free land, 2=grounded-ice, 3=floating-ice, 4=non-Greenland land
        Mask = (Mask.==1) .| (Mask.==2)
        dx, dy = step.(dims(Zbed))
        xc, yc =  Zbed.dims[1].val, Zbed.dims[2].val
    end

    return Zbed, Hice, Mask, abs(dx), abs(dy), xc, yc
end
## to save the two JLD files, run:
# Zbed, Hice, Mask, dx, dy, xc, yc = load_data(; nx=96, use_nc=true)
# save("../data/BedMachineGreenland_$(size(Hice,1))_$(size(Hice,2)).jld", "Mask", Mask, "Hice", Hice, "Zbed", Zbed, "xc", xc, "yc", yc, "dx", dx, "dy", dy)
# Zbed, Hice, Mask, dx, dy, xc, yc = load_data(; nx=160, use_nc=true)
# save("../data/BedMachineGreenland_$(size(Hice,1))_$(size(Hice,2)).jld", "Mask", Mask, "Hice", Hice, "Zbed", Zbed, "xc", xc, "yc", yc, "dx", dx, "dy", dy)

"""
    as_geoarray(A, template::AbstractGeoArray)

Make a GeoArray with the data from `A` and the geo-info from `template`.

Kwargs:
- name: give it a name
- staggered: uses a staggered grid, i.e. of size (nx-1,ny-1)
"""
function as_geoarray(A, template::AbstractGeoArray; name=nothing, staggerd=false)
    @assert length(size(A))==2
    dd = if staggerd
        x, y = dims(template)
        (X( 0.5*(val(x)[2:end] + val(x)[1:end-1]), mode=x.mode, metadata=x.metadata),
         Y( 0.5*(val(y)[2:end] + val(y)[1:end-1]), mode=y.mode, metadata=y.metadata))
    else
        dims(template)
    end
    return GeoArray(A, dims=dd, name=name)
end

"Enumeration of ideal resolutions for different downscalings, nx, and ny"
resolutions =
[(ds = 1, nx32 = 10208, cropx = 10, ny8 = 18344, cropy = 2)
 (ds = 2, nx32 = 5088, cropx = 43, ny8 = 9168, cropy = 11)
 (ds = 3, nx32 = 3392, cropx = 44, ny8 = 6112, cropy = 12)
 (ds = 4, nx32 = 2528, cropx = 109, ny8 = 4584, cropy = 13)
 (ds = 5, nx32 = 2016, cropx = 142, ny8 = 3664, cropy = 30)
 (ds = 6, nx32 = 1696, cropx = 47, ny8 = 3056, cropy = 15)
 (ds = 7, nx32 = 1440, cropx = 144, ny8 = 2616, cropy = 40)
 (ds = 8, nx32 = 1248, cropx = 241, ny8 = 2288, cropy = 49)
 (ds = 9, nx32 = 1120, cropx = 146, ny8 = 2032, cropy = 66)
 (ds = 10, nx32 = 992, cropx = 307, ny8 = 1832, cropy = 35)
 (ds = 11, nx32 = 928, cropx = 20, ny8 = 1664, cropy = 52)
 (ds = 12, nx32 = 832, cropx = 245, ny8 = 1528, cropy = 21)
 (ds = 13, nx32 = 768, cropx = 246, ny8 = 1408, cropy = 54)
 (ds = 14, nx32 = 704, cropx = 375, ny8 = 1304, cropy = 103)
 (ds = 15, nx32 = 672, cropx = 152, ny8 = 1224, cropy = 0)
 (ds = 16, nx32 = 608, cropx = 505, ny8 = 1144, cropy = 57)
 (ds = 17, nx32 = 576, cropx = 442, ny8 = 1080, cropy = 2)
 (ds = 18, nx32 = 544, cropx = 443, ny8 = 1016, cropy = 75)
 (ds = 19, nx32 = 512, cropx = 508, ny8 = 960, cropy = 124)
 (ds = 21, nx32 = 480, cropx = 158, ny8 = 872, cropy = 54)
 (ds = 22, nx32 = 448, cropx = 383, ny8 = 832, cropy = 63)
 (ds = 24, nx32 = 416, cropx = 257, ny8 = 760, cropy = 129)
 (ds = 26, nx32 = 384, cropx = 259, ny8 = 704, cropy = 67)
 (ds = 29, nx32 = 352, cropx = 38, ny8 = 632, cropy = 46)
 (ds = 32, nx32 = 320, cropx = 9, ny8 = 568, cropy = 201)
 (ds = 35, nx32 = 288, cropx = 172, ny8 = 520, cropy = 180)
 (ds = 40, nx32 = 256, cropx = 17, ny8 = 456, cropy = 145)
 (ds = 45, nx32 = 224, cropx = 182, ny8 = 408, cropy = 30)
 (ds = 53, nx32 = 192, cropx = 94, ny8 = 344, cropy = 166)
 (ds = 63, nx32 = 160, cropx = 200, ny8 = 288, cropy = 264)
 (ds = 79, nx32 = 128, cropx = 184, ny8 = 232, cropy = 96)
 (ds = 104, nx32 = 96, cropx = 337, ny8 = 176, cropy = 145)
 (ds = 162, nx32 = 64, cropx = 11, ny8 = 112, cropy = 363)
 (ds = 329, nx32 = 32, cropx = 18, ny8 = 56, cropy = 250)]

"""
     get_resolution(;ds=nothing, nx=nothing, ny=nothing)

Return the resolution for the BedMachine Greenland dataset which
fits best with one of the desired
- downscaling (ds),
- number of x-gridpoints (nx), or
- number of y-gridpoints (ny)
such that the grid is of a good shape for performant GPU code.
Namely nx is divisible by 32 and ny by 8.

Downscale and crop the BedMachine data like so

    downscale_and_crop(A, get_resolution(...))
"""
function get_resolution(;ds=nothing, nx=nothing, ny=nothing)
    @assert sum([ds,nx,ny].==nothing)==2 "Only one option ds, nx, or ny can be used at once."
    i = if ds!=nothing
        findclosest(r -> abs(r.ds-ds), resolutions)
    elseif nx!=nothing
        findclosest(r -> abs(r.nx32-nx), resolutions)
    elseif ny!=nothing
        findclosest(r -> abs(r.ny8-ny), resolutions)
    end
    return resolutions[i]
end

"""
    findclosest(fn, collection)

Find item in `collection` for which `fn` returns
the lowest value. If there is a tie, pick the first item encountered.
Returns index.

    findclosest(x -> abs(3 - x), [0, 1, 6, 5]) # -> 2
"""
function findclosest(fn, collection)
    i = 0
    dist = typemax(fn(collection[1]))
    for (ii,it) in enumerate(collection)
        if fn(it) < dist
            dist = fn(it)
            i = ii
        end
    end
    return i
end

"Apply downscaling and cropping returned from `get_resolution`."
downscale_and_crop(A, res) = A[1 + floor(Int, res.cropx/2):res.ds:end - ceil(Int, res.cropx/2),
                               1 + floor(Int, res.cropy/2):res.ds:end - ceil(Int, res.cropy/2)]

# Convert missing to NaN and also re-do eltype of array
function missing2nan(ar::Array, T=Float32)
    ar = convert(Matrix{Union{T,Missing}}, ar)
    ar[ismissing.(ar)] .= NaN
    return convert(Matrix{T}, ar)
end
function missing2nan(ar::AbstractGeoArray, T=Float32)
    # missing2nan(ar::AbstractGeoArray, T=Float32) = replace_missing(ar, NaN)
    data = missing2nan(ar.data)
    return GeoArray(data, dims=ar.dims, name=ar.name)
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 4, 40, 20, 4, :y, Quadrilateral(8, 2), false),
  # ϕ (phi, friction angle), c (cohesion), ψ (psi, dilation angle),
  # γ (gamma, soil unit weight), E, ν (v, Poisson's ratio)
  # [ϕ, c, ψ, γ, E, ν]
  :properties => [20.0 10.0 20.0 16.0 1.0e5 0.3;],
  :x_coords => 0.0:0.25:10.0,
  :y_coords => 5.0:-0.25:0.0,
  :tol => 0.001,
  :limit => 200,
  :incs => 25,
  :nbo2 => 8,       # No of elements to be rigidly displaced in bearing capacity
  :qs => 20,        # Surface surcharge
  :presc => -0.001  # Magnitude of incremental vertical displacements rigid footing
)

data |> display
println()

res_df, g_coord, g_num, disp = p63_skyline(data)
println()

res_df


text:
using DeepUncertainty

# MC Dropout block
function mc_conv_bn(
    kernelsize,
    inplanes,
    outplanes,
    dropout_rate = 0.2,
    activation = relu;
    rev = false,
    initβ = Flux.zeros32,
    initγ = Flux.ones32,
    ϵ = 1f-5,
    momentum = 1f-1,
    kwargs...,
)
    layers = []

    if rev
        activations = (conv = activation, bn = identity)
        bnplanes = inplanes
    else
        activations = (conv = identity, bn = activation)
        bnplanes = outplanes
    end

    dropout = (x; k...) -> Flux.dropout(x, dropout_rate; k...)
    conv_layer =
        Conv(kernelsize, Int(inplanes) => Int(outplanes), activations.conv; kwargs...)
    mclayer = MCLayer(conv_layer, dropout)
    push!(layers, mclayer)
    push!(
        layers,
        BatchNorm(
            Int(bnplanes),
            activations.bn;
            initβ = initβ,
            initγ = initγ,
            ϵ = ϵ,
            momentum = momentum,
        ),
    )

    return rev ? reverse(layers) : layers
end

cat_channels(x, y) = cat(x, y; dims = 3)


text:
using Caesar, RoME, IncrementalInference
using RoMEPlotting, Gadfly, KernelDensityEstimatePlotting
using KernelDensityEstimate
using DocStringExtensions
using DelimitedFiles, JLD

include(joinpath(@__DIR__,"slamUtils.jl"))
include(joinpath(@__DIR__,"plotSASUtils.jl"))

expnum = 2;
if expnum == 1
    topdir = joinpath(ENV["HOME"],"data","kayaks","07_18")
    trialstr = "07_18_parsed_set1";
    datadir = joinpath(topdir,trialstr);
    ijldname = datadir * "/exp1" * ".jld"
elseif expnum==2
    topdir = joinpath(ENV["HOME"],"data","kayaks","07_18")
    trialstr = "07_18_parsed_set2";
    datadir = joinpath(topdir,trialstr);
    ijldname = datadir * "/exp1" * ".jld"
elseif expnum==3
    topdir = joinpath(ENV["HOME"],"data","kayaks","07_20")
    trialstr = "07_20_parsed_set1";
    datadir = joinpath(topdir,trialstr);
    ijldname = datadir * "/exp1" * ".jld"
end

cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml");
chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");
cfgd=loadConfigFile(cfgFile)
cfgd["range_snr_floor"]=0.95
iload = load(ijldname)
nrange = iload["nrange"]
irange = iload["irange"];
itemp = readdlm(datadir*"/inav.csv",',',Float64,'\n', skipstart=1);

# total length defined here
saswindow = 6;
nsasfac = 15;

wstart = 370; wlen = saswindow*nsasfac;
dataframes = collect(nrange[wstart]:nrange[wstart]+wlen);

sasstart = 1;
sasdataframes = [zeros(Int,saswindow) for _ in 1:nsasfac];
sasposes = [];
irangenew = irange[wstart] : 5 : irange[wstart+nsasfac*saswindow]
irangewindowed = collect(irangenew)[1:saswindow:nsasfac*saswindow+1]
igt = [itemp[irangewindowed,2] itemp[irangewindowed,4]];

poses = [Symbol("x$i") for i in 1:wlen]
posData = importdata_nav(dataframes,datadir=datadir);
navchecked, errorind = sanitycheck_nav(posData)

# dposData = deepcopy(posData)
# cumulativeDrift!(dposData,[0.0;0],[0.1,0.1])

# Check the trajectory
posPl = Gadfly.plot(layer(x=posData[:,1],y=posData[:,2], Geom.point, Geom.path, Theme(default_color=colorant"green",highlight_width = 0pt)),layer(x=igt[:,1],y=igt[:,2], Geom.point, Geom.path, Theme(default_color=colorant"red" ,highlight_width = 0pt))); posPl |> PDF("/tmp/test.pdf");

for i=1:nsasfac
    framestart = nrange[wstart]+sasstart-1+(i-1)*saswindow;
    sasdataframes[i] = collect(framestart:framestart+saswindow-1)
    push!(sasposes,[Symbol("x$j") for j in sasstart+(i-1)*saswindow:sasstart+saswindow-1+(i-1)*saswindow])
end

fg = initfg();
rtkCov= Matrix(Diagonal([0.1;0.1].^2));
for i in 1:wlen
  addVariable!(fg, poses[i], Point2)
  pp = PriorPoint2(MvNormal(posData[i,:], rtkCov))
  addFactor!(fg, [poses[i];], pp, autoinit=true)
  tmpInit = approxConv(fg,Symbol("x$(i)f1"),poses[i]);
  XXkde = manikde!(tmpInit,getManifolds(fg,poses[i]));
  setValKDE!(fg,poses[i],XXkde);
end

beacons = [Symbol("l$j") for j in 1:nsasfac]

for i=1:nsasfac
    addVariable!(fg, beacons[i], Point2)

    waveformData = importdata_waveforms(sasdataframes[i],2, datadir=datadir);
    sas2d = prepareSAS2DFactor(saswindow, waveformData, rangemodel=:Correlator,
                               cfgd=cfgd, chirpFile=chirpFile)
    addFactor!(fg, [beacons[i];sasposes[i]], sas2d, autoinit=false)
    @show currfac = ls(fg,beacons[i])[1]

    tmpInit = approxConv(fg,currfac,beacons[i]);
    XXkde = kde!(tmpInit);
    setValKDE!(fg,beacons[i],XXkde);

    if i>1
        dpμ = [0.0;0.0];
        dpσ = Matrix(Diagonal([2.0;2].^2))
        vp = Point2Point2(MvNormal(dpμ,dpσ))
        addFactor!(fg, [beacons[i-1];beacons[i]], vp, autoinit=false)
    end
end

writeGraphPdf(fg,viewerapp="", engine="neato", filepath = "/tmp/test.pdf")
# getSolverParams(fg).drawtree = true
getSolverParams(fg).limititers=500
#getSolverParams(fg).showtree = true
getSolverParams(fg).maxincidence = 200

tree = solveTree!(fg)

# Some Debuggging
# @time L1 = approxConv(fg,:l1x1x2x3x4x5x6x7x8f1,:l1)
# @time L1 = approxConv(fg,:l2x9x10x11x12x13x14x15x16f1,:l2)

mfIn = importdata_waveforms(sasdataframes[3],2, datadir=datadir);
range_mfData = zeros(Complex{Float64}, 8192, saswindow)
chirpIn = readdlm(chirpFile,',',Float64,'\n')
range_mf = prepMF(chirpIn,8192,saswindow) # MF
range_mf(mfIn,range_mfData) # matched filter
freq_sampling= 37500.0; sound_speed= 1481.0 ;
ranget = [0:8000-500;]*sound_speed/freq_sampling
bss = AliasingScalarSampler(ranget,exp.(norm.(range_mfData[1:7501,3])),SNRfloor=0.95)
out = rand(bss,1000)
Gadfly.plot(layer(x=out, Geom.histogram)) |> PDF("/tmp/rangetest.pdf")

# plk = plotBeaconSolve(fg)
plk=[];
for sym in ls(fg) #plotting all landmarks labeled
    if occursin(r"l",string(sym))
        X1 = getKDEMax(getBelief(fg,sym))
        push!(plk, layer(x=[X1[1];],y=[X1[2];],label=String["$sym";], Geom.point,Geom.label), Theme(default_color=colorant"red",point_size = 1.5pt,highlight_width = 0pt))
        L1pd = predictbelief(fg,sym,ls(fg, sym))
        push!(plk,layer(x=L1pd[1,:],y=L1pd[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))
    end
end

# K1 = plotKDEContour(kde!(getVal(fg,:l1)[1:2,:]),xlbl="X (m)", ylbl="Y (m)",levels=5,layers=true);
# push!(plk,K1...)
# push!(plk,Gadfly.Theme(key_position = :none));

push!(plk,layer(x=posData[:,1],y=posData[:,2], Geom.point, Geom.path, Theme(default_color=colorant"blue")), layer(x=igt[:,1],y=igt[:,2], Geom.point, Geom.path, Theme(default_color=colorant"green" )))
# push!(plk, Coord.cartesian(xmin=0, xmax=120, ymin=-80, ymax=20,fixed=true))
plkplot = Gadfly.plot(plk...); plkplot |> PDF("/tmp/test.pdf")

plk= [];

push!(plk,layer(x=posData[:,1],y=posData[:,2], Geom.path, Theme(default_color=colorant"green")), layer(x=igt[:,1],y=igt[:,2], Geom.point, Geom.path, Theme(default_color=colorant"red" )))

L1 = getVal(getVariable(fg, :l1))
push!(plk,layer(x=L1[1,:],y=L1[2,:],Geom.histogram2d(xbincount=300, ybincount=300)))

# push!(plk, Coord.cartesian(xmin=40, xmax=100, ymin=-120, ymax=-50,fixed=true))
plkplot = Gadfly.plot(plk...); plkplot |> PDF("/tmp/test.pdf");

## Plot the SAS factor pairs
fsym = :l1x1x2x3x4x5x6x7x8f1
pl = plotSASPair(fg, fsym, show=true, filepath="/tmp/testDP_2pp_3vpvp_init.pdf");


text:
function RungeKuttaExplicit!(V,dt,Fcn,CG,Param)
RK=Param.RK;
fV=Param.fV
Vn=Param.Vn

Vn .= V
for iStage=1:RK.nStage
  V .= Vn;
  for jStage=1:iStage-1
    @views V .= V .+ dt*RK.ARKE[iStage,jStage] .* fV[:,:,:,jStage];
  end
  Fcn(view(fV,:,:,:,iStage),V,CG,Param);
end
V .= Vn;
for iStage=1:RK.nStage
  @views V .= V .+ dt*RK.bRKE[iStage] .* fV[:,:,:,iStage];
end
end


text:
export
  local_lind,
  evolve_generator

"""
    local_lind(A[; epsilon])

Split the elements of the matrix `A` into a collection of sparse matrices with
exactly one non-zero element. Martices are created if the absolute value of the
nonzero element is there are not smaller than `epsilon`, where `epsilon` should
be nonnegative. The `epsilon` defaults to `eps()` if not specified.

# Examples
```jldoctest; setup = :(using QSWalk)
julia> A = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> local_lind(A)
4-element Array{SparseArrays.SparseMatrixCSC{Float64,Ti} where Ti<:Integer,1}:

  [1, 1]  =  1.0

  [1, 2]  =  2.0

  [2, 1]  =  3.0

  [2, 2]  =  4.0

julia> local_lind(A, epsilon = 1.5)
3-element Array{SparseArrays.SparseMatrixCSC{Float64,Ti} where Ti<:Integer,1}:

  [1, 2]  =  2.0

  [2, 1]  =  3.0

  [2, 2]  =  4.0

```
"""
function local_lind(A::Matrix{T}; epsilon::Real = eps()) where T<:Number
  @argumentcheck epsilon >= 0 "Epsilon should be nonegative"

  L = SparseMatrixCSC{T}[]
  for i = 1:size(A, 1), j = 1:size(A, 2)
    if abs(A[i, j]) >=  epsilon
        push!(L, A[i, j]*ketbra(i, j, size(A, 1)))
    end
  end
  L
end

function local_lind(A::SparseMatrixCSC{T}; epsilon::Real = eps()) where T<:Number
  @argumentcheck epsilon>= 0 "Epsilon should be nonegative"

  L = SparseMatrixCSC{T}[]
  for i = 1:size(A, 1), j = A[i, :].nzind
    if abs(A[i, j]) >=  epsilon
        push!(L, A[i, j]*ketbra(i, j, size(A, 1)))
    end
  end
  L
end

"""

    evolve_generator_create(H, L, localH, α, β)

Internal function for creating the generator for the evolution superoperator.
Given Hamiltonian `H`, collection of Lindblad operator `L`, local Hamiltonian
`localH` and scaling parameters `α` and `β` the function computes

``-i α (H ⊗ 1 - 1 ⊗ H) + β (-i(localH ⊗ 1-1 ⊗ localH)+∑ (L ⊗ L̄ - 1/2(L^†L ⊗ 1 + 1 ⊗ L^T L̄ )))``

where α and β should sum to one, or both be equal to one. In the later case
there are ignored.

*Note:* The function does not check whether `H` or `localH` are hermitian.

# Arguments
- `H`: Hamiltonian, must be hermitian,
- `L`: collection of Lindblad operators, each must be of the same size as `H`,
- `locH`: local Hamiltonian, suggested for nonmoralized QS walk, must be
hermitian and of the size of `H`,
- `α`: scaling parameter corresponding to Hamiltonian, should be in [0, 1],
- `β`: scaling parameter corresponding to Lindbladian part, should be in [0, 1].

# Return
The function return the generator, which can be used in `evolve` function.
"""
function evolve_generator_create(H::AbstractMatrix{<:Number},
                                 L::AbstractVector{<:AbstractMatrix{<:Number}},
                                 localH::AbstractMatrix{<:Number},
                                 α::Real,
                                 β::Real)
  @argumentcheck size(H) !=  (0, 0) "Matrix H must not be sizeless"
  @argumentcheck size(H, 1) == size(H, 2) "Matrix H must be square"
  @assert all([size(lindbladian) == size(H) for lindbladian in L]) "Lindblad operators must be of the same size as Hamiltonian"
  @assert size(H) == size(localH) "Matrix localH must be of the same size as H"
  @argumentcheck 0 <=  α <=  1 && 0 <=  β <=  1 "Value of ω must be nonngeative and smaller than one"

  F = spzeros(ComplexF64, (size(H).^2)...)
  id = Matrix{ComplexF64}(I, size(H)...)
  for i = 1:length(L)
      F +=  kron(L[i], conj(L[i]))-0.5*kron(L[i]'*L[i], id)-0.5*kron(id, transpose(L[i])*conj(L[i]))
  end
  F +=  im*(kron(id, conj(localH))-kron(localH, id))
  F = α*im*(kron(id, conj(H))-kron(H, id)) + β*F
  F
end

"""
    evolve_generator(H, L[, localH][, ω])

Create the generator for the evolution superoperator. Given Hamiltonian `H`,
collection of Lindblad operator `L`, local Hamiltonian `localH` and scaling
parameter `ω`, the generator is obtained as a sum

``-i(1-ω) (H ⊗ 1 - 1 ⊗ H) + ω (-i(localH ⊗ 1 - 1 ⊗ localH) + ∑(L ⊗ L̄ - 1/2(L^†L ⊗ 1 + 1 ⊗ L^T L̄ )))``

The last two arguments are optional.

If `localH` is not given, it defaults to sparse zero matrix of the size of `H`.

If `ω` is not given, both parts are taken with the same intensity and the global
operator takes the form

``-i(H ⊗ 1 - 1 ⊗ H) + (-i(localH ⊗ 1 - 1 ⊗ localH) + ∑(L ⊗ L̄ - 1/2(L^†L ⊗ 1 + 1 ⊗ L^T L̄ )))``

# Arguments
- `H`: Hamiltonian, must be hermitian,
- `L`: collection of Lindblad operators, each must be of the same size as `H`,
- `localH`: local Hamiltonian, suggested for nonmoralized QS walk, must be hermitian and of the size of `H`,
- `ω`: scaling parameter, should be in [0, 1].

# Return
The generator matrix, which can be used in `evolve` function.

# Examples

```jldoctest; setup = :(using QSWalk)
julia> H, L, localH = [0 1+im; 1-im 0], [0. 1; 0 0], [1. 0.; 0. 1.]
(Complex{Int64}[0 + 0im 1 + 1im; 1 - 1im 0 + 0im], [0.0 1.0; 0.0 0.0], [1.0 0.0; 0.0 1.0])

julia> evolve_generator(H, [L], localH, 1/2)
4×4 Array{Complex{Float64},2}:
  0.0+0.0im    0.5+0.5im    0.5-0.5im   0.5+0.0im
 -0.5+0.5im  -0.25+0.0im    0.0+0.0im   0.5-0.5im
 -0.5-0.5im    0.0+0.0im  -0.25+0.0im   0.5+0.5im
  0.0+0.0im   -0.5-0.5im   -0.5+0.5im  -0.5+0.0im

```
"""
function evolve_generator(H::AbstractMatrix{<:Number},
                          L::AbstractVector{<:AbstractMatrix{<:Number}},
                          localH::AbstractMatrix{<:Number},
                          ω::Real)
  evolve_generator_create(H, L, localH, 1-ω, ω)
end,

function evolve_generator(H::AbstractMatrix{<:Number},
                          L::AbstractVector{<:AbstractMatrix{<:Number}},
                          localH::AbstractMatrix{<:Number} = spzeros(eltype(H), size(H)...))
  evolve_generator_create(H, L, localH, 1., 1.)
end,

function evolve_generator(H::AbstractMatrix{T},
                          L::AbstractVector{<:AbstractMatrix{<:Number}},
                          ω::Real) where T<:Number
  evolve_generator_create(H, L, spzeros(T, size(H)...), 1-ω , ω)
end