text:
"""
    Module Geometry
All mesh information and setup related to element connectivity
"""
module Geometry

import LinearAlgebra.norm
import Basis1D
import Basis2D.dPhi2D

export Mesh
export periodic_square
export setup_nodes!, setup_quads!
export precompute_jacobians!, precompute_ldg_switch!
export output_mesh

" Dimension of overall mesh "
const DIM = 2
" Number of faces per element "
const N_FACES = 2*DIM
" Number of vertices per element "
const N_VERTICES = 2^DIM
" Number of vertices per element face "
const N_FVERTICES = 2^(DIM-1)

"""
    struct Mesh
Contains all mesh information related to domain and element connectivity
"""
mutable struct Mesh
    " Global number of elements in mesh "
    n_elements::Int
    " Global number of vertices in mesh "
    n_vertices::Int
    " Globally minimum edge length across all elements "
    min_dx::Float64
    " Global vector of N-D vertices that form corners of elements "
    vertices
    """ Element-to-Vertex mapping:
    For every element iK, vertex i, `e2v[i,iK]` = vertex index of local vertex i """
    e2v
    """ Face-to-Vertex mapping:
    For every face iF, face vertex i, `f2v[i,iF]` = this element's vertex index """
    f2v
    """ Element-to-Element neighbor mapping:
    For every element iK, face iF, `e2e[iF,iK]` = element index of neighbor on face iF """
    e2e
    """ Element-to-Face mapping:
    For every element iK, face iF, `e2f[iF,iK]` = face index on neighbor's face array """
    e2f
    """ Normals mapping:
    For every element iK, face iF, `normals[:,iF,iK]` = outward normal vector of face iF """
    normals
    """
    Bilinear mapping:

    For every element iK, dimension l,
    ``x_l`` = `sum(bilinear_mapping[:,l,iK]*phi[xi_l,:])`
    """
    bilinear_mapping

    ######################################
    # Initialized after solver is created
    ######################################
    " Order of DG method "
    order::Int
    " Local number of DG nodes per element "
    n_nodes::Int
    " Local number of quadrature points per element "
    n_quads::Int
    " Local number of DG nodes per element face "
    n_face_nodes::Int
    " Local number of quadrature points per element face "
    n_face_quads::Int
    """
    True coordinates of DG nodes for every element

    For every element iK, node i, `coord_nodes[:,i,iK]` = ND coordinate of node i
    """
    coord_nodes
    """
    True coordinates of DG quadrature points for every element

    For every element iK, quad point i, `coord_quads[:,i,iK]` = ND
    coordinate of quad point i
    """
    coord_quads
    """
    Jacobian of transformation evaluated at each quadrature point, element

    For every element iK, quad point iQ, `Jk[:,:,iQ,iK]` ∈ ℜ^(DIM x DIM)
    is the Jacobian of the mapping from the reference element to element iK at
    that quadrature point.
    """
    Jk
    """
    |det Jacobian| of transformation evaluated each quadrature point, element

    For every element iK, quad point iQ, `Jkdet[iQ,iK]` is |det Jac| of the
    mapping from the reference element to element iK at quad point iQ.
    """
    Jkdet
    """
    Jacobian of transformation evaluated at each face quadrature point, element

    For every element iK, face iF, face quad point iQ, `JkF[:,:,iQ,iF,iK]`
    ∈ ℜ^(DIM x DIM) is the Jacobian of the mapping from the reference face
    element to face iF of element iK at face quad point iQ
    """
    JkF
    """
    |det Jacobian| of face integral evaluated at each face quadrature point, element

    For every element iK, face iF, face quad point iQ, `JkFdet[iQ,iF,iK]` is
    |gamma'(t)| from reference face element to face iF at face quad point iQ.
    """
    JkFdet
    """
    Element face to node mapping

    For every face iF (same mapping across all elements), face node iFN,
    `ef2n[iFN,iF]` = local node ID of face node iFN on face iF
    """
    ef2n
    """
    Neighbor face node to face node mapping

    For every face node iFN (same mapping across all faces/elements),
    `nfn2fn[iFN]` = face node ID on neighbor's face

    Neighbor node: `solution[nfn2fn[iFN], e2f[iF,iK], :, e2e[iF,iK]]`
    """
    nfn2fn
    """
    Element face to quadrature point mapping

    For every face iF (same mapping across all elements), face quad iFQ,
    `ef2q[iFQ,iF]` = local quad ID of face quad iFQ on face iF
    """
    ef2q
    """
    Neighbor face quadrature point to quadrature point mapping

    For every face quadrature point ID iFQ (same mapping across all faces/elements),
    `nfq2fq[iFQ]` = face quadrature point ID on neighbor's face

    Neighbor quad: `solution[nfq2fq[iFQ], e2f[iF,iK], :, e2e[iF,iK]]`
    """
    nfq2fq
    """
    Local DG switch (only needed for 2nd-order applications):

    For every element iK, face iF, `ldg_switch[iF,iK]` = true if ldg flux
    should be taken from this element's face values
    """
    ldg_switch

    function Mesh(ne_, nv_, mdx_, v_, e2v_, f2v_, e2e_, e2f_, normals_, bm_)
        return new(ne_, nv_, mdx_, v_, e2v_, f2v_, e2e_, e2f_, normals_, bm_,
            0,0,0,0,0,[],[],[],[],[],[],[],[],[],[],[])
    end
end

"""
    Base.show(io::IO, mesh::Mesh)
Allow for pretty printing of mesh type
"""
function Base.show(io::IO, mesh::Mesh)
    print(io, "$(mesh.n_vertices) vertices connected with $(mesh.n_elements) elements.")
end

"""
    periodic_square(ns::Array)
Create a periodic square [0,1]^2 with prod(ns) global number of elements
"""
function periodic_square(ns)::Mesh
    @assert length(ns) == DIM
    top_right = ones(DIM)
    bot_left  = zeros(DIM)
    n_elements = prod(ns)
    n_vertices = prod(ns.+1)
    min_edge_length = minimum((top_right .- bot_left)./ns)

    # Initialize vertices
    nd_vertices = Array{Float64,DIM+1}(undef,DIM,(ns.+1)...)
    for iy = 1:ns[2]+1
        curry = (iy-1)*(top_right[2]-bot_left[2])/ns[2]+bot_left[2]
        for ix = 1:ns[1]+1
            currx = (ix-1)*(top_right[1]-bot_left[1])/ns[1]+bot_left[1]
            nd_vertices[1,ix,iy] = currx
            nd_vertices[2,ix,iy] = curry
        end
    end
    vertices = reshape(nd_vertices, DIM, n_vertices)

    # Initialize element to vertices array
    e2v = Array{Int,2}(undef, N_VERTICES, n_elements)
    for iy = 1:ns[2]
        yoff1 = (iy-1)*(ns[1]+1)
        yoff2 = (iy  )*(ns[1]+1)
        for ix = 1:ns[1]
            xoff1 = ix # Only fastest direction contains the 1-indexing
            xoff2 = ix+1
            e_index = ix+(iy-1)*ns[1]
            # Stores vertices in CCW order
            e2v[1, e_index] = xoff1+yoff2
            e2v[2, e_index] = xoff1+yoff1
            e2v[3, e_index] = xoff2+yoff1
            e2v[4, e_index] = xoff2+yoff2
        end
    end

    # Initialize element to face arrays
    e2e = Array{Int, 2}(undef, N_FACES, n_elements)
    e2f = Array{Int, 2}(undef, N_FACES, n_elements)
    for iy = 1:ns[2]
        iyM = mod(iy-2,ns[2]) * ns[1]
        iyP = mod(iy,  ns[2]) * ns[1]
        iy0 = (iy-1) * ns[1]
        for ix = 1:ns[1]
            # Only fastest direction contains the 1-indexing
            ixM = mod(ix-2,ns[1])+1
            ixP = mod(ix,  ns[1])+1
            ix0 = ix
            e_index = ix0+iy0

            # Neighbor elements are stored in -x,-y,+x,+y order to be CCW
            e2e[1,e_index] = ixM+iy0
            e2e[2,e_index] = ix0+iyM
            e2e[3,e_index] = ixP+iy0
            e2e[4,e_index] = ix0+iyP

            # Face ID of this element's -x face will be neighbor's +x face
            e2f[1,e_index] = 3
            e2f[2,e_index] = 4
            e2f[3,e_index] = 1
            e2f[4,e_index] = 2
        end
    end

    # Initialize f2v
    f2v = init_f2v()

    # Initialize bilinear_mapping
    bilinear_mapping = init_bilinear_mapping(vertices,e2v)

    # Initialize normals and min_dx
    (normals, min_dx) = init_normals(vertices,f2v,e2v)

    mesh = Mesh(n_elements, n_vertices, min_dx, vertices, e2v, f2v, e2e, e2f, normals, bilinear_mapping)
    return mesh
end

"""
    init_f2v()
Label faces of a rectangle according to CCW vertices
"""
function init_f2v()
    f2v = Array{Int,2}(undef, N_FVERTICES, N_FACES)
    # Mapping assumes input vertices are CCW
    for iF = 1:N_FACES
        f2v[1,iF] = iF
        f2v[2,iF] = mod(iF, N_FACES)+1
    end
    return f2v
end

"""
    init_bilinear_mapping(vertices,e2v)
Compute the bilinear mapping from element coordinates
"""
function init_bilinear_mapping(vertices, e2v)
    map_order = 1
    map_dofs = (map_order+1)^DIM
    # Sanity check that bilinear mapping is the correct size
    @assert map_dofs == N_VERTICES

    # Gmsh =   4---3   Mapping =  3---4
    # CCW      |   |   wants      |   |
    # corners  1---2              1---2
    c2m = [1,2,4,3]
    bilinear_mapping = Array{Float64,3}(undef, N_VERTICES, DIM, size(e2v,2))
    for iK = 1:size(e2v,2)
        for l = 1:DIM
            for iV = 1:N_VERTICES
                bilinear_mapping[c2m[iV],l,iK] = vertices[l,e2v[iV,iK]]
            end
        end
    end
    return bilinear_mapping
end

"""
    init_normals(vertices, f2v, e2v)
Return the normals from the element coordinates, assuming 2D mesh is ordered CCW.
Simultaneously computes the minimum edge length across global mesh.
"""
function init_normals(vertices, f2v, e2v)
    n_elements = size(e2v,2)
    normals = zeros(Float64, DIM, N_FACES, n_elements)
    min_dx = typemax(Float64)
    for iK = 1:n_elements
        for iF = 1:N_FACES
            a = f2v[1,iF]
            b = f2v[2,iF]
            AB = vertices[:,e2v[b,iK]].-vertices[:,e2v[a,iK]]
            len = norm(AB,2)
            # Normal to (dx,dy): (dy,-dx)
            normals[1,iF,iK] = AB[2]/len
            normals[2,iF,iK] = -AB[1]/len

            min_dx = min(min_dx, len) # Compute minimum edge length
        end
    end
    return (normals, min_dx)
end

"""
    setup_nodes!(mesh::Mesh, InterpTk, order; compute_global=false)
Initialize element nodes within Mesh. If `compute_global`, then also compute
the global nodes from bilinear mapping of reference nodes and store in
`coord_nodes`. This would be needed if solutions, for instance, to be stored at
nodes would be needed to view.
"""
function setup_nodes!(mesh::Mesh, InterpTk, order; compute_global=false)
    mesh.order = order
    mesh.n_nodes = size(InterpTk,1)
    if compute_global
        mesh.coord_nodes = Array{Float64,3}(undef, DIM, mesh.n_nodes, mesh.n_elements)
        for iK = 1:mesh.n_elements
            for l = 1:DIM
                mesh.coord_nodes[l,:,iK] = InterpTk*mesh.bilinear_mapping[:,l,iK]
            end
        end
    end

    (mesh.n_face_nodes,mesh.ef2n,mesh.nfn2fn) = init_face_map(mesh, mesh.order+1)
end

"""
    setup_quads!(mesh::Mesh, InterpTkQ; compute_global=false)
Initialize element quadrature points within Mesh. If `compute_global`, then also
compute global quadrature points from bilinear mapping of reference quad points
and store in `coord_quads`. This would be needed if material properties, for instance,
to be stored at quadrature points would be needed.
"""
function setup_quads!(mesh::Mesh, InterpTkQ; compute_global=false)
    if mesh.order <= 0
        error("Mesh must setup nodes before quads!")
    end
    mesh.n_quads = size(InterpTkQ,1)
    if compute_global
        mesh.coord_quads = Array{Float64,3}(undef, DIM, mesh.n_quads, mesh.n_elements)
        for iK = 1:mesh.n_elements
            for l = 1:DIM
                mesh.coord_quads[l,:,iK] = InterpTkQ*mesh.bilinear_mapping[:,l,iK]
            end
        end
    end

    nQ = ceil(Int, mesh.order+1/2) # Quadrature points per face
    (mesh.n_face_quads, mesh.ef2q, mesh.nfq2fq) = init_face_map(mesh, nQ)
end

"""
    init_face_map(mesh::Mesh, size1D)
Initialize face maps for a given number of nodes.
Assumes that each face is a side with exactly size1D nodes.
For every element iK, face iF, face node iFN:

my node       @: `soln[efmap[iFN, iF], :, iK]`

neighbor node @: `soln[efmap[nfmap[iFN], e2f[iF,iK]], :, e2e[iF,iK]]`
"""
function init_face_map(mesh::Mesh, size1D)
    efmap = Array{Int64,2}(undef, size1D, N_FACES)
    # efmap has explicit dependence on f2v. Only handles CCW vertices
    # 4---3
    # |   |
    # 1---2
    for iF = 1:N_FACES
        if mesh.f2v[:,iF] == [4,1]
            # -x direction face
            x_off = 1 # 1-indexed
            y_off = (size1D-1:-1:0).*size1D
            efmap[:,iF] = x_off.+y_off
        elseif mesh.f2v[:,iF] == [1,2]
            # -y direction face
            y_off = 0
            x_off = 1:size1D # 1-indexed
            efmap[:,iF] = x_off.+y_off
        elseif mesh.f2v[:,iF] == [2,3]
            # +x direction face
            x_off = size1D # 1-indexed
            y_off = (0:size1D-1).*size1D
            efmap[:,iF] = x_off.+y_off
        elseif mesh.f2v[:,iF] == [3,4]
            # +y direction face
            y_off = (size1D-1)*size1D
            x_off = size1D:-1:1 # 1-indexed
            efmap[:,iF] = x_off.+y_off
        else
            error("Vertices for f2v are not CCW!")
        end
    end

    # nfmap is reversed version of efmap
    nfmap = collect(size1D:-1:1)

    # Face of 2D mesh is 1D
    fnodes = size1D
    return (fnodes,efmap,nfmap)
end

"""
    precompute_jacobians!(mesh::Mesh, xTk1D, xQ1D)
Precompute Jacobian of mappings from reference element to all elements,
evaluated at quadrature points.

Depending on the specified flags, will
initialize mesh.Jk (Jacobians), mesh.Jkdet (|det Jac|), mesh.JkF (face Jacobians),
and/or mesh.JkFdet(|det face Jac|).

xTk1D is 1D points that bilinear_mapping Tk is defined on
xQ1D is 1D quadrature points, and assumes faces are xQ1D, while volume is (xQ1D,xQ1D)
"""
function precompute_jacobians!(mesh::Mesh, xTk1D, xQ1D)
    # Compute grad of mapping from reference bases, evaluated on element quad points
    xyTk = (xTk1D,xTk1D)
    dPhiTkQ = dPhi2D(xyTk, (xQ1D,xQ1D))

    # Compute grad of mapping from reference bases, evaluated on face quad points
    n_from = prod(size.(xyTk,1))
    n_to = size(xQ1D,1) # 1D face
    dPhiTkQF = Array{Float64,4}(undef, n_to,n_from,DIM,N_FACES)
    for iF = 1:N_FACES
        # xyQF = "kronecker product" that represents face quad points on face iF
        if mesh.f2v[:,iF] == [4,1]
            # -x direction face
            xQF = [xTk1D[1]]
            yQF = reverse(xQ1D)
        elseif mesh.f2v[:,iF] == [1,2]
            # -y direction face
            xQF = xQ1D
            yQF = [xTk1D[1]]
        elseif mesh.f2v[:,iF] == [2,3]
            # +x direction face
            xQF = [xTk1D[2]]
            yQF = xQ1D
        elseif mesh.f2v[:,iF] == [3,4]
            # +y direction face
            xQF = reverse(xQ1D)
            yQF = [xTk1D[2]]
        else
            error("Vertices for f2v are not CCW!")
        end
        xyQF = (xQF,yQF)
        dPhiTkQF[:,:,:,iF] = dPhi2D(xyTk, xyQF)
    end

    # Jacobian = dPhiTkQ[:,:,lj]*bilinear_mapping[:,li,iK]
    mesh.Jk = Array{Float64,4}(undef, DIM,DIM,mesh.n_quads,mesh.n_elements)
    mesh.Jkdet = Array{Float64,2}(undef, mesh.n_quads,mesh.n_elements)
    for iK = 1:mesh.n_elements
        for li = 1:DIM
            for lj = 1:DIM
                mesh.Jk[li,lj,:,iK] = dPhiTkQ[:,:,lj]*mesh.bilinear_mapping[:,li,iK]
            end
        end
        Jac = @view mesh.Jk[:,:,:,iK]
        mesh.Jkdet[:,iK] = Jac[1,1,:].*Jac[2,2,:].-Jac[2,1,:].*Jac[1,2,:]
    end
    # Face Jacobian = dPhiTkQF[:,:,lj,iF]*bilinear_mapping[:,li,iK]
    mesh.JkF = Array{Float64,5}(undef, DIM,DIM,mesh.n_face_quads,N_FACES,mesh.n_elements)
    mesh.JkFdet = Array{Float64,3}(undef, mesh.n_face_quads,N_FACES,mesh.n_elements)
    for iK = 1:mesh.n_elements
        for iF = 1:N_FACES
            for li = 1:DIM
                for lj = 1:DIM
                    mesh.JkF[li,lj,:,iF,iK] = dPhiTkQF[:,:,lj,iF]*mesh.bilinear_mapping[:,li,iK]
                end
            end
            # 1D integrals => only needs |gamma'(t)| term along correct x_i_{l_i}
            lj = 0
            if mesh.f2v[:,iF] == [4,1]
                # -x direction face
                lj = 2
            elseif mesh.f2v[:,iF] == [1,2]
                # -y direction face
                lj = 1
            elseif mesh.f2v[:,iF] == [2,3]
                # +x direction face
                lj = 2
            elseif mesh.f2v[:,iF] == [3,4]
                # +y direction face
                lj = 1
            else
                error("Vertices for f2v are not CCW!")
            end
            for iQ = 1:mesh.n_face_quads
                mesh.JkFdet[iQ,iF,iK] = norm(mesh.JkF[:,lj,iQ,iF,iK],2)
            end
        end
    end
end

"""
    function precompute_ldg_switch!(mesh::Mesh)
Precompute the Local DG switch for 2nd-order applications. Opposite faces of
each element are guaranteed to have opposing switch values.
"""
function precompute_ldg_switch!(mesh::Mesh)
    # Initialize opposite_face array
    opposite_face = zeros(Int64,Geometry.N_FACES)
    for iF = 1:Geometry.N_FACES
        for iF2 = 1:Geometry.N_FACES
            # Opposite faces share no common vertices
            if isempty(intersect(mesh.f2v[:,iF],mesh.f2v[:,iF2]))
                opposite_face[iF] = iF2
                break
            end
        end
    end

    visited = falses(Geometry.N_FACES, mesh.n_elements)
    mesh.ldg_switch = zeros(Bool, Geometry.N_FACES, mesh.n_elements)
    while true
        index = findfirst(x -> !x, visited)
        if isnothing(index)
            break
        end
        iF = index[1]; iK = index[2];
        init_iF = iF; init_iK = iK
        choose_this_element = true
        while true
            mesh.ldg_switch[iF,iK] = choose_this_element
            visited[iF,iK] = true
            nK = mesh.e2e[iF,iK]
            if nK < 0 # Boundary
                if choose_this_element # Now reverse direction from initial choice
                    choose_this_element = !choose_this_element
                    iK = init_iK
                    iF = opposite_face[init_iF]
                else # We have already reversed once
                    break
                end
            else # Inter-element border
                nF = mesh.e2f[iF,iK]
                mesh.ldg_switch[nF,nK] = !choose_this_element
                visited[nF,nK] = true
                iK = nK
                iF = opposite_face[nF]
                if visited[iF,iK]
                    if choose_this_element # Now reverse direction from initial choice
                        choose_this_element = !choose_this_element
                        iK = init_iK
                        iF = opposite_face[init_iF]
                    else # We have already reversed once
                        break
                    end
                end
            end
        end
    end
end

"""
    output_mesh(mesh::Mesh, filename)
Outputs the given mesh as a Gmsh 4.1 file format .msh file
"""
function output_mesh(mesh::Mesh, filename)
    file = open(filename, "w")
    # Uses 4.1 MSH file format from https://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
    println(file, "\$MeshFormat")
    println(file, "4.1 0 8") # version 4.1, ASCII
    println(file, "\$EndMeshFormat")

    # Vertices (mesh.vertices)
    println(file, "\$Nodes")
    #     num_entities num_vertices start_vertex end_vertex
    println(file, "1 $(mesh.n_vertices) 1 $(mesh.n_vertices)")
    #     2D entity (surface) #1, no parametric coords, num_vertices in entity
    println(file, "2 1 0 $(mesh.n_vertices)")
    for iV = 1:mesh.n_vertices
        println(file, "$iV")
    end
    for iV = 1:mesh.n_vertices
        # Vertices must be in 3D
        println(file, "$(mesh.vertices[1,iV]) $(mesh.vertices[2,iV]) 0.0")
    end
    println(file, "\$EndNodes")

    # Elements (mesh.e2v)
    println(file, "\$Elements")
    #     num_entities num_elements start_element end_element
    println(file, "1 $(mesh.n_elements) 1 $(mesh.n_elements)")
    #     2D entity (surface) #1, element type 3 (4-coord quad), num_elements
    println(file, "2 1 3 $(mesh.n_elements)")
    for iK = 1:size(mesh.e2v,2)
        print(file, "$iK")
        for iV = 1:size(mesh.e2v,1)
            print(file, " $(mesh.e2v[iV,iK])")
        end
        println(file)
    end
    println(file, "\$EndElements")
    # TODO: add periodicity if applicable

    close(file)
end

end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
#=
   From
   "First-Order Probabilistic Inference"
   https://www.youtube.com/watch?v=leIqVD4-Fks
   Time 45:06

   (The shown model in the talks is a BLOG model).

   numSoldierMean: 1095.9922

   Distributions of variable numBatalion (num:0)
   6.00000 =>    1058  (0.105800)
   2.00000 =>    1028  (0.102800)
   5.00000 =>    1021  (0.102100)
   1.00000 =>    1020  (0.102000)
   7.00000 =>    1019  (0.101900)
   10.00000 =>    1018  (0.101800)
   9.00000 =>     969  (0.096900)
   4.00000 =>     958  (0.095800)
   3.00000 =>     955  (0.095500)
   8.00000 =>     954  (0.095400)

   cf ~/blog/batallion.blog
      ~/webppl/batallion.wppl

=#

using Turing
include("jl_utils.jl")

@model function batalion()
    numBatalion ~ DiscreteUniform(1,10)
    large ~ filldist(flip(0.6),numBatalion)
    region ~ filldist(Poisson(numBatalion),numBatalion)
    numSoldier = tzeros(numBatalion)
    for bat in 1:numBatalion
        if large[bat]
            numSoldier[bat] ~ Poisson(1500)
        else
            if region[bat] == 2
                numSoldier[bat] ~ Uniform(300,301) # Dirac(300)
            else
                numSoldier[bat] ~ Poisson(500)
            end
        end
    end
    numSoldierMean ~ Dirac(sum(numSoldier) / numBatalion)
end

model = batalion()

num_chains = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 10_000)
chs = sample(model, IS(), 10_000)
#chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chs)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

display(chs)
# display(plot(chs))

show_var_dist_pct(chs,:numBatalion)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
include("./DOF.jl")
include("./Family.jl")
include("./MutableMaxHeap.jl")

# Revisit later
# function r_maximin(dofs; r=1.1, n_bins=100, ℓ0=fill(Inf, length(dofs)))
#     # Initializing vector containing the reordering
#     P = Vector{Int}(undef, 0)
#     # Checklist to see which element has been added to orderigng already
#     checklist = trues(length(dofs))
#     # choosing the first pivot as having minimal distance according to the "hard coded" distance
#     push!(P, argmax(ℓ0))
#     # vector that stores the length scales of the dofs
#     ℓ_out = similar(ℓ0)
#     ℓ_out[P[1]] = ℓ0[P[1]]
#     checklist[P[1]] = false
#     # updating hte distances according to first pivot
#     ℓ = min.(distance.(dofs, Ref(dofs[P[1]])), ℓ0)
#
#     # compute upper bound of diameter of domain
#     diam = 2 * maximum(ℓ)
#
#     f = Family(dofs)
#     # Setting the descendants of the first node to be all dofs, ordering the distance
#     descendants(f)[P[1]] = sort([(ℓ[k], k) for k = 1 : length(ℓ)])
#     parents(f) .= [(ℓ[k], P[1]) for k = 1 : length(ℓ)]
#
#     # instantiate rack
#     rack = Rack(n_bins, r, diam, ℓ)
#
#     # starting by looking in first bin
#     first_bin = 1
#     while ~iszero(first_bin)
#         # obtain new pivot
#         first_bin, pivot_id = next_dof(first_bin, rack)
#         # if pivot was not in last bin
#
#         if !iszero(pivot_id)
#             new_pivot!(P, ℓ_out, rack, f, pivot_id, dofs, ℓ, checklist)
#             checklist[pivot_id] = false
#         end
#     end
#     # Once all remaining dofs are in the smallest bin, we simply add them to the ordering in arbitrary order
#     # First set the length-scale of the remaining dofs to the minimal resolution
#     ℓ_out[setdiff(1 : length(dofs), P)] = ℓ[setdiff(1 : length(dofs), P)]
#     P = unique(vcat(P, 1 : length(dofs)))
#     return P, ℓ_out
# end

function new_pivot!(ℓ, h::MutableMaxHeap, f::Family, buffer, pivot, dofs)
    # add pivot to ordering. Note that this tags the pivot and leads to it NOT being added automatically to its descendants.
    pivot_id, pivot_ℓ = pivot.id, pivot.val
    new_column!(f, pivot_id)
    pivot_number = f.n_columns[1]
    pivot = dofs[pivot_id]

    # add the pivot to it's own list of children
    buffer_iterator = 1
    buffer[buffer_iterator] = Member(zero(typeof(pivot_ℓ)), pivot_id)
    # and downgrade its distance in the heap
    decrease_key!(h, pivot_id, 0.0)

    # load the distance from the pivot to its parent
    distance_pivot_parent = f.parents[pivot_id].val
    for index in column_iterator(f, f.parents[pivot_id].id)
        parent_descendant = f.rowval[index]
        distance_parent_descendant = parent_descendant.val

        # break loop if descendant is too far from parent to be descendant of pivot
        if distance_parent_descendant >  distance_pivot_parent + 2 * pivot_ℓ
            break
        # else compute distance and possibly add as descendant or parent
        elseif f.revP[parent_descendant.id] == zero(eltype(f.revP))
            distance_descendant_pivot = distance(dofs[parent_descendant.id], pivot)

            # If distance between pivot and descendant is small enough, add descendent to children of pivot
            # In this case, we also update the length scale of the pivot
            if distance_descendant_pivot < 2 * pivot_ℓ
                buffer_iterator += 1
                buffer[buffer_iterator] = Member(distance_descendant_pivot, parent_descendant.id)
                decrease_key!(h, parent_descendant.id, distance_descendant_pivot)
                ℓ[parent_descendant.id] = min(distance_descendant_pivot, ℓ[parent_descendant.id])
            end
            # if the descendant was not added to the ordering yet, update distance and check if pivot should become new parent of descendant.
            # Check if pivot should become the descendant's next parent
            if distance_descendant_pivot + 2 * ℓ[parent_descendant.id] < pivot_ℓ && distance_descendant_pivot < distance_parent_descendant
                f.parents[parent_descendant.id] = Member(distance_descendant_pivot, pivot_number)
            end
        end
    end
    # Sorts the descendants of the new pivot by distance to the pivot
    viewBuffer = view(buffer, 1 : buffer_iterator)
    sort!(viewBuffer, alg=QuickSort)
    new_children!(f, viewBuffer)
end

# computing the exact naive minimax ordering
function naive_maximin(dofs; ℓ0=fill(Inf, length(dofs)))
    P = Vector{Int}(undef, length(dofs))
    ℓ = Vector{Float64}(undef, length(dofs))
    P[1] = 1
    ℓ[1] = Inf
    for k = 2 : length(dofs)
        dists = fill(Inf, length(dofs))
        for (l, distval) in enumerate(dists)
            for m = 1 : (k - 1)
                dists[l] = min(dists[l], distance(dofs[P[m]], dofs[l]))
            end
        end
        P[k] = argmax(dists)
        ℓ[P[k]] = maximum(dists)
    end
    return P, ℓ
end

function maximin(dofs; ℓ0=fill(Inf, length(dofs)))
    N = length(dofs)

    # instantiate the family and choose first pivot
    f = Family(dofs)
    new_column!(f, argmax(ℓ0))
    # updating the distances according to first pivot and add corresponding descendants to f
    ℓ = min.(distance.(dofs, Ref(dofs[f.P[1]])), ℓ0)
    new_children!(f, sort([Member(ℓ[id], id) for id = 1 : length(dofs)]))
    # setting all parents to be given by the first node
    f.parents .= [Member(ℓ[id], 1) for id = 1 : N]

    # Instantiate buffer to be used later.
    buffer = Vector{eltype(f.rowval)}(undef, length(dofs))
    # Instantiate the mutable binary heap with all distances set to maximum
    h = MutableMaxHeap{Float64,Int,Int}([Node{Float64,Int,Int}(typemax(Float64), i, zero(Int)) for i = 1 : length(dofs)], 1 : N)
    # adapt keys in heap
    for k = 1 : N
        decrease_key!(h, k, ℓ[k])
    end
    # set the length scale of the first entry to the output value.
    ℓ[f.P[1]] = ℓ0[f.P[1]]

    # compute the remaining parts of the ordering.
    # Right now the ordering over all dofs is computed but in the future we could add additional criteria for the breaking the for loop, for instance when pivot_ℓ falls below a treshhold.
    for k = 2 : length(ℓ)
        # obtain new pivot
        pivot = top_node!(h)

        # if pivot was not in last bin
        new_pivot!(ℓ, h, f, buffer, pivot, dofs)
    end
    # First set the length-scale of the remaining dofs to the minimal resolution
    return f.P, ℓ
end


text:
#=

  Map coloring in Julia ConstraintSolver.jl

  Simple map coloring problem.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function map_coloring(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    n, _ = size(problem)
    num_colors = 4
    # println("n:$n num_colors:$num_colors")
    @variable(model, 1 <= x[1:n] <= num_colors, Int)

    for c1 in 1:n, c2 in 1:n
        if problem[c1,c2] == 1
            @constraint(model,x[c1] != x[c2])
        end
    end

    # symmertry breaking
    for c in 1:min(n,num_colors)
        @constraint(model,x[c] <= c)
    end

    # @objective(model,Max,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

map_coloring_problems = Dict(
    # Connections between these countries:
    # [belgium, denmark, france, germany, netherlands, luxembourg]
    :1 => resize_matrix([[0, 0, 1, 1, 1, 1],
                         [0, 0, 0, 1, 0, 0],
                         [1, 0, 0, 1, 1, 0],
                         [1, 1, 1, 0, 1, 1],
                         [1, 0, 1, 1, 0, 0],
                         [1, 0, 0, 1, 0, 0]])

)

@time map_coloring(map_coloring_problems[:1])

# random instance
function random_map(n=10,lim=0.5)
    m = zeros(Int8,n,n)
    for i in 1:n, j in 1:i-1
       if rand() >= lim
        m[i,j] = 1
        m[j,i] = 1
       end
    end
    return m
end

for n in 10:20
    println("\nn:$n")
    m = random_map(n,0.5)
    # display(m)
    @time map_coloring(m,true,false)
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using Dates
using JSON
using NMAP
using Test
using XMLDict

const dir = joinpath(dirname(pathof(NMAP)), "..", "test")
const resources = joinpath(dir, "resources")

@testset "Bulk Options" begin
    scanner = NMAP.Scanner()
    for name in names(NMAP, all=true, imported=false)
        f = getfield(NMAP, name)
        if f isa Function && NMAP.Option in Base.return_types(f)
            ms = methods(f).ms
            for method in ms
                test_scanner = NMAP.Scanner()
                if method.nargs < 2
                    opt = f()
                    test_scanner(opt)
                    len_diff = length(test_scanner.args)
                    cur_len  = length(scanner.args)
                    scanner(opt)
                    @test cur_len + len_diff == length(scanner.args)
                end
            end
        end
    end
end

@testset "XML to Scan" begin
    for (root, dir, files) in walkdir(resources)
        for file in files
            input = joinpath(root, file)
            scan  = NMAP.Scan(read(input))
            @test scan isa Scan
            xml   = XMLDict.parse_xml(String(read(joinpath(root, file))))
            @test xml[:args] == scan.args
        end
    end
end

include("utils.jl")
include("test_unmarshalling.jl")


text:
@testset "UnitCell" begin
    cell = UnitCell([5.0, 3.0, 1.0], [110.0, 120.0, 80.0])

    @test lengths(cell) ≈ [5.0, 3.0, 1.0]
    @test angles(cell) ≈ [110.0, 120.0, 80.0]

    cell = UnitCell([2.0, 3.0, 4.0])

    @test lengths(cell) == [2.0, 3.0, 4.0]
    @test angles(cell) == [90.0, 90.0, 90.0]

    @test volume(cell) == 2.0 * 3.0 * 4.0

    set_lengths!(cell, [10.0, 20.0, 30.0])
    @test lengths(cell) == [10.0, 20.0, 30.0]

    remove_chemfiles_warning() do
        @test_throws ChemfilesError set_angles!(cell, [80.0, 89.0, 100.0])
    end

    cell_matrix = [
        10.0 0.0  0.0
        0.0  20.0 0.0
        0.0  0.0  30.0
    ]
    @test matrix(cell) ≈ cell_matrix atol=1e-10

    @test shape(cell) == Chemfiles.Orthorhombic

    set_shape!(cell, Chemfiles.Triclinic)
    @test shape(cell) == Chemfiles.Triclinic

    set_angles!(cell, [80.0, 89.0, 100.0])
    @test angles(cell) == [80.0, 89.0, 100.0]

    copy = deepcopy(cell)
    @test lengths(copy) ≈ [10.0, 20.0, 30.0]

    set_lengths!(copy, [10.0, 10.0, 10.0])
    @test lengths(copy) ≈ [10.0, 10.0, 10.0]
    @test lengths(cell) ≈ [10.0, 20.0, 30.0]

    cell = UnitCell([10.0, 20.0, 30.0])
    @test wrap!(cell, [12.0, 22.0, -5.0]) == [2.0, 2.0, -5.0]

    cell = UnitCell(cell_matrix)
    @test lengths(cell) ≈ [10.0, 20.0, 30.0]
    @test angles(cell) ≈ [90.0, 90.0, 90.0]
end


text:
using ParallelAnalysis
using Random, Distributions, LinearAlgebra

Random.seed!(1234)
a = rand(Uniform(0.5, 2.0), 30)
b = [sort(rand(Uniform(-3, 3), 4); rev = false) for i in 1:30]
θ = rand(Normal(0, 1), 3000)

resp = generate_response(θ, a, b)
heuIRT = @time heuristicIRT(resp; method = :em)
@code_warntype heuristicIRT(resp)
heuIRT.a
heuIRT.d
heuIRT.b

test = polycor(resp)
ParallelAnalysis.replace_diagonal!(test)
eigvals(test)


text:
# Matplotlib utility functions.

export add_arrow

"""
    add_arrow(line, size=12)

Add an arrow of `size` to the object `line` from Matplotlib. This requires importing PyPlot,
and only works for Line2D.
"""
function add_arrow(line, size=12)

   color = line.get_color()
   xdata, ydata = line.get_data()

   for i = 1:2
      start_ind = length(xdata) ÷ 3 * i
      end_ind = start_ind + 1
      line.axes.annotate("",
         xytext=(xdata[start_ind], ydata[start_ind]),
         xy=(xdata[end_ind], ydata[end_ind]),
         arrowprops=Dict("arrowstyle"=>"-|>", "color"=>color),
         size=size
      )
   end

end


text:
using ComputerVisionMetrics
using ImageEdgeDetection
using Distances
using Test
using TestSetExtensions


text:
include("../imports.jl")

@testset ExtendedTestSet "fitonecycle!" begin
    learner = testlearner(Recorder())
    @test_nowarn fitonecycle!(learner, 5)
end

@testset ExtendedTestSet "decay_optim" begin
    optim = ADAM()
    @test decay_optim(optim, 0.1) isa Optimiser
    @test decay_optim(Optimiser(ADAM(), ADAM()), 0.1) isa Optimiser
    @test decay_optim(optim, 0.1).os[1] isa WeightDecay
    o = decay_optim(Optimiser(ADAM(), WeightDecay(.5)), .1)
    @test o.os[1] isa WeightDecay
    @test o.os[2] isa ADAM
end


text:
using Distributions
using Test

# Test for gradlogpdf on univariate distributions

@test isapprox(gradlogpdf(Beta(1.5, 3.0), 0.7)    , -5.9523809523809526 , atol=1.0e-8)
@test isapprox(gradlogpdf(Chi(5.0), 5.5)          , -4.7727272727272725 , atol=1.0e-8)
@test isapprox(gradlogpdf(Chisq(7.0), 12.0)       , -0.29166666666666663, atol=1.0e-8)
@test isapprox(gradlogpdf(Exponential(2.0), 7.0)  , -0.5                , atol=1.0e-8)
@test isapprox(gradlogpdf(Gamma(9.0, 0.5), 11.0)  , -1.2727272727272727 , atol=1.0e-8)
@test isapprox(gradlogpdf(Gumbel(3.5, 1.0), 4.0)  , -1.6065306597126334 , atol=1.0e-8)
@test isapprox(gradlogpdf(Laplace(7.0), 34.0)     , -1.0                , atol=1.0e-8)
@test isapprox(gradlogpdf(Logistic(-6.0), 1.0)    , -0.9981778976111987 , atol=1.0e-8)
@test isapprox(gradlogpdf(LogNormal(5.5), 2.0)    ,  1.9034264097200273 , atol=1.0e-8)
@test isapprox(gradlogpdf(Normal(-4.5, 2.0), 1.6) , -1.525              , atol=1.0e-8)
@test isapprox(gradlogpdf(TDist(8.0), 9.1)        , -0.9018830525272548 , atol=1.0e-8)
@test isapprox(gradlogpdf(Weibull(2.0), 3.5)      , -6.714285714285714  , atol=1.0e-8)
@test isapprox(gradlogpdf(Uniform(-1.0, 1.0), 0.3),  0.0                , atol=1.0e-8)

# Test for gradlogpdf on multivariate distributions

@test isapprox(gradlogpdf(MvNormal([1., 2.], [1. 0.1; 0.1 1.]), [0.7, 0.9])   ,
    [0.191919191919192, 1.080808080808081]   ,atol=1.0e-8)
@test isapprox(gradlogpdf(MvTDist(5., [1., 2.], [1. 0.1; 0.1 1.]), [0.7, 0.9]),
    [0.2150711513583442, 1.2111901681759383] ,atol=1.0e-8)


text:
export IntervalCurve

struct IntervalCurve{S<:Space,T,VT} <: SegmentDomain{T}
    curve::Fun{S,T,VT}
end

==(a::IntervalCurve, b::IntervalCurve) = a.curve == b.curve

isempty(::IntervalCurve) = false

points(c::IntervalCurve, n::Integer) = c.curve.(points(domain(c.curve),n))

checkpoints(d::IntervalCurve) = fromcanonical.(Ref(d),checkpoints(domain(d.curve)))

for op in (:(leftendpoint),:(rightendpoint),:(rand))
    @eval $op(c::IntervalCurve) = c.curve($op(domain(c.curve)))
end

canonicaldomain(c::IntervalCurve) = domain(c.curve)

fromcanonical(c::IntervalCurve{S,T},x) where {S<:Space,T<:Number} = c.curve(x)
function tocanonical(c::IntervalCurve,x)
    rts=roots(c.curve-x)
    @assert length(rts)==1
    first(rts)
end

fromcanonicalD(c::IntervalCurve,x)=differentiate(c.curve)(x)

function indomain(x,c::IntervalCurve)
    rts=roots(c.curve-x)
    if length(rts) ≠ 1
        false
    else
        in(first(rts),canonicaldomain(c))
    end
end

isambiguous(d::IntervalCurve) = ncoefficients(d.curve)==0 && isambiguous(domain(d.curve))

reverseorientation(d::IntervalCurve) = IntervalCurve(reverseorientation(d.curve))
convert(::Type{IntervalCurve{S,T}},::AnyDomain) where {S,T}=Fun(S(AnyDomain()),[NaN])

arclength(d::IntervalCurve) = linesum(ones(d))


text:
module ConstantArrays

export ConstantArray, getindex!

"""
    ConstantArray{T,D} <: AbstractArray{T,D}
- `value::T`
- `dim::Dims{D}`
"""
struct ConstantArray{T,D} <: AbstractArray{T,D}
    value::T
    dim::Dims{D}
    ConstantArray(value::T, dim::Dims{D}) where {T,D} =
        new{T, D}(value, dim)
end

"""
    ConstantArray(value, dim::Dims)
    ConstantArray(dim::Dims) (default value=1)

Return a read-only `AbstractArray` of dimension `dim`
as if it is filled with the constant `value`,
i.e., equivalent to `fill(value, dim)`
except it uses much less memory.

The `setindex!` operation is not allowed for `ConstantArray`

# Examples
```jldoctest
julia> x = ConstantArray(7, (2,3))
2×3 ConstantArray{Int64,2}:
 7  7  7
 7  7  7

julia> x[1]
7

julia> x[1] = 10
ERROR: setindex! not defined for ConstantArray{Int64,2}
[...]
```
"""
ConstantArray(dim::Dims{D}) where {D} = ConstantArray(one(Bool), dim)

Base.IteratorSize(::Type{<:ConstantArray{T,N}}) where {T,N} = Base.HasLength()
Base.IteratorEltype(::Type{<:ConstantArray{T,N}}) where {T,N} = T
Base.eltype(::Type{<:ConstantArray{T,N}}) where {T,N} = T
Base.size(x::ConstantArray) = x.dim
Base.size(x::ConstantArray, args...) = x.dim[args...]

# probably pointless because isbitstype(ConstantArray{T,D}) is true:
Base.copy(x::ConstantArray) = ConstantArray(x.value, x.dim)

Base.@propagate_inbounds function Base.getindex(x::ConstantArray, i::Int)
    @boundscheck 1 <= i <= prod(x.dim) || throw(BoundsError("$i"))
    return x.value
end

# harder to handle ':' etc. so let Base do it
#Base.@propagate_inbounds Base.getindex(x::ConstantArray, I...) = fill(x.value, ?)

Base.firstindex(::ConstantArray) = 1
Base.lastindex(x::ConstantArray) = prod(x.dim)
Base.IndexStyle(::Type{<:ConstantArray{T,N}}) where {T,N} = IndexLinear()
#Base.iterate(x::ConstantArray, args...) = iterate(x.value, args...)
Base.length(x::ConstantArray) = prod(x.dim)

Base.axes(x::ConstantArray) = ntuple(i -> Base.OneTo(x.dim[i]), length(x.dim))

"""
    getindex(x::AbstractArray{T,D}, mask::ConstantArray{Bool,D})

Ensure that `x[mask]` is efficient.
"""
Base.getindex(x::AbstractArray{T,D}, mask::ConstantArray{Bool,D}) where {T,D} =
mask.value ? vec(x) : ones(T, 0)

"""
    getindex!(y::AbstractVector, x::AbstractArray{T,D}, mask::ConstantArray{Bool,D})
Equivalent to the in-place `y .= x[mask]` but is non-allocating.

For non-Boolean indexing, just use `@views y .= A[index]`, per
https://discourse.julialang.org/t/efficient-non-allocating-in-place-getindex-for-bitarray/42268
"""
@inline function getindex!(
    y::AbstractVector,
    x::AbstractArray{T,D},
    mask::ConstantArray{Bool,D},
) where {T,D}
    axes(y) == (Base.OneTo(length(mask)),) || throw("y axes $(axes(y))")
    axes(mask) == axes(x) || throw(DimensionMismatch("x vs mask"))
    mask[1] == true || fail("only valid for mask with trues")
    y .= vec(x)
    return y
end

"""
    show(io::IO, x::ConstantArray)
    show(io::IO, ::MIME"text/plain", x::ConstantArray)
Concise description.
"""
Base.show(io::IO, x::ConstantArray) =
    print(io, "$(summary(x)) = $(x.value)")

function Base.show(io::IO, ::MIME"text/plain", x::ConstantArray)
    show(io, x)
end

end # module


text:
# ---
# title: 318. Maximum Product of Word Lengths
# id: problem318
# author: Indigo
# date: 2021-06-10
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/maximum-product-of-word-lengths/description/>
# hidden: true
# ---
#
# Given a string array `words`, find the maximum value of `length(word[i]) *
# length(word[j])` where the two words do not share common letters. You may
# assume that each word will contain only lower case letters. If no such two
# words exist, return 0.
#
# **Example 1:**
#
#
#
#     Input: ["abcw","baz","foo","bar","xtfn","abcdef"]
#     Output:16
#     Explanation:The two words can be "abcw", "xtfn".
#
# **Example 2:**
#
#
#
#     Input: ["a","ab","abc","d","cd","bcd","abcd"]
#     Output:4
#     Explanation:The two words can be "ab", "cd".
#
# **Example 3:**
#
#
#
#     Input: ["a","aa","aaa","aaaa"]
#     Output:0
#     Explanation:No such pair of words.
#
#
#
#
# **Constraints:**
#
#   * `0 <= words.length <= 10^3`
#   * `0 <= words[i].length <= 10^3`
#   * `words[i]` consists only of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function max_product(words::Vector{String})
    mp = Dict{UInt32, Int}()
    for word in words
        init = zero(UInt32)
        for c in word
            init |= UInt32(1) << (c - 'a')
        end
        mp[init] = max(length(word), get(mp, init, 0))
    end
    res = 0
    for x in keys(mp), y in keys(mp)
        (x & y == 0) && (res = max(res, mp[x] * mp[y]))
    end
    return res
end
## @lc code=end


text:
export FunctionalCoder, InputsCoder, MLPCoder

"""
    struct FunctionalCoder

A coder that codes to a discretisation for a functional representation.

# Fields
- `disc::Discretisation`: Discretisation for the functional representation.
- `coder`: Coder.
"""
struct FunctionalCoder
    disc::Discretisation
    coder
end

@Flux.functor FunctionalCoder

code(c::FunctionalCoder, xz, z, x; kws...) =
    code(c.coder, xz, z, c.disc(xz, x; kws...); kws...)

function code_track(c::FunctionalCoder, xz, z, x, h; kws...)
    x_disc = c.disc(xz, x; kws...)
    return code_track(c.coder, xz, z, x_disc, vcat(h, [x_disc]); kws...)
end

recode(c::FunctionalCoder, xz, z, h; kws...) =
    recode(c.coder, xz, z, h[2:end]; kws...)

"""
    struct InputsCoder

Code with the target inputs.
"""
struct InputsCoder end

@Flux.functor InputsCoder

code(encoder::InputsCoder, xz, z, x::AA; kws...) = x, x

"""
    struct MLPCoder

Code with an MLP.

# Fields
- `mlp₁`: Pre-pooling MLP.
- `mlp₂`: Post-pooling MLP.
"""
struct MLPCoder
    mlp₁
    mlp₂
end

@Flux.functor MLPCoder

code(encoder::MLPCoder, xz::AA, z::AA, x::AA; kws...) =
    x, encoder.mlp₂(mean(encoder.mlp₁(cat(xz, z, dims=2)), dims=1))

function code(encoder::MLPCoder, xz::Nothing, z::Nothing, x::AA; kws...)
    batch_size = size(x, 3)
    r = zeros_gpu(Float32, 1, encoder.mlp₁.dim_out, batch_size)
    return x, encoder.mlp₂(r)
end


text:
# Contains evaluation functions
# ----------------------------------
"""
	 Evaluates the FIS

	eval_fis(fis, input_values,defuzz_method = "WTAV")

	 Parameters
	 ----------
	 `fis` is the inference system to evaluate
	 `input_values` is a Vector of inputs
	 `defuzz_method` is the method for defuzzification, see defuzz function definition

"""
function eval_fis(fis::FISMamdani, input_values::Vector{<:AbstractFloat},defuzz_method = "WTAV")

	firing_strengths = AbstractFloat[]
	for rule in fis.rules
		tmp_strengths = AbstractFloat[]
		for i in 1:length(rule.input_mf_names)
			push!(tmp_strengths, fis.input_mfs_dicts[i][rule.input_mf_names[i]].eval(input_values[i]))
		end
			push!(firing_strengths, firing(tmp_strengths,rule.firing_method ))
	end
	defuzz(firing_strengths, fis.rules, fis.output_mfs_dict, defuzz_method)

end

"""
	 Evaluates the FIS

	eval_fis(fis, input_values)

	 Parameters
	 ----------
	 `fis` is the inference system to evaluate
	 `input_values` is a Vector of inputs

"""
function eval_fis(fis::FISSugeno, input_values::Vector{<:AbstractFloat})

	firing_strengths = AbstractFloat[]
	for rule in fis.rules
		tmp_strengths = AbstractFloat[]
		for i in 1:length(rule.input_mf_names)
			push!(tmp_strengths, fis.input_mfs_dicts[i][rule.input_mf_names[i]].eval(input_values[i]))
		end
			push!(firing_strengths, firing(tmp_strengths,rule.firing_method ))
	end

	push!(input_values, 1)
	n_firing_strengths = firing_strengths / sum(firing_strengths)
	out = 0.0
	for i = 1:length(fis.rules)
		out += n_firing_strengths[i] * (fis.rules[i].output_mf' * input_values)[1]
	end
	pop!(input_values)
	return out

end

"""
	 Defuzzifies the output using the given firing strengths

	defuzz(firing_strengths, rules, output_mfs_dict, defuzz_method)

	 Parameters
	 ----------
	 `firing_strengths` is a Vector of firing strengths
	 		one for each output membership function
	 `rules` is a Vector of Rule
	 `output_mfs_dict` is a Dict of output membership functions
	 `defuzz_method` is the method for defuzzification
	 		"MOM" - Mean of Maximum
	 		"WTAV" - Weighted Average
"""
function defuzz(firing_strengths::Vector{AbstractFloat}, rules::Vector{Rule},	output_mfs_dict::Dict{AbstractString, MF}, defuzz_method::AbstractString)

	if defuzz_method == "MOM"
		max_firing_index = indmax(firing_strengths)
		max_fired_mf_name = rules[max_firing_index].output_mf
		output_mfs_dict[max_fired_mf_name].mean_at(maximum(firing_strengths))
	elseif defuzz_method == "WTAV"
		mean_vec = AbstractFloat[]
		for i in 1:length(rules)
			push!(mean_vec, output_mfs_dict[rules[i].output_mf].mean_at(firing_strengths[i]))
		end
		(mean_vec' * firing_strengths)[1] / sum(firing_strengths)
	end

end

"""
	 Docs goes here
"""
function firing(tmp_strengths::Vector{<:AbstractFloat},firing_method::AbstractString)
	if firing_method == "MIN"
		return	minimum_value(tmp_strengths)
	elseif firing_method == "A-PROD"
		return algebraic_product(tmp_strengths)
	elseif firing_method == "B-DIF"
		return bounded_difference(tmp_strengths)
	elseif firing_method == "D-PROD"
		return drastic_product(tmp_strengths)
	elseif firing_method == "E-PROD"
		return einstein_product(tmp_strengths)
	elseif firing_method == "H-PROD"
		return hamacher_product(tmp_strengths)
	elseif firing_method == "MAX"
		return	maximum_value(tmp_strengths)
	elseif firing_method == "A-SUM"
		return algebraic_sum(tmp_strengths)
	elseif firing_method == "B-SUM"
		return bounded_sum(tmp_strengths)
	elseif firing_method == "D-SUM"
		return drastic_sum(tmp_strengths)
	elseif firing_method == "E-SUM"
		return einstein_sum(tmp_strengths)
	elseif firing_method == "H-SUM"
		return hamacher_sum(tmp_strengths)
	end
end


text:
"""
    test_linalg(A, b, c, uind, uval, h)

Verifies that the linear algebra operations are properly defined for the matrix
    structure given in input.
"""
function test_linalg(
    A::AbstractMatrix{Tv},
    b::AbstractVector{Tv},
    c::AbstractVector{Tv},
    uind::AbstractVector{Ti},
    uval::AbstractVector{Tv},
    h::Tulip.PrimalDualPoint{Tv}
) where{Tv<:Real, Ti<:Integer}

    # Dimension check
    m = size(A, 1)
    n = size(A, 2)
    n == size(c, 1) || throw(DimensionMismatch(""))
    m == size(b, 1) || throw(DimensionMismatch(""))
    n == size(h.x, 1) || throw(DimensionMismatch(""))
    m == size(h.y, 1) || throw(DimensionMismatch(""))
    n == size(h.s, 1) || throw(DimensionMismatch(""))
    size(uind, 1) == size(uval, 1) || throw(DimensionMismatch(""))

    # matrix-vector multiplication
    A * h.x;
    Base.LinAlg.At_mul_B(A, b)
    Base.LinAlg.A_mul_B!(h.y, A, h.x)

    # Cholesky factorization
    F = Tulip.symbolic_cholesky(A)
    d = 1.0 + rand(n)
    Tulip.Cholesky.cholesky!(A, d, F)  # in-place update

    # solve linear system
    y = F \ b

    return nothing
end

# SparseCSC matrix
srand(0)
m = 4
n = 16
c = rand(n)
b = rand(m)
u = sprand(n, 0.75)
uind = u.nzind
uval = u.nzval
p = nnz(u)
A = sprand(m, n, 1.0)
h = Tulip.PrimalDualPoint(zeros(n), zeros(p), zeros(m), zeros(n), zeros(p))
test_linalg(A, b, c, uind, uval, h)

# BlockAngular matrix
srand(0)
nblocks = 4
cols = [rand(m, 4) for _ in 1:nblocks]
A = Tulip.Cholesky.DenseBlockAngular(cols)
(m, n) = size(A)
c = rand(n)
b = rand(m)
u = sprand(n, 0.75)
uind = u.nzind
uval = u.nzval
p = nnz(u)

h = Tulip.PrimalDualPoint(zeros(n), zeros(p), zeros(m), zeros(n), zeros(p))

test_linalg(A, b, c, uind, uval, h)


text:
using sem, Arrow, ModelingToolkit, LinearAlgebra,
    SparseArrays, DataFrames, Optim, LineSearches,
    Statistics

cd("test")

## Observed Data
dat = DataFrame(Arrow.Table("comparisons/data_dem.arrow"))
par_ml = DataFrame(Arrow.Table("comparisons/par_dem_ml.arrow"))
par_ls = DataFrame(Arrow.Table("comparisons/par_dem_ls.arrow"))

dat =
    select(
        dat,
        [:x1, :x2, :x3, :y1, :y2, :y3, :y4, :y5, :y6, :y7, :y8])
# observed
semobserved = SemObsCommon(data = Matrix{Float64}(dat))

#diff
diff_fin = SemFiniteDiff(
    BFGS(),
    Optim.Options())

diff_fin_new = SemFiniteDiff(
    BFGS(),
    Optim.Options())

## Model definition
@ModelingToolkit.variables x[1:31]

#x = rand(31)

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14)]

start_val_ml = Vector{Float64}(par_ml.start[par_order])
start_val_ls = Vector{Float64}(par_ls.start[par_order])
start_val_snlls = Vector{Float64}(par_ls.start[par_order][21:31])

# loss
loss_ml = Loss([SemML(semobserved, [0.0], similar(start_val_ml))])
loss_ls = Loss([sem.SemWLS(semobserved, [0.0], similar(start_val_ml))])
loss_snlls = Loss([sem.SemSWLS(semobserved, [0.0], similar(start_val_ml))])

# start_val_ml = [fill(1.0, 14); fill(0, 6); fill(1, 8); fill(0, 3)]
# start_val_ls = [fill(1.0, 14); fill(0, 6); fill(1, 8); fill(0, 3)]
# start_val_snlls = [fill(1.0, 8); fill(0, 3)]

# imply
imply_ml = ImplySymbolic(A, S, F, x, start_val_ml)
imply_ls = sem.ImplySymbolicWLS(A, S, F, x, start_val_ls)
imply_snlls = sem.ImplySymbolicSWLS(A, S, F, x[21:31], start_val_snlls)

# model
model_ml = Sem(semobserved, imply_ml, loss_ml, diff_fin)
model_ls = Sem(semobserved, imply_ls, loss_ls, diff_fin)
model_snlls = Sem(semobserved, imply_snlls, loss_snlls, diff_fin_new)

# fit
solution_ml = sem_fit(model_ml)
solution_ls = sem_fit(model_ls)
solution_snlls = sem_fit(model_snlls)

all(#
        abs.(solution_ml.minimizer .- par_ml.est[par_order]
            ) .< 0.05*abs.(par_ml.est[par_order]))

all(#
            abs.(solution_ls.minimizer .- par_ls.est[par_order]
                ) .< 0.05*abs.(par_ls.est[par_order]))

all(#
            abs.(solution_snlls.minimizer .- par_ls.est[par_order][21:31]
                ) .< 0.05*abs.(par_ls.est[par_order][21:31]))

model_snlls(solution_ls.minimizer[21:31])

t = solution_ls.minimizer[21] .+ collect(-0.5:0.05:0.5)

t_vec = [[t_i; solution_ls.minimizer[22:31]...] for t_i in t]

y = model_snlls.(t_vec)

using Plots

plot(t, y, xticks = t[1:5:21])

function get_fvec(model, pos, solution, iterator)
    t = solution[pos] .+ collect(iterator)
    t_vec = [copy(solution) for i in t]
    for i = 1:size(t, 1)
        t_vec[i][pos] = t[i]
    end
    y = [model(t_vec_i) for t_vec_i in t_vec]
    return plot(t, y, xticks = t[1:5:21])
end

get_fvec(model_snlls, 10, solution_snlls.minimizer, -0.5:0.05:0.5)

plots = [get_fvec(model_snlls, i, solution_snlls.minimizer, -0.5:0.05:0.5) for i in 1:11]

plots[11]

plot(plots...)

solution_snlls.minimizer - par_ls.est[par_order][21:31]

using BenchmarkTools

@benchmark solution_ml = sem_fit(model_ml)
@benchmark solution_ls = sem_fit(model_ls)
@benchmark solution_snlls = sem_fit(model_snlls)

@benchmark model_ml(start_val_ml)
@benchmark model_ls(start_val_ls)
@benchmark model_snlls(start_val_snlls)

############### NLopt

diff_nlopt = SemFiniteDiff(
    :LD_LBFGS,
    nothing)

start_val_snlls = [fill(1.0, 8); fill(0, 3)]
imply_snlls = sem.ImplySymbolicSWLS(A, S, F, x[21:31], start_val_snlls)

model_ml = Sem(semobserved, imply_ml, loss_ml, diff_nlopt)
model_ls = Sem(semobserved, imply_ls, loss_ls, diff_nlopt)
model_snlls = Sem(semobserved, imply_snlls, loss_snlls, diff_nlopt)

solution_ml =
    sem.sem_fit_nlopt(
        model_ml;
        lower = fill(-11.0, 31),
        upper = fill(10.0, 31),
        local_algo = :LD_LBFGS)[2]
solution_ls =
    sem.sem_fit_nlopt(
        model_ls;
        lower = fill(-11.0, 31),
        upper = fill(10.0, 31))[2]
solution_snlls =
    sem.sem_fit_nlopt(
        model_snlls;
        lower = fill(-11.0, 11),
        upper = fill(10.0, 11))[2]

all(#
    abs.(solution_ml .- par_ml.est[par_order]
        ) .< 0.05*abs.(par_ml.est[par_order]))

all(#
    abs.(solution_ls .- par_ls.est[par_order]
        ) .< 0.05*abs.(par_ls.est[par_order]))

all(#
    abs.(solution_snlls .- par_ls.est[par_order][21:31]
        ) .< 0.05*abs.(par_ls.est[par_order][21:31]))

@benchmark solution_ls =
    sem.sem_fit_nlopt(
        model_ls;
        lower = fill(-11.0, 31),
        upper = fill(10.0, 31))[2]


text:
"""
$(FUNCTIONNAME)(grid_x::AbstractRange, ils_in::Array{FT}, ils_Δ::Array{FT}, extended_dims::Array{Int}=[]) where {FT <: AbstractFloat}

Pre-compute the ILS table input as function of spectral distance from center converted to the modeling grid
Input: grid_x, ils_in, ils_Δ, extended_dims
Output: Offset-Array with tabulated responses interpolated to grid_x
"""
function prepare_ils_table(grid_x::AbstractRange, ils_in::Array{FT}, ils_Δ::Array{FT}, extended_dims::Array{Int}=[]) where {FT <: AbstractFloat}
    @assert minimum(abs.(grid_x)) < eps(FT) "grid_x must include 0 (center pixel)"
    ind_0 = argmin(abs.(grid_x))
    axis_pixel = (-ind_0 + 1):(grid_x.len - ind_0)
    # Dimension of ILS table (at least 1D, first dimension needs to be across wavenumber/wavelength)
    dims  = size(ils_in);
    # number of spectral positions of ILS table
    n_x = dims[1]
    # Number of ILS per detector position (can be 1 if ILS is constant across detector grid)
    n_pos = dims[2]
    ils    = view(ils_in, :, :,  extended_dims...);
    ils_Δ_ = view(ils_Δ,  :, :,  extended_dims...);
    ils_pixel = zeros(FT, grid_x.len, n_pos);

    for i = 1:n_pos
        ind = findall(minimum(ils_Δ_[:,i]) .< grid_x .< maximum(ils_Δ_[:,i]));
        interp = Interpolations.LinearInterpolation(ils_Δ_[:,i], ils[:,i])
        ils_pixel[ind,i] = interp.(grid_x[ind]);
    end
    # normalize here
    #return OffsetArray(ils_pixel ./ sum(ils_pixel, dims=1), axis_pixel, 1:n_pos)
    return OffsetArray(ils_pixel * FT(grid_x.step) , axis_pixel, 1:n_pos)
end
# This can still derive min/max range automatically, need to double check

function create_instrument_kernel(FTS::FTSInstrument, grid_x::AbstractRange, ν̅)
    @unpack MOPD, β, FOV = FTS
    # Need to make sure this is 0!
    i₀ = argmin(abs.(grid_x))
    axis_pixel = (-i₀ + 1):(grid_x.len - i₀)
    @info "Center of OffsetArrays is " minimum(abs.(grid_x))
    #Δσ = 1/(2FTS.OPD)

    # Create sinc FTS kernel
    XX = 2MOPD * grid_x
    sinc_kernel = OffsetArray(2MOPD * sinc.(XX),axis_pixel)

    # Now to FOV box:
    Θ = ν̅/2 * FOV^2
    β = β * Θ
    @show Θ, β
    # Create empty box_kernel:
    box = collect(0*grid_x);
    for i in eachindex(grid_x)
        if grid_x[i] ≤ 0
            if Θ > β && grid_x[i] ≥ -(Θ - β)
                box[i] = 1
            elseif grid_x[i] ≥ -(Θ + β)
                box[i] = acos((grid_x[i]^2 + β^2 - Θ^2)/(2grid_x[i] * β)) / π
            end
        end
    end
    box_kernel = OffsetArray(box,axis_pixel)
    box_kernel /= sum(box_kernel)
    sinc_kernel /= sum(sinc_kernel)
    return imfilter(box_kernel,sinc_kernel)
end

"Create kernel from a single continous distribution function"
function create_instrument_kernel(di::ContinuousUnivariateDistribution, grid_x::AbstractRange)
    ils = pdf.(di,grid_x)
    i₀ = argmin(abs.(grid_x))
    axis_pixel = (-i₀ + 1):(grid_x.len - i₀)
    ils /= sum(ils)
    return OffsetArray(ils,axis_pixel)
end

"Create kernel from multiple single continous distribution function (convolves those)"
function create_instrument_kernel(di::Array{T}, grid_x::AbstractRange) where T <: ContinuousUnivariateDistribution
    ils = pdf.(di[1],grid_x)
    for i=2:length(di)
        # might change this in future to just do multiplications in fft space and then ifft back
        ils = imfilter(ils, pdf.(di[i],grid_x))
    end
    i₀ = argmin(abs.(grid_x))
    axis_pixel = (-i₀ + 1):(grid_x.len - i₀)
    ils /= sum(ils)
    return OffsetArray(ils,axis_pixel)
end


text:
#=
Created on 20/01/2021 16:45:00
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Illustration of the composite and parametric types:
- composite types
    - mutable/immuatble
- unions
- Parametric types
=#

#=
Composite types, sometimes call records, structs or object, can store several values in its *fields*.

When defining a new composite type, we can choose them to be mutable or immuatble:
- mutable types are defined using `mutable struct ... end`, they allow the fields to be changed after the object is created;
- immuatble types are defined similarly using `mutable struct ... end`, after creating the object its fields cannot be changed.
Mutable types are a bit more flexible, though might be a somewhat less safe and are more difficult to work with. As the compiler
knows everything in advance, it might better optimize for immutable types. Which one you choose depends on your application, though
generally immuatable types are the better choice!

As an example, let us define an agent type for an ecological individual-based model (IBM). We create the abstract type `Agent` for which we
can then specify several children types.
=#

abstract type Agent end

#=
The concrete types in such an IBM might represent an animal type you want to model, for example preys and predators. Making a concrete type of
a prey animal, we want each to have an unique identifier (represented by an integer) and a position. As we expect the agent to move, hence changing
its position when our simulation runs, we choose a mutable type.
=#

mutable struct Prey <: Agent
    id::Int
    pos
end
# notice the type annotation for `id`, which we choose to always reprsenent by an integer.
# the position might be represented by (x,y) coordinates, or as a location ID, or a position on a grid. We don't know, so we leave it untyped.

# defining a composite type immediately a constructor available.
deer = Prey(1, (0.5, 1.9))

# you can always check which field names are available
fieldnames(Prey)

# the fields can be accessed simply
deer.id, deer.pos

# FYI: This is just syntactic sugar for `getfield`, e.g. `getfield(deer, :id)`

# Similarly, a predator type can be defined. In addition to an id and position, which each agent has, they also have a size, determining its mobility.

mutable struct Predator <: Agent
    id::Int
    pos
    size::Float64
end

wolf = Predator(2, (0.0, 0.0), 40.0)  # 40 kg wolf

#=
Using the `.` syntax for accessing the fields is not very tidy! We should define custom getter functions
for the user to access the relevant fields. We could define `id` and `position` methods to get the respective
fields for the two agents. However, since these fields should be defined for every `Agent` type, we can just create
these for the Agent type!
=#

id(agent::Agent) = agent.id
position(agent::Agent) = agent.pos

# Here, we could theoretically have ommited the type annotation in the function. Then the function would accept
# objects of the non-agent type and likely yield an error because they don't have the `id` or `pos` field. Now,
# these functions will return a `MethodError` when given a non-`Agent` input.

# A slightly more interesting example is by extending `size`.

Base.size(agent::Predator) = agent.size

# Here, we had to import `size` because we are extending a function from the `Base` library to work with a new type (doing something vastly different
# than its original function).

# similarly, we can program behaviour between the agents

interact(agent1::Agent, agent2::Agent) = nothing
interact(agent1::Predator, agent2::Prey) = "eat"
interact(agent1::Prey, agent2::Predator) = "run"

# We have chosen the default behaviour that two Agents of unspecified types do not interact at al,
# this will now be the case when a prey meets other prey, a predator an other predator or a new third type comes into the equation.

# FYI: Since in these simple examples, the `interact` methods do not use their arguments, merely perform type checking, we could have written this as `interact(::Agent,::Agent) = ...` etc.

# PARAMETRIC TYPES

#=
Sometimes we want more flexiblility in defining types. Think of designing a new type of matrix. Here you would like to work them for all
numeric datatypes, Int, Int8, Float6, Rational, in addition to new datatypes that might not even be defined yet! To this end, we use
*parametric types*, types that depend on another type.

For example, consider a 2-dimensional coordinate:
=#

struct Point{T}
    x::T
    y::T
end

# Here, each coordinate of the type `Point` has two attributes, `x` and `y`, of the same type. The specific type of Point can vary.

p = Point(1.0, 2.0)

Point(1, 2)

# note that

p isa Point

# But what will happen if you evaluate `Point(1, 2.0)`?

Point(1, 2.0)

# Parametric types can be used in dispatch. For example, if we want to compute the norm of a Point, this would only make sense if Point is a number.

norm(p::Point{T} where {T<:Number}) = sqrt(p.x^2 + p.y^2)

norm(p)

# Constructors

## Outer constructors

# Constructors are functions that create new objects. We have already seen that when creating a new `struc`, this immediately initiates the constructor (e.g., `Point(1.0, 2.0)`). These can also be made explicitly:

Point(x::T, y::T) where {T<:Real} = Point{T}(x,y)

# Constructors, however, allow us to have custom behavior when initializing types. For example, we have seen that `Point(1, 2.0)` won't work, because the two inputs are of the same type.
# In this case, we can make the rule that one of the inputs has to be promoted to a more general type.

Point(x::Real, y::Real) = Point(promote(x, y)...)

Point(1, 2.0)

# We can write other constructors just like functions. For example, support that when we provide a single `x`, we want to create a point (x, y):

Point(x) = Point(x, x);

Point(1)

## Inner constructors

#The above examples show *outer constructors*. These are defined outside the structure. We can also use *inner constructors*, which are declared within the definition of the type. These make use of the keyword `new`. For example, let us define an ordered pair.

struct OrderedPair
  x
  y
  function OrderedPair(x, y)
    if x < y
      new(x, y)
    else
      new(y, x)
    end
  end
end

OrderedPair(18, 23)

OrderedPair(8, 2)

# FYI: for parametric types, the `new` keyword should be type annotated. So, for in the `Point` example one would use `new{T}(x,y)`.


text:
using Neighborhood
using Distances

export uzal_cost
export uzal_cost_local

"""
    uzal_cost(Y::Dataset; kwargs...) → L
Compute the L-statistic `L` for input dataset `Y` according to Uzal et al.[^Uzal2011], based on
theoretical arguments on noise amplification, the complexity of the
reconstructed attractor and a direct measure of local stretch which constitutes
an irrelevance measure. It serves as a cost function of a state space
trajectory/embedding and therefore allows to estimate a "goodness of a
embedding" and also to choose proper embedding parameters, while minimizing
`L` over the parameter space. For receiving the local cost function `L_local`
(for each point in state space - not averaged), use `uzal_cost_local(...)`.

## Keyword arguments

* `samplesize = 0.5`: Number of considered fiducial points v as a fraction of
  input state space trajectory `Y`'s length, in order to average the conditional
  variances and neighborhood sizes (read algorithm description) to produce `L`.
* `K = 3`: the amount of nearest neighbors considered, in order to compute σ_k^2
  (read algorithm description).
  If given a vector, minimum result over all `k ∈ K` is returned.
* `metric = Euclidean()`: metric used for finding nearest neigbhors in the input
  state space trajectory `Y.
* `w = 1`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `Tw = 40`: The time horizon (in sampling units) up to which E_k^2 gets computed
  and averaged over (read algorithm description).

## Description
The `L`-statistic is based on theoretical arguments on noise amplification, the
complexity of the reconstructed attractor and a direct measure of local stretch
which constitutes an irrelevance measure. Technically, it is the logarithm of
the product of `σ`-statistic and a normalization statistic `α`:

L = log10(σ*α)

The `σ`-statistic is computed as follows. `σ = √σ² = √(E²/ϵ²)`.
`E²` approximates the conditional variance at each point in state space and
for a time horizon `T ∈ Tw`, using `K` nearest neighbors. For each reference
point of the state space trajectory, the neighborhood consists of the reference
point itself and its `K+1` nearest neighbors. `E²` measures how strong
a neighborhood expands during `T` time steps. `E²` is averaged over many time
horizons `T = 1:Tw`. Consequently, `ϵ²` is the size of the neighborhood at the
reference point itself and is defined as the mean pairwise distance of the
neighborhood. Finally, `σ²` gets averaged over a range of reference points on
the attractor, which is controlled by `samplesize`. This is just for performance
reasons and the most accurate result will obviously be gained when setting
`samplesize=1.0`

The `α`-statistic is a normalization factor, such that `σ`'s from different
embeddings can be compared. `α²` is defined as the inverse of the sum of
the inverse of all `ϵ²`'s for all considered reference points.

[^Uzal2011]: Uzal, L. C., Grinblat, G. L., Verdes, P. F. (2011). [Optimal reconstruction of dynamical systems: A noise amplification approach. Physical Review E 84, 016223](https://doi.org/10.1103/PhysRevE.84.016223).
"""
function uzal_cost(Y::AbstractDataset{D, ET};
        Tw::Int = 40, K::Int = 3, w::Int = 1, samplesize::Real = 0.5,
        metric = Euclidean()
    ) where {D, ET}

    # select a random state space vector sample according to input samplesize
    NN = length(Y)-Tw;
    NNN = floor(Int, samplesize*NN)
    ns = sample(1:NN, NNN; replace=false) # the fiducial point indices

    vs = Y[ns] # the fiducial points in the data set

    vtree = KDTree(Y[1:end-Tw], metric)
    allNNidxs, allNNdist = all_neighbors(vtree, vs, ns, K, w)
    ϵ² = zeros(NNN)             # neighborhood size
    E²_avrg = zeros(NNN)        # averaged conditional variance
    E² = zeros(Tw)
    ϵ_ball = zeros(ET, K+1, D) # preallocation
    u_k = zeros(ET, D)

    # loop over each fiducial point
    for (i,v) in enumerate(vs)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v
        # pairwise distance of fiducial points and `v`
        pdsqrd = fiducial_pairwise_dist_sqrd(view(Y.data, NNidxs), v, metric)
        ϵ²[i] = (2/(K*(K+1))) * pdsqrd  # Eq. 16
        # loop over the different time horizons
        for T = 1:Tw
            E²[T] = comp_Ek2!(ϵ_ball, u_k, Y, ns[i], NNidxs, T, K, metric) # Eqs. 13 & 14
        end
        # Average E²[T] over all prediction horizons
        E²_avrg[i] = mean(E²)                   # Eq. 15
    end
    σ² = E²_avrg ./ ϵ² # noise amplification σ², Eq. 17
    σ²_avrg = mean(σ²) # averaged value of the noise amplification, Eq. 18
    α² = 1 / mean(ϵ².^(-1)) # for normalization, Eq. 21
    L = log10(sqrt(σ²_avrg)*sqrt(α²))
end

function fiducial_pairwise_dist_sqrd(fiducials, v, metric)
    pd = zero(eltype(fiducials[1]))
    pd += evaluate(metric, v, v)^2
    for (i, v1) in enumerate(fiducials)
        pd += evaluate(metric, v1, v)^2
        for j in i+1:length(fiducials)
            @inbounds pd += evaluate(metric, v1, fiducials[j])^2
        end
    end
    return sum(pd)
end

"""
    comp_Ek2!(ϵ_ball,u_k,Y,v,NNidxs,T,K,metric) → E²(T)
Returns the approximated conditional variance for a specific point in state space
`ns` (index value) with its `K`-nearest neighbors, which indices are stored in
`NNidxs`, for a time horizon `T`. This corresponds to Eqs. 13 & 14 in [^Uzal2011].
The norm specified in input parameter `metric` is used for distance computations.
"""
function comp_Ek2!(ϵ_ball, u_k, Y, ns::Int, NNidxs, T::Int, K::Int, metric)
    # determine neighborhood `T` time steps ahead
    ϵ_ball[1, :] .= Y[ns+T]
    @inbounds for (i, j) in enumerate(NNidxs)
        ϵ_ball[i+1, :] .= Y[j + T]
    end

    # compute center of mass
    @inbounds for i in 1:size(Y)[2]; u_k[i] = sum(view(ϵ_ball, :, i))/(K+1); end # Eq. 14

    E²_sum = 0
    @inbounds for j = 1:K+1
        E²_sum += (evaluate(metric, view(ϵ_ball, j, :), u_k))^2
    end
    E² = E²_sum / (K+1)         # Eq. 13
end

"""
    uzal_cost_local(Y::Dataset; kwargs...) → L_local
Compute the local L-statistic `L_local` for input dataset `Y` according to
Uzal et al.[^Uzal2011]. The length of `L_local` is `length(Y)-Tw` and
denotes a value of the local cost-function to each of the points of the
state space trajectory.

In contrast to [`uzal_cost`](@ref) `σ²` here does not get averaged over all the
state space reference points on the attractor. Therefore, the mean of 'L_local'
is different to `L`, when calling `uzal_cost`, since the averaging is performed
before logarithmizing.

Keywords as in [`uzal_cost`](@ref).

[^Uzal2011]: Uzal, L. C., Grinblat, G. L., Verdes, P. F. (2011). [Optimal reconstruction of dynamical systems: A noise amplification approach. Physical Review E 84, 016223](https://doi.org/10.1103/PhysRevE.84.016223).
"""
function uzal_cost_local(Y::AbstractDataset{D, ET};
        Tw::Int = 40, K::Int = 3, w::Int = 1, samplesize::Real = 0.5,
        metric = Euclidean()
    ) where {D, ET}

    # select a random state space vector sample according to input samplesize
    NN = length(Y)-Tw;
    NNN = floor(Int, samplesize*NN)
    ns = sample(1:NN, NNN; replace=false) # the fiducial point indices

    vs = Y[ns] # the fiducial points in the data set

    vtree = KDTree(Y[1:end-Tw], metric)
    allNNidxs, allNNdist = all_neighbors(vtree, vs, ns, K, w)
    ϵ² = zeros(NNN)             # neighborhood size
    E²_avrg = zeros(NNN)        # averaged conditional variance
    E² = zeros(Tw)
    ϵ_ball = zeros(ET, K+1, D) # preallocation
    u_k = zeros(ET, D)

    # loop over each fiducial point
    for (i,v) in enumerate(vs)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v
        # pairwise distance of fiducial points and `v`
        pdsqrd = fiducial_pairwise_dist_sqrd(view(Y.data, NNidxs), v, metric)
        ϵ²[i] = (2/(K*(K+1))) * pdsqrd  # Eq. 16
        # loop over the different time horizons
        for T = 1:Tw
            E²[T] = comp_Ek2!(ϵ_ball, u_k, Y, ns[i], NNidxs, T, K, metric) # Eqs. 13 & 14
        end
        # Average E²[T] over all prediction horizons
        E²_avrg[i] = mean(E²)                   # Eq. 15
    end
    σ² = E²_avrg ./ ϵ² # noise amplification σ², Eq. 17
    α² = 1 / mean(ϵ².^(-1)) # for normalization, Eq. 21
    L_local = log10.(sqrt.(σ²).*sqrt(α²))
    return L_local
end


text:
struct SRKParam <: EoSParam
    a::PairParam{Float64}
    b::PairParam{Float64}
    acentricfactor::SingleParam{Float64}
    Tc::SingleParam{Float64}
    Pc::SingleParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type SRKModel <: ABCubicModel end
@newmodel SRK SRKModel SRKParam

export SRK
function SRK(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","SAFT/PCSAFT/PCSAFT_unlike.csv"]; userlocations=userlocations, verbose=verbose)
    Mw = params["Mw"]
    k  = params["k"]
    _pc = params["pc"]
    pc = _pc.values

    Tc = params["Tc"]
    Tc_ = Tc.values
    acentricfactor = params["w"]
    a = epsilon_LorentzBerthelot(SingleParam(params["pc"], @. 1/(9*(2^(1/3)-1))*R̄^2*Tc_^2/pc), k)
    b = sigma_LorentzBerthelot(SingleParam(params["pc"], @. (2^(1/3)-1)/3*R̄*Tc_/pc))

    packagedparams = SRKParam(a, b, acentricfactor, Tc,_pc,Mw)
    model = SRK(packagedparams, idealmodel; ideal_userlocations=ideal_userlocations, verbose=verbose)
    return model
end

function ab_consts(::Type{<:SRKModel})
    Ωa =  1/(9*(2^(1/3)-1))
    Ωb = (2^(1/3)-1)/3
    return Ωa,Ωb
end

function cubic_ab(model::SRKModel,T,z=SA[1.0],n=sum(z))
    a = model.params.a.values
    b = model.params.b.values
    ω = model.params.acentricfactor.values
    Tc = model.params.Tc.values
    _1 = one(T+n)
    invn = one(n)/n
    #root of α
    αx = @. (_1+(0.480+1.547*ω-0.176*ω^2)*(1-√(T/Tc))) * z * invn
    #αx = (αi*xi for (αi,xi) in zip(α,x)) #lazy iterator
    #αx .*= x
    #āᾱ = sum(a .* .√(α * α') .* (x * x'))
    āᾱ =dot(αx, Symmetric(a), αx)
    b̄ = dot(z, Symmetric(b), z) * invn * invn
    return āᾱ ,b̄
end

function cubic_abp(model::SRKModel, V, T, x)
    a,b = cubic_ab(model,T,x)
    p =  R̄*T/(V-b) - a/((V+b)*V)
    return a,b,p
end

function cubic_poly(model::SRKModel,p,T,z)
    x = z/sum(z)
    a,b = cubic_ab(model,T,x)
    RT⁻¹ = 1/(R̄*T)
    A = a*p* RT⁻¹* RT⁻¹
    B = b*p* RT⁻¹
    _1 = one(a)
    return [-A*B, -B*(B+_1) + A, -_1, _1]
end

function a_res(model::SRKModel, V, T, z=SA[1.0])
    n = sum(z)
    ā,b̄ = cubic_ab(model,T,z,n)
    ρ = n/V
    RT⁻¹ = 1/(R̄*T)

    return -log(1-b̄*ρ) - ā*RT⁻¹*log(b̄*ρ+1)/b̄
    #+ log(V)
    #=
    x = z/sum(z)
    n = sum(z)
    a = model.params.a.values
    b = model.params.b.values
    ω = model.params.acentricfactor.values
    Tc = model.params.Tc.values

    α = @. (1+(0.480+1.547*ω-0.176*ω^2)*(1-√(T/Tc)))^2

    āᾱ = sum(a .* .√(α * α') .* (x * x'))
    b̄ = sum(b .* (x * x'))
    return -log(V-n*b̄) - āᾱ/(R̄*T*b̄)*log(1+n*b̄/V)
    =#
end

cubic_zc(::SRKModel) = 1/3


text:
using Documenter, AutoMLPipeline

makedocs(
  source = "src",
  build = "build",
  modules = [AutoMLPipeline],
  clean = true,
  sitename = "AutoMLPipeline Documentation",
  doctest = false,
  pages = Any[
    "HOME" => "index.md",
    "Tutorial" => Any[
       "tutorial/pipeline.md",
       "tutorial/preprocessing.md",
       "tutorial/learning.md",
       "tutorial/extending.md"
    ],
    "Manual" => Any[
      "Pipeline" => "man/pipeline.md",
      "Preprocessors" => "man/preprocessors.md",
      "Learners" => "man/learners.md",
      "Meta-Ensembles" => "man/metaensembles.md"
    ],
    "Library" => Any[
      "Types and Functions" => "lib/typesfunctions.md"
    ]
  ],
  format = Documenter.HTML(
     prettyurls = get(ENV, "CI", nothing) == "true"
  )
)

deploydocs(
    repo   = "github.com/IBM/AutoMLPipeline.jl.git",
)


text:
# matrix_eqn.jl

"""
Solves the discrete lyapunov equation.

The problem is given by

    AXA' - X + B = 0

`X` is computed by using a doubling algorithm. In particular, we iterate to
convergence on `X_j` with the following recursions for j = 1, 2,...
starting from X_0 = B, a_0 = A:

    a_j = a_{j-1} a_{j-1}
    X_j = X_{j-1} + a_{j-1} X_{j-1} a_{j-1}'

##### Arguments

- `A::Matrix{Float64}` : An n x n matrix as described above.  We assume in order
for  convergence that the eigenvalues of `A` have moduli bounded by unity
- `B::Matrix{Float64}` :  An n x n matrix as described above.  We assume in order
for convergence that the eigenvalues of `B` have moduli bounded by unity
- `max_it::Int(50)` :  Maximum number of iterations

##### Returns

- `gamma1::Matrix{Float64}` Represents the value X

"""
function solve_discrete_lyapunov(A::ScalarOrArray,
                                 B::ScalarOrArray,
                                 max_it::Int=50)
    # TODO: Implement Bartels-Stewardt
    n = size(A, 2)
    alpha0 = reshape([A;], n, n)
    gamma0 = reshape([B;], n, n)

    alpha1 = zeros(alpha0)
    gamma1 = zeros(gamma0)

    diff = 5
    n_its = 1

    while diff > 1e-15

        alpha1 = alpha0*alpha0
        gamma1 = gamma0 + alpha0*gamma0*alpha0'

        diff = maximum(abs(gamma1 - gamma0))
        alpha0 = alpha1
        gamma0 = gamma1

        n_its += 1

        if n_its > max_it
            error("Exceeded maximum iterations, check input matrices")
        end
    end

    return gamma1
end

"""
Solves the discrete-time algebraic Riccati equation

The prolem is defined as

    X = A'XA - (N + B'XA)'(B'XB + R)^{-1}(N + B'XA) + Q

via a modified structured doubling algorithm.  An explanation of the algorithm
can be found in the reference below.

##### Arguments

- `A` : k x k array.
- `B` : k x n array
- `R` : n x n, should be symmetric and positive definite
- `Q` : k x k, should be symmetric and non-negative definite
- `N::Matrix{Float64}(zeros(size(R, 1), size(Q, 1)))` : n x k array
- `tolerance::Float64(1e-10)` Tolerance level for convergence
- `max_iter::Int(50)` : The maximum number of iterations allowed

Note that `A, B, R, Q` can either be real (i.e. k, n = 1) or matrices.

##### Returns
- `X::Matrix{Float64}` The fixed point of the Riccati equation; a  k x k array
representing the approximate solution

##### References

Chiang, Chun-Yueh, Hung-Yuan Fan, and Wen-Wei Lin. "STRUCTURED DOUBLING
ALGORITHM FOR DISCRETE-TIME ALGEBRAIC RICCATI EQUATIONS WITH SINGULAR CONTROL
WEIGHTING MATRICES." Taiwanese Journal of Mathematics 14, no. 3A (2010): pp-935.

"""
function solve_discrete_riccati(A::ScalarOrArray, B::ScalarOrArray,
                                Q::ScalarOrArray,
                                R::ScalarOrArray,
                                N::ScalarOrArray=zeros(size(R, 1), size(Q, 1));
                                tolerance::Float64=1e-10,
                                max_it::Int=50)
    # Set up
    dist = tolerance + 1
    best_gamma = 0.0

    n = size(R, 1)
    k = size(Q, 1)
    I = eye(k)

    current_min = Inf
    candidates = [0.0, 0.01, 0.1, 0.25, 0.5, 1.0, 2.0, 10.0, 100.0, 10e5]
    BB = B' * B
    BTA = B' * A

    for gamma in candidates
        Z = getZ(R, gamma, BB)
        cn = cond(Z)
        if isfinite(cn)
            Q_tilde = -Q + N' * (Z \ (N + gamma .* BTA)) + gamma .* I
            G0 = B * (Z \ B')
            A0 = (I - gamma .* G0) * A - B * (Z \ N)
            H0 = gamma .* (A' * A0) - Q_tilde
            f1 = cond(Z, Inf)
            f2 = gamma .* f1
            f3 = cond(I + G0 * H0)
            f_gamma = max(f1, f2, f3)

            if f_gamma < current_min
                best_gamma = gamma
                current_min = f_gamma
            end
        end
    end

    if isinf(current_min)
        msg = "Unable to initialize routine due to ill conditioned args"
        error(msg)
    end

    gamma = best_gamma
    R_hat = R + gamma .* BB

    # Initial conditions
    Q_tilde = - Q + N' * (R_hat\(N + gamma .* BTA)) + gamma .* I
    G0 = B * (R_hat\B')
    A0 = (I - gamma .* G0) * A - B * (R_hat\N)
    H0 = gamma .* A'*A0 - Q_tilde
    i = 1

    # Main loop
    while dist > tolerance

        if i > max_it
            msg = "Maximum Iterations reached $i"
            error(msg)
        end

        A1 = A0 * ((I + G0 * H0)\A0)
        G1 = G0 + A0 * G0 * ((I + H0 * G0)\A0')
        H1 = H0 + A0' * ((I + H0*G0)\(H0*A0))

        dist = Base.maxabs(H1 - H0)
        A0 = A1
        G0 = G1
        H0 = H1
        i += 1
    end

    return H0 + gamma .* I  # Return X
end

"""
Simple method to return an element Z in the Riccati equation solver whose type is Float64 (to be accepted by the cond() function)

##### Arguments

- `BB::Float64` : result of B' * B
- `gamma::Float64` : parameter in the Riccati equation solver
- `R::Float64`

##### Returns
- `::Float64` : element Z in the Riccati equation solver

"""
getZ(R::Float64, gamma::Float64, BB::Float64) = R + gamma * BB

"""
Simple method to return an element Z in the Riccati equation solver whose type is Float64 (to be accepted by the cond() function)

##### Arguments

- `BB::Union{Vector, Matrix}` : result of B' * B
- `gamma::Float64` : parameter in the Riccati equation solver
- `R::Float64`

##### Returns
- `::Float64` : element Z in the Riccati equation solver

"""
getZ(R::Float64, gamma::Float64, BB::Union{Vector, Matrix}) = R + gamma * BB[1]

"""
Simple method to return an element Z in the Riccati equation solver whose type is Matrix (to be accepted by the cond() function)

##### Arguments

- `BB::Matrix` : result of B' * B
- `gamma::Float64` : parameter in the Riccati equation solver
- `R::Matrix`

##### Returns
- `::Matrix` : element Z in the Riccati equation solver

"""
getZ(R::Matrix, gamma::Float64, BB::Matrix) = R + gamma .* BB


text:
module UMVETests

using LinearAlgebra
using Subsystems
using Test

@testset "umve estimator" begin
    A = Matrix{Float64}(I, 2, 2)
    B = reshape([1., 0.], :, 1)
    C = I
    Aij = [0.1 0; 0 0]
    W = Diagonal([0.1, 0.1])
    V = W

    @testset "handles empty neighbourhoods" begin
        s = LdtiSubsystem(1, A, B, C)
        @testset "with svd" begin
            o = UMVE(s.A, s.B, s.C, s.E, W, V)
            Rk = Subsystems.R(o.A, o.C, o.Pi, o.W, o.V)
            M∅ = Subsystems.M(o.G, o.C, Rk)
            Kk = Subsystems.K(o.A, o.C, o.Pi, o.W, Rk)
            Ā∅ = Subsystems.barA(Kk, o.C, o.G, M∅)
            L̄∅ = Subsystems.barL(Kk, o.C, o.G, M∅)

            @test M∅ == zeros(0,2)
            @test Ā∅ == (I - Kk*o.C)
            @test L̄∅ == Kk
        end
        @testset "without svd" begin
            o = UMVE(s.A, s.B, s.C, s.E, W, V, false)
            Rk = Subsystems.R(o.A, o.C, o.Pi, o.W, o.V)
            M∅ = Subsystems.M(o.G, o.C, Rk)
            Kk = Subsystems.K(o.A, o.C, o.Pi, o.W, Rk)
            Ā∅ = Subsystems.barA(Kk, o.C, o.G, M∅)
            L̄∅ = Subsystems.barL(Kk, o.C, o.G, M∅)

            @test M∅ == zeros(0,2)
            @test Ā∅ == (I - Kk*o.C)
            @test L̄∅ == Kk
        end
    end

    @testset "Constructor" begin
        @testset "builds correctly if exists" begin
            N = 4
            LSS = [LdtiSubsystem(i, A, B, C) for i in 1:N]
            # make ring topology
            [addNeighbour(LSS[i], LSS[i%N + 1], Aij) for i in 1:N]

            @testset "with svd" begin
                O = [UMVE(s.A, s.B, s.C, s.E, W, V) for s in LSS]
                for o in O
                    @test o.G ≈ reshape([1. 0], :, 1)
                end
            end
            @testset "without svd" begin
                O = [UMVE(s.A, s.B, s.C, s.E, W, V, false) for s in LSS]
                for o in O
                    @test o.G ≈ Matrix{eltype(o.G)}(I, o.nxi, o.nxi)
                end
            end
        end
        @testset "can build with no neighbours" begin
            s = LdtiSubsystem(1, A, B, C)
            @testset "with svd" begin
                o = UMVE(s.A, s.B, s.C, s.E, W, V)
                @test o.G == zeros(s.nx, 0)
            end
            @testset "without svd" begin
                o = UMVE(s.A, s.B, s.C, s.E, W, V, false)
                @test o.G == zeros(s.nx, 0)
                @test o.Ebar == zeros(s.nx, 0)
            end
        end
    end

end

end