text:
"""
    write_opf(file, opf)

Write an MTG with explicit geometry to disk as an OPF file.

# Notes

Node attributes `:ref_meshes` and `:geometry` are treated as reserved keywords and
should not be used without knowing their meaning:

- `:ref_meshes`: a `RefMeshes` structure that holds the MTG reference meshes.
- `:geometry`: a [`geometry`](@ref) instance

# Examples

```julia
using PlantGeom
file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf")
opf = read_opf(file)
write_opf("test.opf", opf)
opf2 = read_opf("test.opf")
viz(opf2)
```
"""
function write_opf(file, mtg)
    doc = XMLDocument()
    opf_elm = ElementNode("opf")
    setroot!(doc, opf_elm)
    opf_elm["version"] = 2.0
    opf_elm["editable"] = true

    # Writing the reference meshes (meshBDD):
    meshBDD = addelement!(opf_elm, "meshBDD")

    if mtg[:ref_meshes] === nothing
        error("No reference meshes found in the MTG.")
    end

    for (key, mesh) in enumerate(mtg[:ref_meshes].meshes)
        mesh_elm = addelement!(meshBDD, "mesh")
        mesh_elm["name"] = mesh.name
        mesh_elm["shape"] = ""
        mesh_elm["Id"] = key - 1 # opf uses 0-based indexing
        mesh_elm["enableScale"] = mesh.taper

        addelement!(
            mesh_elm,
            "points",
            string("\n", join([vcat([p.coords for p in mesh.mesh.points]...)...], "\t"), "\n")
        )

        if length(mesh.normals) == Meshes.nelements(mesh) && length(mesh.normals) != Meshes.nvertices(mesh)
            # If the normals are per triangle, re-compute them per vertex:
            vertex_normals = normals_vertex(mesh)
        else
            vertex_normals = mesh.normals
        end

        norm_elm = addelement!(
            mesh_elm,
            "normals",
            string("\n", join(vcat([p.coords for p in vertex_normals]...), "\t"), "\n")
        )

        if mesh.texture_coords !== nothing && length(mesh.texture_coords) > 0
            # texture_coords are optional
            norm_elm = addelement!(
                mesh_elm,
                "textureCoords",
                string("\n", join(vcat([[p.coords...] for p in mesh.texture_coords]...), "\t"), "\n")
            )
        end

        faces_elm = addelement!(mesh_elm, "faces")

        face_id = [0]
        for i = 1:Meshes.nelements(Meshes.topology(mesh.mesh))
            face_elm = addelement!(
                faces_elm,
                "face",
                string("\n", join(Meshes.topology(mesh.mesh).connec[i].indices .- 1, "\t"), "\n")
            )
            #? NB: we remove one because face index are 0-based in the opf
            face_elm["Id"] = face_id[1]
            face_id[1] = face_id[1] + 1
        end

    end

    # Parsing the materialBDD section.
    materialBDD = addelement!(opf_elm, "materialBDD")
    for (key, mesh) in enumerate(mtg[:ref_meshes].meshes)
        mat_elm = addelement!(materialBDD, "material")
        mat_elm["Id"] = key - 1 # opf uses 0-based indexing

        mat = material_to_opf_string(mesh.material)
        addelement!(
            mat_elm,
            "emission",
            mat[:emission]
        )

        addelement!(
            mat_elm,
            "ambient",
            mat[:ambient]
        )

        addelement!(
            mat_elm,
            "diffuse",
            mat[:diffuse]
        )

        addelement!(
            mat_elm,
            "specular",
            mat[:specular]
        )

        addelement!(
            mat_elm,
            "shininess",
            mat[:shininess]
        )
    end

    # Parsing the shapeBDD section.
    shapeBDD = addelement!(opf_elm, "shapeBDD")
    for (key, mesh) in enumerate(mtg[:ref_meshes].meshes)
        shape_elm = addelement!(shapeBDD, "shape")
        shape_elm["Id"] = key - 1 # opf uses 0-based indexing

        addelement!(
            shape_elm,
            "name",
            mesh.name
        )

        addelement!(
            shape_elm,
            "meshIndex",
            string(key - 1)
        )

        addelement!(
            shape_elm,
            "materialIndex",
            string(key - 1)
        )
    end

    # Parsing the attributeBDD section:
    attrBDD = addelement!(opf_elm, "attributeBDD")
    attrs = MultiScaleTreeGraph.get_features(mtg)
    for i = 1:size(attrs, 1)
        shape_elm = addelement!(attrBDD, "attribute")
        shape_elm["name"] = string(attrs[i, 1])
        attr_type = attrs[i, 2]

        if attr_type == "STRING"
            attr_type_opf = "String"
        elseif attr_type == "REAL"
            attr_type_opf = "Double"
        elseif attr_type == "INT"
            attr_type_opf = "Integer"
        end

        shape_elm["class"] = attr_type_opf
    end

    # Parsing the topology section:
    mtg_topology_to_xml!(mtg, opf_elm)

    # prettyprint(doc)
    write(file, doc)
end

"""
    mtg_to_opf_link(link)

Takes an MTG link as input ("/", "<" or "+") and outputs its corresponding link as declared
in the OPF format ("decomp", "follow" or "branch")
"""
function mtg_to_opf_link(link)
    if link == "/"
        "decomp"
    elseif link == "<"
        "follow"
    elseif link == "+"
        "branch"
    end
end

"""
    mtg_topology_to_xml!(node, xml_parent)

Write the MTG topology, attributes and geometry into XML format. This function is used to
write the "topology" section of the OPF.
"""
function mtg_topology_to_xml!(node, xml_parent, xml_gtparent=nothing, ref_meshes=get_ref_meshes(node))

    if isroot(node)
        xml_parent = attributes_to_xml(node, xml_parent, xml_gtparent, ref_meshes)
    end

    if !isleaf(node)
        for chnode in MultiScaleTreeGraph.ordered_children(node)
            xml_node = attributes_to_xml(chnode, xml_parent, xml_gtparent, ref_meshes)
            mtg_topology_to_xml!(chnode, xml_node, xml_parent, ref_meshes)
        end
    end
end

"""
    attributes_to_xml(node, xml_parent)

Write an MTG node into an XML node.
"""
function attributes_to_xml(node, xml_parent, xml_gtparent, ref_meshes)
    opf_link = isroot(node) ? "topology" : mtg_to_opf_link(node.MTG.link)

    if opf_link == "follow"
        xml_gtparent === nothing ? error("Root node should start with a '/' link") : nothing
        xml_node = addelement!(xml_gtparent, opf_link)
    else
        xml_node = addelement!(xml_parent, opf_link)
    end

    xml_node["class"] = node.MTG.symbol
    xml_node["scale"] = node.MTG.scale
    xml_node["id"] = node.id #! maybe this should be `node.MTG.index` instead ? But I think is is unique

    for key in keys(node.attributes)
        if key == :geometry
            geom = addelement!(xml_node, string(key))
            geom["class"] = "Mesh"

            if node[key].ref_mesh_index === nothing
                get_ref_mesh_index!(node, ref_meshes)
            end
            addelement!(geom, "shapeIndex", string(node[key].ref_mesh_index - 1))
            # NB: opf uses 0-based indexing, that's why we use ref_mesh_index - 1

            # Make the homogeneous matrix from the transformations:

            mat4x4 = mat_to_opf_string(node[key].transformation)

            addelement!(
                geom,
                "mat",
                string(
                    "\n",
                    join(mat4x4[1, :], "\t"),
                    "\n",
                    join(mat4x4[2, :], "\t"),
                    "\n",
                    join(mat4x4[3, :], "\t"),
                    "\n"
                )
            )
            #? NB: Only the three first rows are written as the fourth is always the same
            addelement!(geom, "dUp", string(node[key].dUp))
            addelement!(geom, "dDwn", string(node[key].dDwn))
        elseif key == :ref_meshes
            # We don't write the reference meshes here but before in the opf
            continue
        else
            addelement!(xml_node, string(key), string(node[key]))
        end
    end

    return xml_node
end

function mat_to_opf_string(trans::T) where {T<:UniformScaling}
    Matrix{Float64}(I, 3, 4)
end

function mat_to_opf_string(trans::IdentityTransformation)
    Matrix{Float64}(I, 3, 4)
end

function mat_to_opf_string(trans)
    hcat(trans.linear, trans.translation)
end


text:
using Documenter, OTSM

makedocs(
    format = Documenter.HTML(),
    sitename = "OTSM.jl",
    authors = "Hua Zhou, Joong-Ho Won",
    clean = true,
    debug = true,
    pages = [
        "index.md"
    ]
)

deploydocs(
    repo   = "github.com/Hua-Zhou/OTSM.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)

ext: jl


text:
export DuckietownEnv, DuckietownLF, DuckietownNav

##
#Wrapper to control the simulator using velocity and steering angle
#instead of differential drive motor velocities
##

struct DuckietownEnv
    sim::Simulator
    gain::Float32
    trim::Float32
    radius::Float32
    k::Float32
    limit::Float32
end

function DuckietownEnv(; gain::Float32=1f0, trim::Float32=0f0, radius::Float32=0.0318f0,
                       k::Float32=27f0, limit::Float32=1f0, kwargs...)
    sim = Simulator(; kwargs...)
    #logger.info('using DuckietownEnv')
    return DuckietownEnv(sim, gain, trim, radius, k, limit)
end

reset!(dt_env::DuckietownEnv) = reset!(dt_env.sim)

function step!(dt_env::DuckietownEnv, action::Vector{Float32})
     vel, angle = action

     # Distance between the wheels
     baseline = dt_env.sim.fixedparams.wheel_dist

     # assuming same motor constants k for both motors
     k_r = dt_env.k
     k_l = dt_env.k

     # adjusting k by gain and trim
     k_r_inv = (dt_env.gain + dt_env.trim) / k_r
     k_l_inv = (dt_env.gain - dt_env.trim) / k_l

     omega_r = (vel + 0.5f0 * angle * baseline) / dt_env.radius
     omega_l = (vel - 0.5f0 * angle * baseline) / dt_env.radius

     # conversion from motor rotation rate to duty cycle
     u_r = omega_r * k_r_inv
     u_l = omega_l * k_l_inv

     # limiting output to limit, which is 1.0 for the duckiebot
     u_r_limited = max(min(u_r, dt_env.limit), -dt_env.limit)
     u_l_limited = max(min(u_l, dt_env.limit), -dt_env.limit)

     vels = [u_l_limited, u_r_limited]

     obs, v, reward, done, info = step!(dt_env.sim, vels)
     mine = Dict{String,Any}()
     mine["k"] = dt_env.k
     mine["gain"] = dt_env.gain
     mine["train"] = dt_env.trim
     mine["radius"] = dt_env.radius
     mine["omega_r"] = omega_r
     mine["omega_l"] = omega_l
     info["DuckietownEnv"] = mine
     return obs, action, reward, done, info
end

##
#Environment for the Duckietown lane following task with
#and without obstacles (LF and LFV tasks)
##
struct DuckietownLF
    dt_env::DuckietownEnv
end

DuckietownLF(; kwargs...) = DuckietownLF(DuckietownEnv(; kwargs...))

reset!(dt_lf::DuckietownLF) = reset!(dt_lf.dt_env)

step!(dt_lf::DuckietownLF, action::Vector{Float32}) = step!(dt_lf.dt_env, action)

##
#Environment for the Duckietown navigation task (NAV)
##
struct DuckietownNav
    dt_env::DuckietownEnv
    goal_tile::Union{Nothing, Dict{String,Any}}
end

function DuckietownNav(; kwargs...)
    goal_tile = nothing
    DuckietownNav(DuckietownEnv(; kwargs...), goal_tile)
end

function reset!(dt_nav::DuckietownNav)
    sim = dt_nav.dt_env.sim
    obs = reset!(sim)

    # Find the tile the agent starts on
    start_tile_pos = get_grid_coords(sim.cur_pos)
    start_tile = self._get_tile(start_tile_pos...)

    # Select a random goal tile to navigate to
    @assert length(_drivable_tiles(sim)) > 1
    while true
        tile_idx = rand(sim.fixedparams.rng, 1:length(_drivable_tiles(sim)))
        self.goal_tile = _drivable_tiles(sim)[tile_idx]
        dt_nav.goal_tile != sim.fixedparams._map.start_tile && break
    end

    return obs
end

function step!(dt_nav::DuckietownNav, action::Vector{Float32})
    obs, a, reward, done, info = step!(dt_env, action)

    info["goal_tile"] = dt_nav.goal_tile

    # TODO: add term to reward based on distance to goal?

    cur_tile_coords = get_grid_coords(dt_nav.dt_env.sim.cur_pos)
    cur_tile = _get_tile(cur_tile_coords)

    if cur_tile == dt_nav.goal_tile
        done = true
        reward = 1000f0
    end

    return obs, a, reward, done, info
end


text:
module GooglePlaces

include("loginitialiser.jl")

include("functools.jl")

include("api.jl")

include("nearbysearch.jl")

include("placedetails.jl")

end # module GooglePlaces


text:
using RTRRT
using Test

init_pose = [0.5, 2.0, 0]
params = RTRRT.parameters(1.0, [0.0, 10.0], [0.0, 10.0], 5000, 0.4, 0.1, 1.0, 20)

tree = RTRRT.Tree(init_pose, params)
for i = 1:1
    RTRRT.expansion(tree)
    # RTRRT.rewire(tree)
    RTRRT.rewire_from_root(tree)
end
RTRRT.plot_tree(tree.all_nodes)
RTRRT.visualize_map(tree.map)
@testset "RTRRT.jl" begin
    # Write your tests here
end


text:
module Test_linsystools

using LinearAlgebra
using MatrixPencils
using Test
using Polynomials

@testset "Linear System Tools" begin

@testset "lseval, lpseval & lps2ls" begin

# g = [z]
A = [0.0 1.0; -1.0 0.0];
E = [0.0 0.0; 0.0 -1.0];
B = [-1.0; 0.0];
C = [1.0 0.0];
D = [0.0];

@time  G = lseval(A,E,B,C,D,1)
@test G ≈ ones(1,1)

@time  G = lseval(A,E,B,C,D,Inf)
@test G ≈ Inf*ones(1,1)

# g = [z 1/(z-1)]
A = [ 1     0     0; 0     1     0; 0     0     1];
E = [ 0     1     0; 0     0     0; 0     0     1];
B = [ 0     0; -1     0; 0     1];
C = [1     0     1];
D = [0 0];
@time  G = lseval(A,E,B,C,D,1)
@test G ≈ [1 Inf]

@time  G = lseval(A,E,B,C,D,Inf)
@test G ≈ [Inf 0]

@time  G = lseval(A,E,B,C,D,0)
@test G ≈ [0 -1]

# g = [1/z (z-1); z z/(1-z)]
A =  [    0     0     0     0     0     0
0     1     0     0     0     0
0     0     1     0     0     0
0     0     0     1     0     0
0     0     0     0     1     0
0     0     0     0     0     1];
E = [1     0     0     0     0     0
0     0     1     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     1];
B = [ 1     0
0     0
-1     0
0     0
0    -2
0     1];
C = [    1.0000         0         0    0.5000   -0.5000         0
0    1.0000         0         0         0   -1.0000];
D = [     0     0
0    -1];

@time  G = lseval(A,E,B,C,D,0)
@test ≈(G,[Inf -1; 0 0],atol=1.e-7)

@time  G = lseval(A,E,B,C,D,1)
@test ≈(G,[1 0; 1 Inf],atol=1.e-7)

@time  G = lseval(A,E,B,C,D,2)
@test ≈(G,[0.5   1.0; 2.0  -2.0],atol=1.e-7)

@time  G = lseval(A,E,B,C,D,Inf)
@test ≈(G,[  0.0  Inf; Inf   -1.0],atol=1.e-7)

@time  G = lseval(A,E,B,C,D,im)
@test ≈(G,[ -1.0im  -1.0+1.0im; 1.0im  -0.5+0.5im],atol=1.e-7)

@time  G = lseval(A,E,B,C,D,exp(im))
@test G ≈ [5.403023058681398e-01 - 8.414709848078965e-01im -4.596976941318602e-01 + 8.414709848078965e-01im;
           5.403023058681398e-01 + 8.414709848078965e-01im -5.000000000000000e-01 + 9.152438608562260e-01im]

A = rand(2,2); E = rand(2,2); B = rand(2,3); C = rand(1,2); D = rand(1,3);
@time  G = lseval(A,E,B,C,D,Inf)
@test G ≈ D

A2 = [
    -2    -3     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0
     0     0    -2    -3     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     1     0     0     0
     0     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     0     1];
 E2 = [
     1     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     1     0     0     0     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     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0];
 B2 = [
     1     0
     0     0
     0     1
     0     0
    -1     0
     0     0
     0    -1
     0     0
     0     0];
C2 = [
     1     0     1    -3     0     1     0     2     0
     0     1     1     3     0     1     0     0     1];
D2 = zeros(Int,2,2);
F2 = ones(size(B2)...);
G2 = ones(size(C2)...);
H2 = ones(size(D2)...);
sys = (A2,E2,B2,F2,C2,G2,D2,H2);

sys1 = lps2ls(A2,I,B2,missing,C2,G2,D2,H2);
sys2 = lps2ls(A2,I,B2,missing,C2,G2,D2,H2,compacted = true);
@test lsequal(sys1...,sys2...,atol1=1.e-7,atol2=1.e-7)

sys1 = lps2ls(sys...);
sys2 = lps2ls(sys...,compacted = true);
@test lsequal(sys1...,sys2...,atol1=1.e-7,atol2=1.e-7)

@test lseval(sys2...,1) ≈ (C2-G2)*((E2-A2)\(B2-F2)) + D2-H2

g1 = lseval(sys1...,1)
@test g1 ≈ (C2-G2)*((E2-A2)\(B2-F2)) + D2-H2
g2 = lseval(sys2...,1)
@test g2 ≈ (C2-G2)*((E2-A2)\(B2-F2)) + D2-H2

sys = (A2,E2,B2,F2,C2,G2,D2,H2);
@time  G = lpseval(sys...,0)
@test G ≈ lseval(A2,E2,B2,C2,D2,0)
@test G ≈ lseval(lps2ls(sys...)...,0)

@time  G = lpseval(A2,E2,view(B2,:,1),view(F2,:,1),C2,G2,view(D2,:,1),view(H2,:,1),0)
@test G ≈ lseval(A2,E2,view(B2,:,1),C2,view(D2,:,1),0)

@time  G = lpseval(A2,E2,B2,F2,view(C2,1:1,:),view(G2,1:1,:),view(D2,1:1,:),view(H2,1:1,:),0)
@test G ≈ lseval(A2,E2,B2,view(C2,1:1,:),view(D2,1:1,:),0)

@time  G = lpseval(sys...,Inf)
@test G ≈ [Inf Inf; Inf Inf]

@time  G = lpseval(A2,E2,view(B2,:,1),view(F2,:,1),C2,G2,view(D2,:,1),view(H2,:,1),Inf)
@test G ≈ [Inf; Inf]

# Example Van Dooren & Dewilde, LAA 1983.
# strongly controllable realization

sys = ([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.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; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[1.0 3.0 0.0; 1.0 4.0 2.0; 0.0 -1.0 -2.0],
[-1.0 -4.0 -2.0; -0.0 -0.0 -0.0; -1.0 -4.0 -2.0],
[1.0 2.0 -2.0; 0.0 -1.0 -2.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]);
@time  G = lpseval(sys...,0)
@test G ≈ [1.0 2.0 -2.0; 0.0 -1.0 -2.0; 0.0 0.0 0.0]

@time  G = lpseval(sys...,Inf)
@test G ≈ [Inf  Inf  Inf;
Inf  Inf  Inf;
Inf  Inf  Inf]

end

@testset "lsequal and lpsequal" begin

A = zeros(0,0); E = zeros(0,0); C = zeros(2,0); B = zeros(0,3); D = zeros(2,3);
@test lsequal(A,E,B,C,D,A,E,B,C,D)

@test !lsequal(A,E,B,C,D,A,E,C,B,D)

A = zeros(0,0); E = zeros(0,0); C = zeros(2,0); G = zeros(2,0);  B = zeros(0,3); F = zeros(0,3); D = zeros(2,3); H = zeros(2,3);
@test lpsequal(A,E,B,F,C,G,D,H,A,E,B,F,C,G,D,H)

@test !lpsequal(A,E,B,F,C,G,D,H,A,E,C,F,B,G,D,H)

end

@testset "lsminreal and lsminreal2" begin

for fast in (true, false)

A2 = zeros(0,0); E2 = zeros(0,0); C2 = zeros(0,0); B2 = zeros(0,0); D2 = zeros(0,0);
sys = (A2,E2,B2,C2,D2);

sys1 = lsminreal(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,0,0)

sys1 = lsminreal2(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,0,0)

A2 = rand(3,3); E2 = zeros(3,3); C2 = zeros(0,3); B2 = zeros(3,0); D2 = zeros(0,0);
sys = (A2,E2,B2,C2,D2);

sys1 = lsminreal(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (3,0,0)

sys1 = lsminreal(sys..., fast = fast, contr = false)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,3,0)

sys1 = lsminreal2(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (3,0,0)

sys1 = lsminreal2(sys..., fast = fast, contr = false)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,3,0)

# B and D vectors
A2 = rand(3,3); E2 = zeros(3,3); C2 = zeros(1,3); B2 = zeros(3,1); D2 = zeros(1);
sys = (A2,E2,B2,C2,D2);

sys1 = lsminreal(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (3,0,0)

sys1 = lsminreal(sys..., fast = fast, contr = false)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,3,0)

sys1 = lsminreal2(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (3,0,0)

sys1 = lsminreal2(sys..., fast = fast, contr = false)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,3,0)

# Example 1: DeTeran, Dopico, Mackey, ELA 2009

A2 = [1.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  1.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  1.0  0.0
0.0  0.0  0.0  0.0  0.0  1.0];

E2 = [0.0  0.0  1.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  1.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.0  0.0  0.0  0.0  0.0  0.0];

B2 = [0.0  0.0
0.0  0.0
0.0  1.0
1.0  0.0
1.0  0.0
0.0  0.0];

C2 = [-1.0   0.0  0.0  0.0  0.0  0.0
0.0  -1.0  0.0  0.0  0.0  0.0];

D2 = [0.0  0.0
0.0  1.0];

sys = (A2,E2,B2,C2,D2);

# compute minimal realization
sys1 = lsminreal(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,1)  # nuc == 2 && nuo == 0 && nse == 1
# an order reduction without enforcing controllability and observability may not be possible
sys1 = lsminreal(sys...,contr=false,obs=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,0,0)  # nuc == 0 && nuo == 0 && nse == 0 !!
# compute an irreducible realization which still contains a non-dynamic mode
sys1 = lsminreal(sys...,noseig=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,0)  # nuc == 2 && nuo == 0 && nse == 0

sys = (E2,A2,B2,C2,D2);
# compute minimal realization for a standard system (i.e., irreducible realization)
sys1 = lsminreal(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,0)  # nuc == 2 && nuo == 0 && nse == 0
# an order reduction without enforcing controllability and observability may not be possible
sys1 = lsminreal(sys...,contr=false,obs=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,0,0)  # nuc == 0 && nuo == 0 && nse == 0
# compute an irreducible realization which still contains a non-dynamic mode
sys1 = lsminreal(sys...,noseig=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,0)  # nuc == 2 && nuo == 0 && nse == 0

sys = (A2,E2,B2,C2,D2);
# compute minimal realization
sys1 = lsminreal2(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,1)  # nuc == 2 && nuo == 0 && nse == 1
# minimal realization is possible only applying the infinite controllability/observability algorithm
sys1 = lsminreal2(sys...,finite = false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,1)  # nuc == 2 && nuo == 0 && nse == 0
# order reduction may results even when applying the finite controllability/observability algorithm
sys1 = lsminreal2(sys...,infinite = false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (1,0,0)  # nuc == 1 && nuo == 0 && nse == 0
# an order reduction without enforcing controllability and observability may not be possible
sys1 = lsminreal2(sys...,contr=false,obs=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,0,0)  # nuc == 0 && nuo == 0 && nse == 0 !!
# compute an irreducible realization which still contains a non-dynamic mode
sys1 = lsminreal2(sys...,noseig=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,0)  # nuc == 2 && nuo == 0 && nse == 0

sys = (E2,A2,B2,C2,D2);
# compute minimal realization for a standard system (i.e., irreducible realization)
sys1 = lsminreal2(sys..., fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,0)  # nuc == 2 && nuo == 0 && nse == 0
# an order reduction without enforcing controllability and observability may not be possible
sys1 = lsminreal2(sys...,contr=false,obs=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (0,0,0)  # nuc == 0 && nuo == 0 && nse == 0
# compute an irreducible realization which still contains a non-dynamic mode
sys1 = lsminreal2(sys...,noseig=false, fast = fast)
@test lsequal(sys...,sys1[1:5]...) && sys1[6:8] == (2,0,0)  # nuc == 2 && nuo == 0 && nse == 0

# Example Van Dooren & Dewilde, LAA 1983.
# P = zeros(3,3,3)
# P[:,:,1] = [1 2 -2; 0 -1 -2; 0 0 0]
# P[:,:,2] = [1 3 0; 1 4 2; 0 -1 -2]
# P[:,:,3] = [1 4 2; 0 0 0; 1 4 2]

# observable realization with A2 = I
A2 = [ 1.0  0.0  0.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.0
0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0];

E2 = [ 0.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0
0.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.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  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0];

B2 = [ 0.0   0.0   0.0
0.0   0.0   0.0
0.0   0.0   0.0
1.0   3.0   0.0
1.0   4.0   2.0
0.0  -1.0  -2.0
1.0   4.0   2.0
0.0   0.0   0.0
1.0   4.0   2.0];

C2 = [ -1.0   0.0   0.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.0
0.0   0.0  -1.0  0.0  0.0  0.0  0.0  0.0  0.0];

D2 = [ 1.0   2.0  -2.0
0.0  -1.0  -2.0
0.0   0.0   0.0];

# build a strong (least order) structured linearization which preserves the finite eigenvalues,
# infinite zeros (with multiplicities), the left and right Kronecker structure
#A2,E2,B2,C2,D2 = pol2lp(P)
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 5 && nuo == 0 && nse == 1
val, iz, info = spzeros(A1,E1,B1,C1,D1)
@test val ≈ [1.0] && iz == [] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)
val, info = speigvals(A1,E1,B1,C1,D1)
@test val ≈ [1., Inf, Inf, Inf] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,obs=false, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 5 && nuo == 0 && nse == 1
val, iz, info = spzeros(A1,E1,B1,C1,D1)
@test val ≈ [1.0] && iz == [] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)
val, info = speigvals(A1,E1,B1,C1,D1)
@test val ≈ [1., Inf, Inf, Inf] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 5 && nuo == 0 && nse == 1
val, iz, info = spzeros(A1,E1,B1,C1,D1)
@test val ≈ [1.0] && iz == [] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)
val, info = speigvals(A1,E1,B1,C1,D1)
@test val ≈ [1., Inf, Inf, Inf] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,obs=false,finite=false, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 5 && nuo == 0 && nse == 1
val, iz, info = spzeros(A1,E1,B1,C1,D1)
@test val ≈ [1.0] && iz == [] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)
val, info = speigvals(A1,E1,B1,C1,D1)
@test val ≈ [1., Inf, Inf, Inf] &&  (info.rki, info.lki, info.id, info.nf) == ([0], [1], [1, 1, 1], 1)

# perform lsminreal with A and E interchanged
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
E1, A1, B1, C1, D1, nuc, nuo, nse  = lsminreal(E,A,B,C,D,obs = false, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 5 && nuo == 0 && nse == 0

# perform lsminreal2 with A and E interchanged
A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
E1, A1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(E,A,B,C,D,obs = false, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 5 && nuo == 0 && nse == 0

# controllable realization with A2 = I
A2 = [ 1.0  0.0  0.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.0
0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0];

E2 = [ 0.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0  0.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.0  0.0  0.0  1.0  0.0
0.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.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  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0];

B2 = [   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.0   0.0   0.0
0.0   0.0   0.0
-1.0   0.0   0.0
0.0  -1.0   0.0
0.0   0.0  -1.0];

C2 = [  1.0  4.0  2.0  1.0   3.0   0.0  0.0  0.0  0.0
0.0  0.0  0.0  1.0   4.0   2.0  0.0  0.0  0.0
1.0  4.0  2.0  0.0  -1.0  -2.0  0.0  0.0  0.0
];

D2 = [ 1.0   2.0  -2.0
0.0  -1.0  -2.0
0.0   0.0   0.0];

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,contr = false, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,contr = false,finite=false, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A2 = [
    1.     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     0     0     0
    0     0     1     0     0     0     0     0     0
    0     0     0     1     0     0     0     0     0
    0     0     0     0     1     0     0     0     0
    0     0     0     0     0     1     0     0     0
    0     0     0     0     0     0     1     0     0
    0     0     0     0     0     0     0     1     0
    0     0     0     0     0     0     0     0     1];
E2 = [
    0.     0     0     0     0     0     0     0     0
    1     0     0     0     0     0     0     0     0
    0     1     0     0     0     0     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     1     0     0     0     0
    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     1     0];
B2 = [
    -1.     0     0
     0     0     0
     0     0     0
     0    -1     0
     0     0     0
     0     0     0
     0     0    -1
     0     0     0
     0     0     0];
C2 = [
    0.     1     1     0     3     4     0     0     2
    0     1     0     0     4     0     0     2     0
    0     0     1     0    -1     4     0    -2     2];

D2 = zeros(Float64,3,3);

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  =
     lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  =
     lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,contr=false, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,I,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,I,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 6 && nuo == 3 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  =
      lsminreal2(A,E,B,C,D,contr=false,finite=false,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast);
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 5 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,I,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,I,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 6 && nuo == 3 && nse == 0

A2 = [
    -2    -3     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0
     0     0    -2    -3     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     1     0     0     0
     0     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     0     1];
 E2 = [
     1     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     1     0     0     0     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     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0];
 B2 = [
     1     0
     0     0
     0     1
     0     0
    -1     0
     0     0
     0    -1
     0     0
     0     0];
 C2 = [
     1     0     1    -3     0     1     0     2     0
     0     1     1     3     0     1     0     0     1];
D2 = zeros(Int,2,2);

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,noseig=true, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 2 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,obs=false,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 0 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,contr=false,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 2 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,noseig=true, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 2 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,obs=false,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 0 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,contr=false,finite=false,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 2 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') &&
      nuc == 2 && nuo == 0 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') &&
      nuc == 2 && nuo == 0 && nse == 0

# Example 1 - (Varga, Kybernetika, 1990)
A2 = [
1 0 0 0 -1 0 0 0
0 1 0 0 0 -1 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
0 0 0 0 -1 0 0 0
0 0 0 0 0 -1 0 0
0 0 0 0 3 0 1 0
0 0 0 0 0 2 0 1
]
E2 = [
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 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0
0 0 0 1 0 0 0 0
]
B2 = [
      -1 1
      0 0
      0 0
      0 0
      1 -2
      -2 3
      0 0
      3 -3
]
C2 = [
      0 0 0 0 0 0 -1 0
      0 0 0 0 0 0 0 -1
]
D2 = zeros(Int,2,2);

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 6 && nuo == 0 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time @time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 6 && nuo == 0 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') &&
      nuc == 2 && nuo == 4 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A',E',C',B',D',atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') &&
      nuc == 2 && nuo == 4 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A',E',C',B',D',contr=false, atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') &&
      nuc == 0 && nuo == 6 && nse == 1

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A',E',C',B',D',contr=false, atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1',E1',C1',B1',D1') &&
      nuc == 0 && nuo == 6 && nse == 1

# B and D vectors
A2 = [
      0 0 0 -24 0 0 0 0 0 0 0
      1 0 0 -50 0 0 0 0 0 0 0
      0 1 0 -35 0 0 0 0 0 0 0
      0 0 1 -10 0 0 0 0 0 0 0
      0 0 0 0 0 0 0 -30 0 0 0
      0 0 0 0 1 0 0 -61 0 0 0
      0 0 0 0 0 1 0 -41 0 0 0
      0 0 0 0 0 0 1 -11 0 0 0
      0 0 0 0 0 0 0 0 0 0 -15
      0 0 0 0 0 0 0 0 1 0 -23
      0 0 0 0 0 0 0 0 0 1 -9
]
E2 = I;
#B2 = reshape([18; 42; 30; 6; 10;17;8;1;0;-10;-2;],11,1)
B2 = [18; 42; 30; 6; 10;17;8;1;0;-10;-2;]
C2 = [0 0 0 0 0 0 0 1 0 0 0]
D2 = [0]

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 7 && nuo == 3 && nse == 0

@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A1,E1,B1,C1,D1,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 0 && nse == 0

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 7 && nuo == 3 && nse == 0

@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A1,E1,B1,C1,D1,atol1 = 1.e-7, atol2 = 1.e-7, fast = fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 0 && nse == 0

for Ty in (Float64, Complex{Float64})

#fast = true; Ty = Complex{Float64}
n = 10; m = 5; p = 6;
rE = 5;
rA22 = min(2,n-rE)
n2 = rE+rA22
n3 = n-n2
A2 = [rand(Ty,rE,n);
rand(Ty,rA22,n2) zeros(Ty,rA22,n3);
rand(Ty,n3,rE) zeros(Ty,n3,n-rE) ];
E2 = [rand(Ty,rE,rE) zeros(Ty,rE,n-rE);
zeros(Ty,n-rE,n)]
Q = qr(rand(Ty,n,n)).Q
Z = qr(rand(Ty,n,n)).Q
A2 = Q*A2*Z; E2 = Q*E2*Z;
B2 = rand(Ty,n,m); C2 = rand(Ty,p,n);
D2 = zeros(Ty,p,m)

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,fast=fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 0 && nse == rA22

A = copy(A2); E = copy(E2); C = copy(C2); B = copy(B2); D = copy(D2);
@time A1, E1, B1, C1, D1, nuc, nuo, nse  = lsminreal2(A,E,B,C,D,atol1 = 1.e-7, atol2 = 1.e-7,fast=fast)
@test lsequal(A2,E2,B2,C2,D2,A1,E1,B1,C1,D1) &&
      nuc == 0 && nuo == 0 && nse == rA22

end

end
end

@testset "lpsminreal" begin

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,0); F2 = zeros(0,0); C2 = zeros(0,0); G2 = zeros(0,0); D2 = zeros(0,0); H2 = zeros(0,0);
sys = (A2,E2,B2,F2,C2,G2,D2,H2);

@time sys1 = lpsminreal(sys...)
@test lpsequal(sys...,sys1[1:8]...) && sys1[11:12] == (0,0)

A2 = zeros(0,0); E2 = zeros(0,0); B2 = zeros(0,3); F2 = zeros(0,3); C2 = zeros(2,0); G2 = zeros(2,0); D2 = rand(2,3); H2 = rand(2,3);
sys = (A2,E2,B2,F2,C2,G2,D2,H2);

@time sys1 = lpsminreal(sys...)
@test lpsequal(sys...,sys1[1:8]...) && sys1[11:12] == (0,0)

# Example Van Dooren & Dewilde, LAA 1983.
# strongly controllable realization

sys = ([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.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; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[1.0 3.0 0.0; 1.0 4.0 2.0; 0.0 -1.0 -2.0],
[-1.0 -4.0 -2.0; -0.0 -0.0 -0.0; -1.0 -4.0 -2.0],
[1.0 2.0 -2.0; 0.0 -1.0 -2.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]);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (0,2)

# Example Van Dooren & Dewilde, LAA 1983.
# strongly observable realization
sys = ([1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 3.0 0.0; 1.0 4.0 2.0; 0.0 -1.0 -2.0],
[-1.0 -4.0 -2.0; -0.0 -0.0 -0.0; -1.0 -4.0 -2.0],
[0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0],
[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0],
[1.0 2.0 -2.0; 0.0 -1.0 -2.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]);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (2,0)

# Example 1: DeTeran, Dopico, Mackey, ELA 2009
# P = zeros(2,2,3);
# P[:,:,1] = [0 0; 0 1.];
# P[:,:,2] = [0 1.; 1. 0];
# P[:,:,3] = [1. 0; 0 0];

# strongly controllable realization
sys = ([1.0 0.0; 0.0 1.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0], [1.0 0.0; 0.0 1.0],
[0.0 1.0; 1.0 0.0], [-1.0 -0.0; -0.0 -0.0], [0.0 0.0; 0.0 1.0], [0.0 0.0; 0.0 0.0])
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (0, 1)

# strongly observabable realization
sys = ([1.0 0.0; 0.0 1.0], [0.0 0.0; 0.0 0.0], [0.0 1.0; 1.0 0.0], [-1.0 -0.0; -0.0 -0.0],
[0.0 0.0; 0.0 0.0], [1.0 0.0; 0.0 1.0], [0.0 0.0; 0.0 1.0], [0.0 0.0; 0.0 0.0])
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (1,0)

A2 = [
    -2    -3     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0
     0     0    -2    -3     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     1     0     0     0
     0     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     0     1];
 E2 = [
     1     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0
     0     0     0     1     0     0     0     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     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     1     0];
 B2 = [
     1     0
     0     0
     0     1
     0     0
    -1     0
     0     0
     0    -1
     0     0
     0     0];
C2 = [
     1     0     1    -3     0     1     0     2     0
     0     1     1     3     0     1     0     0     1];
D2 = zeros(Int,2,2);
F2 = zeros(size(B2)...);
G2 = zeros(size(C2)...);
H2 = zeros(size(D2)...);

sys = (A2,E2,B2,F2,C2,G2,D2,H2);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (2,4)

sys = (A2,E2,B2,F2,C2,G2,D2,H2);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...,contr=false);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (0,4)

sys = (A2,E2,B2,F2,C2,G2,D2,H2);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...,obs=false);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (2,0)

A2 = [
      0 0 0 -24 0 0 0 0 0 0 0
      1 0 0 -50 0 0 0 0 0 0 0
      0 1 0 -35 0 0 0 0 0 0 0
      0 0 1 -10 0 0 0 0 0 0 0
      0 0 0 0 0 0 0 -30 0 0 0
      0 0 0 0 1 0 0 -61 0 0 0
      0 0 0 0 0 1 0 -41 0 0 0
      0 0 0 0 0 0 1 -11 0 0 0
      0 0 0 0 0 0 0 0 0 0 -15
      0 0 0 0 0 0 0 0 1 0 -23
      0 0 0 0 0 0 0 0 0 1 -9
]
E2 = Matrix{Float64}(I,11,11);
B2 = reshape([18; 42; 30; 6; 10;17;8;1;0;-10;-2;],11,1)
C2 = reshape([0 0 0 0 0 0 0 1 0 0 0],1,11)
D2 = zeros(Int,1,1)
F2 = zeros(size(B2)...);
G2 = zeros(size(C2)...);
H2 = zeros(size(D2)...);
sys = (A2,E2,B2,F2,C2,G2,D2,H2);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (7, 3)

fast = true; Ty = Float64
for fast in (true, false)

for Ty in (Float64, Complex{Float64})

#fast = true; Ty = Complex{Float64}
n = 10; m = 5; p = 6;
rE = 5;
rA22 = min(2,n-rE)
n2 = rE+rA22
n3 = n-n2
A2 = [rand(Ty,rE,n);
rand(Ty,rA22,n2) zeros(Ty,rA22,n3);
rand(Ty,n3,rE) zeros(Ty,n3,n-rE) ];
E2 = [rand(Ty,rE,rE) zeros(Ty,rE,n-rE);
zeros(Ty,n-rE,n)]
Q = qr(rand(Ty,n,n)).Q
Z = qr(rand(Ty,n,n)).Q
A2 = Q*A2*Z; E2 = Q*E2*Z;
B2 = rand(Ty,n,m); C2 = rand(Ty,p,n);
D2 = zeros(Ty,p,m)

F2 = zeros(Ty,size(B2)...);
G2 = zeros(Ty,size(C2)...);
H2 = rand(Ty,size(D2)...);
sys = (A2,E2,B2,F2,C2,G2,D2,H2);
@time A1, E1, B1, F1, C1, G1, D1, H1, V, W, nuc, nuo  = lpsminreal(sys...,atol1 = 1.e-7, atol2 = 1.e-7,fast=fast);
@test lpsequal(sys..., A1, E1, B1/W, F1/W, V'\C1, V'\G1, V'\D1/W, V'\H1/W) && (nuc, nuo) == (5, 3)

end

end

end

end # regular testset
end # module


text:
# Target block layout in multi-field computations
#
# u [c][1 , 1][1 ,bu][np,1 ,nu]
# v [c][1 ,  ][bv,  ][np,nv,  ]
# u [f][lc][1 ,lc][1 ,bu][np,1 ,nu]
# v [f][lc][lc,  ][bv,  ][np,nv,  ]
# u⁺[f][1 ,lc][1 ,bu][np,1 ,nu]
# v⁺[f][lc,  ][bv,  ][np,nv,  ]
# u [c][lf][1 ,lF][1 ,bu][np,1 ,nu]
# v [c][lf][lF,  ][bv,  ][np,nv,  ]
#
# λ [f][1 , 1][1 ,bλ][np,1 ,nλ]
# μ [f][1 ,  ][bμ,  ][np,nμ,  ]
# λ [c][lf][1 ,lf][1 ,bλ][np,1 ,nλ]
# μ [c][lf][lf,  ][bμ,  ][np,nμ,  ]
#
# uh [c][np]
# uh [f][lc][np]
# uh⁺[f][np]
# uh [c][lf][np]

struct ArrayBlock{A,N}
  array::Array{A,N}
  touched::Array{Bool,N}
  function ArrayBlock(array::Array{A,N},touched::Array{Bool,N}) where {A,N}
    @check size(array) == size(touched)
    new{A,N}(array,touched)
  end
end

const VectorBlock = ArrayBlock{A,1} where A
const MatrixBlock = ArrayBlock{A,2} where A

Base.size(b::ArrayBlock) = size(b.array)
Base.length(b::ArrayBlock) = length(b.array)
Base.eltype(::Type{<:ArrayBlock{A}}) where A = A
Base.eltype(b::ArrayBlock{A}) where A = A
Base.ndims(b::ArrayBlock{A,N}) where {A,N} = N
Base.ndims(::Type{ArrayBlock{A,N}}) where {A,N} = N
function Base.getindex(b::ArrayBlock,i...)
  if !b.touched[i...]
    return nothing
  end
  b.array[i...]
end
function Base.setindex!(b::ArrayBlock,v,i...)
 @check b.touched[i...] "Only touched entries can be set"
 b.array[i...] = v
end
function Base.show(io::IO,o::ArrayBlock)
  print(io,"ArrayBlock($(o.array), $(o.touched))")
end
function Base.show(io::IO,k::MIME"text/plain",o::ArrayBlock)
  function _sub_block_str(a)
    function _nice_typeof(v)
      "$v"
    end
    function _nice_typeof(::Type{<:VectorBlock{A}}) where A
      "VectorBlock{$(_nice_typeof(A))}"
    end
    function _nice_typeof(::Type{<:MatrixBlock{A}}) where A
      "MatrixBlock{$(_nice_typeof(A))}"
    end
    if isa(a,AbstractArray) || isa(a,ArrayBlock)
      n = replace(_nice_typeof(typeof(a)),"Gridap.Fields."=>"")
      if ndims(a) == 1
        "$(length(a))-element $n"
      else
        s = replace(replace(replace("$(size(a))",", "=>"x"),"("=>""),")"=>"")
        "$s $n"
      end
    else
      "$a"
    end
  end
  println(io,"$(_sub_block_str(o)):")
  lis = LinearIndices(o.array)
  for i in CartesianIndices(o.array)
    if ndims(o) == 1
      s = "[$(lis[i])]"
    else
      s = replace(replace("$(Tuple(i))","("=>"["),")"=>"]")
    end
    if lis[i] != 1
      print(io,"\n")
    end
    print(io," $s = $(_sub_block_str(o[i]))")
  end
end

function Arrays.testitem(f::ArrayBlock{A}) where A
  @notimplementedif !isconcretetype(A)
  i = findall(f.touched)
  if length(i) != 0
    f.array[i[1]]
  else
    testvalue(A)
  end
end

function Arrays.testvalue(::Type{ArrayBlock{A,N}}) where {A,N}
  s = ntuple(i->0,Val(N))
  array = Array{A,N}(undef,s)
  touched = Array{Bool,N}(undef,s)
  ArrayBlock(array,touched)
end

function Arrays.CachedArray(a::ArrayBlock)
  ai = testitem(a)
  ci = CachedArray(ai)
  array = Array{typeof(ci),ndims(a)}(undef,size(a))
  for i in eachindex(a.array)
    if a.touched[i]
      array[i] = CachedArray(a.array[i])
    end
  end
  ArrayBlock(array,a.touched)
end

function unwrap_cached_array(a::CachedArray)
  a.array
end

function unwrap_cached_array(a::ArrayBlock)
  cache = return_cache(unwrap_cached_array,a)
  evaluate!(cache,unwrap_cached_array,a)
end

function return_cache(::typeof(unwrap_cached_array),a::ArrayBlock)
  ai = testitem(a)
  ci = return_cache(unwrap_cached_array,ai)
  ri = evaluate!(ci,unwrap_cached_array,ai)
  c = Array{typeof(ci),ndims(a)}(undef,size(a))
  array = Array{typeof(ri),ndims(a)}(undef,size(a))
  for i in eachindex(a.array)
    if a.touched[i]
      c[i] = return_cache(unwrap_cached_array,a.array[i])
    end
  end
  ArrayBlock(array,a.touched), c
end

function evaluate!(cache,::typeof(unwrap_cached_array),a::ArrayBlock)
  r,c = cache
  for i in eachindex(a.array)
    if a.touched[i]
      r.array[i] = evaluate!(c[i],unwrap_cached_array,a.array[i])
    end
  end
  r
end

#LinearAlgebra.promote_leaf_eltypes(a::ArrayBlock) = LinearAlgebra.promote_leaf_eltypes(a.array)

function Base.:≈(a::AbstractArray{<:ArrayBlock},b::AbstractArray{<:ArrayBlock})
  all(z->z[1]≈z[2],zip(a,b))
end

function Base.:≈(a::ArrayBlock,b::ArrayBlock)
  if size(a) != size(b)
    return false
  end
  if a.touched != b.touched
    return false
  end
  for i in eachindex(a.array)
    if a.touched[i]
      if !(a.array[i] ≈ b.array[i])
        return false
      end
    end
  end
  true
end

function Base.:(==)(a::ArrayBlock,b::ArrayBlock)
  if size(a) != size(b)
    return false
  end
  for i in eachindex(a.array)
    if a.touched[i] && b.touched[i]
      if a.array[i] != b.array[i]
        return false
      end
    elseif a.touched[i]
      return false
    elseif b.touched[i]
      return false
    end
  end
  true
end

Base.copy(a::ArrayBlock) = ArrayBlock(copy(a.array),copy(a.touched))
Base.eachindex(a::ArrayBlock) = eachindex(a.array)

struct ZeroBlockMap <: Map end

function return_cache(::ZeroBlockMap,a::AbstractArray,b::AbstractArray)
  CachedArray(similar(a,eltype(a),size(b)))
end

function evaluate!(cache,::ZeroBlockMap,a,b::AbstractArray)
  setsize!(cache,size(b))
  r = cache.array
  fill!(r,zero(eltype(r)))
  r
end

function return_cache(::ZeroBlockMap,a::ArrayBlock,b::ArrayBlock)
  A = eltype(a)
  N = ndims(b)
  array = Array{A,N}(undef,size(b))
  touched = fill(false,size(b))
  ArrayBlock(array,touched)
end

function evaluate!(cache,::ZeroBlockMap,a,b::ArrayBlock)
  @check size(cache) == size(b)
  cache
end

struct BlockMap{N} <: Map
  size::NTuple{N,Int}
  indices::Vector{CartesianIndex{N}}
end

function BlockMap(l::Integer,i::Integer)
  s = (l,)
  ci = CartesianIndex((i,))
  BlockMap(s,[ci])
end

function BlockMap(l::Integer,inds::Vector{<:Integer})
  s = (l,)
  cis = [CartesianIndex((i,)) for i in inds]
  BlockMap(s,cis)
end

function BlockMap(s::NTuple,i::Integer)
  cis = CartesianIndices(s)
  ci = cis[i]
  BlockMap(s,[ci])
end

function BlockMap(s::NTuple,inds::Vector{<:NTuple})
  cis = [CartesianIndex(i) for i in inds]
  BlockMap(s,cis)
end

function return_cache(k::BlockMap{N},a::A...) where {A,N}
  array = Array{A,N}(undef,k.size)
  touched = fill(false,k.size)
  for (t,i) in enumerate(k.indices)
    array[i] = a[t]
    touched[i] = true
  end
  ArrayBlock(array,touched)
end

function evaluate!(cache,k::BlockMap{N},a::A...) where {A,N}
  @check size(cache) == k.size
  for (t,i) in enumerate(k.indices)
    cache.array[i] = a[t]
  end
  cache
end

function lazy_map(k::BlockMap,a::LazyArray{<:Fill{<:Broadcasting{typeof(∘)}}})
  b = a.args[1]
  ϕ = a.args[2]
  c = lazy_map(k,b)
  lazy_map(Broadcasting(∘),c,ϕ)
end

function lazy_map(k::BlockMap,a::MemoArray)
  b = a.parent
  c = lazy_map(k,b)
  MemoArray(c)
end

function lazy_map(k::BlockMap,a::LazyArray{<:Fill{typeof(transpose)}})
  @notimplementedif k.size[1] != 1
  lis = LinearIndices(k.size)
  k2 = BlockMap(k.size[2],[ lis[i] for i in k.indices])
  b = a.args[1]
  c = lazy_map(k2,b)
  lazy_map(transpose,c)
end

function lazy_map(::typeof(evaluate),a::LazyArray{<:Fill{<:BlockMap}},x::AbstractArray)
  args = map(i->lazy_map(evaluate,i,x),a.args)
  k = a.maps.value
  lazy_map(k,args...)
end

# This lazy_map is triggered from function gradient(a::CellField) with optimization
# purposes. See https://github.com/gridap/Gridap.jl/pull/638 for more details.
function lazy_map(k::Broadcasting{typeof(gradient)},a::LazyArray{<:Fill{<:BlockMap}})
  args = map(i->lazy_map(k,i),a.args)
  bm = a.maps.value
  lazy_map(bm,args...)
end

function return_cache(f::ArrayBlock{A,N},x) where {A,N}
  fi = testitem(f)
  li = return_cache(fi,x)
  fix = evaluate!(li,fi,x)
  l = Array{typeof(li),N}(undef,size(f.array))
  g = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(f.array[i],x)
    end
  end
  ArrayBlock(g,f.touched),l
end

function evaluate!(cache,f::ArrayBlock,x)
  g,l = cache
  @check g.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      g.array[i] = evaluate!(l[i],f.array[i],x)
    end
  end
  g
end

function linear_combination(u::ArrayBlock,f::ArrayBlock)
  i::Int = findfirst(f.touched)
  fi = f.array[i]
  ui = u.array[i]
  ufi = linear_combination(ui,fi)
  g = Vector{typeof(ufi)}(undef,length(f.touched))
  for i in eachindex(f.touched)
    if f.touched[i]
      g[i] = linear_combination(u.array[i],f.array[i])
    end
  end
  ArrayBlock(g,f.touched)
end

function return_cache(k::LinearCombinationMap,u::ArrayBlock,fx::ArrayBlock)
  i::Int = findfirst(fx.touched)
  fxi = fx.array[i]
  ui = u.array[i]
  li = return_cache(k,ui,fxi)
  ufxi = evaluate!(li,k,ui,fxi)
  l = Vector{typeof(li)}(undef,size(fx.array))
  g = Vector{typeof(ufxi)}(undef,size(fx.array))
  for i in eachindex(fx.array)
    if fx.touched[i]
      l[i] = return_cache(k,u.array[i],fx.array[i])
    end
  end
  ArrayBlock(g,fx.touched),l
end

function evaluate!(cache,k::LinearCombinationMap,u::ArrayBlock,fx::ArrayBlock)
  g,l = cache
  @check g.touched == fx.touched
  for i in eachindex(fx.array)
    if fx.touched[i]
      g.array[i] = evaluate!(l[i],k,u.array[i],fx.array[i])
    end
  end
  g
end

function Base.transpose(f::ArrayBlock{A,1} where A)
  fi = testitem(f)
  fit = transpose(fi)
  g = Matrix{typeof(fit)}(undef,(1,length(f.touched)))
  for i in eachindex(f.touched)
    if f.touched[i]
      g[i] = transpose(f.array[i])
    end
  end
  ArrayBlock(g,collect(transpose(f.touched)))
end

function Base.transpose(f::ArrayBlock{A,2} where A)
  fi = testitem(f)
  fit = transpose(fi)
  ni,nj = size(f)
  g = Matrix{typeof(fit)}(undef,(nj,ni))
  for i in 1:ni
    for j in 1:nj
      if f.touched[i,j]
        g[j,i] = transpose(f.array[i,j])
      end
    end
  end
  ArrayBlock(g,collect(transpose(f.touched)))
end

function  return_cache(k::TransposeMap,f::ArrayBlock{A,1} where A)
  fi = testitem(f)
  li = return_cache(k,fi)
  fix = evaluate!(li,k,fi)
  l = Matrix{typeof(li)}(undef,(1,length(f.array)))
  g = Matrix{typeof(fix)}(undef,(1,length(f.array)))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,f.array[i])
    end
  end
  ArrayBlock(g,collect(transpose(f.touched))),l
end

function  evaluate!(cache,k::TransposeMap,f::ArrayBlock{A,1} where A)
  g,l = cache
  @check g.touched == transpose(f.touched)
  for i in eachindex(f.array)
    if f.touched[i]
      g.array[i] = evaluate!(l[i],k,f.array[i])
    end
  end
  g
end

function integrate(f::ArrayBlock{A,N} where A,args...) where N
  fi = testitem(f)
  intfi = integrate(fi,args...)
  g = Array{typeof(intfi),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      g[i] = integrate(f.array[i],args...)
    end
  end
  ArrayBlock(g,f.touched)
end

function return_cache(
  k::IntegrationMap,fx::ArrayBlock{A,N} where A,args...) where N
  fxi = testitem(fx)
  li = return_cache(k,fxi,args...)
  ufxi = evaluate!(li,k,fxi,args...)
  l = Array{typeof(li),N}(undef,size(fx.array))
  g = Array{typeof(ufxi),N}(undef,size(fx.array))
  for i in eachindex(fx.array)
    if fx.touched[i]
      l[i] = return_cache(k,fx.array[i],args...)
    end
  end
  ArrayBlock(g,fx.touched),l
end

function evaluate!(cache,k::IntegrationMap,fx::ArrayBlock,args...)
  g,l = cache
  @check g.touched == fx.touched
  for i in eachindex(fx.array)
    if fx.touched[i]
      g.array[i] = evaluate!(l[i],k,fx.array[i],args...)
    end
  end
  g
end

function return_value(k::Broadcasting,f::ArrayBlock)
  evaluate(k,f)
end

function return_cache(k::Broadcasting,f::ArrayBlock{A,N}) where {A,N}
  fi = testitem(f)
  li = return_cache(k,fi)
  fix = evaluate!(li,k,fi)
  l = Array{typeof(li),N}(undef,size(f.array))
  g = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,f.array[i])
    end
  end
  ArrayBlock(g,f.touched),l
end

function evaluate!(cache,k::Broadcasting,f::ArrayBlock)
  g,l = cache
  @check g.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      g.array[i] = evaluate!(l[i],k,f.array[i])
    end
  end
  g
end

function return_value(k::Broadcasting{typeof(∘)},f::ArrayBlock,h::Field)
  evaluate(k,f,h)
end

function return_cache(k::Broadcasting{typeof(∘)},f::ArrayBlock{A,N},h::Field) where {A,N}
  fi = testitem(f)
  li = return_cache(k,fi,h)
  fix = evaluate!(li,k,fi,h)
  l = Array{typeof(li),N}(undef,size(f.array))
  g = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,f.array[i],h)
    end
  end
  ArrayBlock(g,f.touched),l
end

function evaluate!(cache,k::Broadcasting{typeof(∘)},f::ArrayBlock,h::Field)
  g,l = cache
  @check g.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      g.array[i] = evaluate!(l[i],k,f.array[i],h)
    end
  end
  g
end

function return_value(k::Broadcasting{<:Operation},f::ArrayBlock,h::Field)
  evaluate(k,f,h)
end

function return_cache(k::Broadcasting{<:Operation},f::ArrayBlock{A,N},h::Field) where {A,N}
  fi = testitem(f)
  li = return_cache(k,fi,h)
  fix = evaluate!(li,k,fi,h)
  l = Array{typeof(li),N}(undef,size(f.array))
  g = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,f.array[i],h)
    end
  end
  ArrayBlock(g,f.touched),l
end

function evaluate!(cache,k::Broadcasting{<:Operation},f::ArrayBlock,h::Field)
  g,l = cache
  @check g.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      g.array[i] = evaluate!(l[i],k,f.array[i],h)
    end
  end
  g
end

function return_value(k::Broadcasting{<:Operation},h::Field,f::ArrayBlock)
  evaluate(k,h,f)
end

function return_cache(k::Broadcasting{<:Operation},h::Field,f::ArrayBlock{A,N}) where {A,N}
  fi = testitem(f)
  li = return_cache(k,h,fi)
  fix = evaluate!(li,k,h,fi)
  l = Array{typeof(li),N}(undef,size(f.array))
  g = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,h,f.array[i])
    end
  end
  ArrayBlock(g,f.touched),l
end

function evaluate!(cache,k::Broadcasting{<:Operation},h::Field,f::ArrayBlock)
  g,l = cache
  @check g.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      g.array[i] = evaluate!(l[i],k,h,f.array[i])
    end
  end
  g
end

function return_value(k::Broadcasting{<:Operation},h::ArrayBlock,f::ArrayBlock)
  evaluate(k,h,f)
end

function return_cache(k::Broadcasting{<:Operation},h::ArrayBlock,f::ArrayBlock{A,N}) where {A,N}
  @notimplemented
end

function evaluate!(cache,k::Broadcasting{<:Operation},h::ArrayBlock,f::ArrayBlock)
  @notimplemented
end

function return_value(k::BroadcastingFieldOpMap,f::ArrayBlock,g::AbstractArray)
  evaluate(k,f,g)
end

function return_cache(
  k::BroadcastingFieldOpMap,f::ArrayBlock{A,N},g::AbstractArray) where {A,N}
  fi = testitem(f)
  li = return_cache(k,fi,g)
  fix = evaluate!(li,k,fi,g)
  l = Array{typeof(li),N}(undef,size(f.array))
  h = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,f.array[i],g)
    end
  end
  ArrayBlock(h,f.touched),l
end

function evaluate!(cache,k::BroadcastingFieldOpMap,f::ArrayBlock,g::AbstractArray)
  h,l = cache
  @check h.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      h.array[i] = evaluate!(l[i],k,f.array[i],g)
    end
  end
  h
end

function return_value(k::BroadcastingFieldOpMap,g::AbstractArray,f::ArrayBlock)
  evaluate(k,g,f)
end

function return_cache(
  k::BroadcastingFieldOpMap,g::AbstractArray,f::ArrayBlock{A,N}) where {A,N}
  fi = testitem(f)
  li = return_cache(k,g,fi)
  fix = evaluate!(li,k,g,fi)
  l = Array{typeof(li),N}(undef,size(f.array))
  h = Array{typeof(fix),N}(undef,size(f.array))
  for i in eachindex(f.array)
    if f.touched[i]
      l[i] = return_cache(k,g,f.array[i])
    end
  end
  ArrayBlock(h,f.touched),l
end

function evaluate!(cache,k::BroadcastingFieldOpMap,g::AbstractArray,f::ArrayBlock)
  h,l = cache
  @check h.touched == f.touched
  for i in eachindex(f.array)
    if f.touched[i]
      h.array[i] = evaluate!(l[i],k,g,f.array[i])
    end
  end
  h
end

for op in (:+,:-,:*)
  @eval begin

    function return_value(k::Broadcasting{typeof($op)},f::ArrayBlock,g::ArrayBlock)
      evaluate(k,f,g)
    end

    function return_cache(k::Broadcasting{typeof($op)},f::ArrayBlock,g::ArrayBlock)
      return_cache(BroadcastingFieldOpMap($op),f,g)
    end

    function evaluate!(cache,k::Broadcasting{typeof($op)},f::ArrayBlock,g::ArrayBlock)
      evaluate!(cache,BroadcastingFieldOpMap($op),f,g)
    end

  end
end

function return_value(k::Broadcasting{typeof(*)},f::Number,g::ArrayBlock)
  evaluate(k,f,g)
end

function return_cache(k::Broadcasting{typeof(*)},f::Number,g::ArrayBlock)
  gi = testitem(g)
  ci = return_cache(k,f,gi)
  hi = evaluate!(ci,k,f,gi)
  array = Array{typeof(hi),ndims(g.array)}(undef,size(g.array))
  c = Array{typeof(ci),ndims(g.array)}(undef,size(g.array))
  for i in eachindex(g.array)
    if g.touched[i]
      c[i] = return_cache(k,f,g.array[i])
    end
  end
  ArrayBlock(array,g.touched), c
end

function evaluate!(cache,k::Broadcasting{typeof(*)},f::Number,g::ArrayBlock)
  r, c = cache
  for i in eachindex(g.array)
    if g.touched[i]
      r.array[i] = evaluate!(c[i],k,f,g.array[i])
    end
  end
  r
end

function return_value(k::Broadcasting{typeof(*)},f::ArrayBlock,g::Number)
  evaluate(k,f,g)
end

function return_cache(k::Broadcasting{typeof(*)},f::ArrayBlock,g::Number)
  return_cache(k,g,f)
end

function evaluate!(cache,k::Broadcasting{typeof(*)},f::ArrayBlock,g::Number)
  evaluate!(cache,k,g,f)
end

function return_value(k::BroadcastingFieldOpMap,f::ArrayBlock,g::ArrayBlock)
  evaluate(k,f,g)
end

function return_cache(k::BroadcastingFieldOpMap,f::ArrayBlock,g::ArrayBlock)
  @notimplemented
end

function evaluate!(cache,k::BroadcastingFieldOpMap,f::ArrayBlock,g::ArrayBlock)
  @notimplemented
end

function return_cache(
  k::BroadcastingFieldOpMap,f::ArrayBlock{A,N},g::ArrayBlock{B,N}) where {A,B,N}
  @notimplementedif size(f) != size(g)
  fi = testvalue(A)
  gi = testvalue(B)
  ci = return_cache(k,fi,gi)
  hi = evaluate!(ci,k,fi,gi)
  m = ZeroBlockMap()
  a = Array{typeof(hi),N}(undef,size(f.array))
  b = Array{typeof(ci),N}(undef,size(f.array))
  zf = Array{typeof(return_cache(m,fi,gi))}(undef,size(f.array))
  zg = Array{typeof(return_cache(m,gi,fi))}(undef,size(f.array))
  t = map(|,f.touched,g.touched)
  for i in eachindex(f.array)
    if f.touched[i] && g.touched[i]
      b[i] = return_cache(k,f.array[i],g.array[i])
    elseif f.touched[i]
      _fi = f.array[i]
      zg[i] = return_cache(m,gi,_fi)
      _gi = evaluate!(zg[i],m,gi,_fi)
      b[i] = return_cache(k,_fi,_gi)
    elseif g.touched[i]
      _gi = g.array[i]
      zf[i] = return_cache(m,fi,_gi)
      _fi = evaluate!(zf[i],m,fi,_gi)
      b[i] = return_cache(k,_fi,_gi)
    end
  end
  ArrayBlock(a,t), b, zf, zg
end

function evaluate!(
  cache,k::BroadcastingFieldOpMap,f::ArrayBlock{A,N},g::ArrayBlock{B,N}) where {A,B,N}
  a,b, zf, zg = cache
  @check size(f) == size(g)
  @check size(a) == size(g)
  m = ZeroBlockMap()
  for i in eachindex(f.array)
    if f.touched[i] && g.touched[i]
      a.array[i] = evaluate!(b[i],k,f.array[i],g.array[i])
    elseif f.touched[i]
      fi = f.array[i]
      gi = evaluate!(zg[i],m,nothing,fi)
      a.array[i] = evaluate!(b[i],k,fi,gi)
    elseif g.touched[i]
      gi = g.array[i]
      fi = evaluate!(zf[i],m,nothing,gi)
      a.array[i] = evaluate!(b[i],k,fi,gi)
    end
  end
  a
end

function return_cache(
  k::BroadcastingFieldOpMap,f::ArrayBlock{A,1},g::ArrayBlock{B,2}) where {A,B}
  fi = testvalue(A)
  gi = testvalue(B)
  ci = return_cache(k,fi,gi)
  hi = evaluate!(ci,k,fi,gi)
  @check size(g.array,1) == 1 || size(g.array,2) == 0
  s = (size(f.array,1),size(g.array,2))
  a = Array{typeof(hi),2}(undef,s)
  b = Array{typeof(ci),2}(undef,s)
  t = fill(false,s)
  for j in 1:s[2]
    for i in 1:s[1]
      if f.touched[i] && g.touched[1,j]
        t[i,j] = true
        b[i,j] = return_cache(k,f.array[i],g.array[1,j])
      end
    end
  end
  ArrayBlock(a,t), b
end

function return_cache(
  k::BroadcastingFieldOpMap,f::ArrayBlock{A,2},g::ArrayBlock{B,1}) where {A,B}
  fi = testvalue(A)
  gi = testvalue(B)
  ci = return_cache(k,fi,gi)
  hi = evaluate!(ci,k,fi,gi)
  @check size(f.array,1) == 1 || size(f.array,2) == 0
  s = (size(g.array,1),size(f.array,2))
  a = Array{typeof(hi),2}(undef,s)
  b = Array{typeof(ci),2}(undef,s)
  t = fill(false,s)
  for j in 1:s[2]
    for i in 1:s[1]
      if f.touched[1,j] && g.touched[i]
        t[i,j] = true
        b[i,j] = return_cache(k,f.array[1,j],g.array[i])
      end
    end
  end
  ArrayBlock(a,t), b
end

function evaluate!(
  cache,k::BroadcastingFieldOpMap,f::ArrayBlock{A,1},g::ArrayBlock{B,2}) where {A,B}
  a,b = cache
  s = size(a.array)
  for j in 1:s[2]
    for i in 1:s[1]
      if f.touched[i] && g.touched[1,j]
        a.array[i,j] = evaluate!(b[i,j],k,f.array[i],g.array[1,j])
      end
    end
  end
  a
end

function evaluate!(
  cache,k::BroadcastingFieldOpMap,f::ArrayBlock{A,2},g::ArrayBlock{B,1}) where {A,B}
  a,b = cache
  s = size(a.array)
  for j in 1:s[2]
    for i in 1:s[1]
      if f.touched[1,j] && g.touched[i]
        a.array[i,j] = evaluate!(b[i,j],k,f.array[1,j],g.array[i])
      end
    end
  end
  a
end

function return_value(
  k::BroadcastingFieldOpMap,a::(ArrayBlock{A,N} where A)...) where N
  evaluate(k,a...)
end

function return_cache(
  k::BroadcastingFieldOpMap,a::(ArrayBlock{A,N} where A)...) where N
  a1 = first(a)
  @notimplementedif any(ai->size(ai)!=size(a1),a)
  ais = map(ai->testvalue(eltype(ai)),a)
  ci = return_cache(k,ais...)
  bi = evaluate!(ci,k,ais...)
  c = Array{typeof(ci),N}(undef,size(a1))
  array = Array{typeof(bi),N}(undef,size(a1))
  for i in eachindex(a1.array)
    @notimplementedif any(ai->ai.touched[i]!=a1.touched[i],a)
    if a1.touched[i]
      _ais = map(ai->ai.array[i],a)
      c[i] = return_cache(k,_ais...)
    end
  end
  ArrayBlock(array,a1.touched), c
end

function evaluate!(
  cache,k::BroadcastingFieldOpMap,a::(ArrayBlock{A,N} where A)...) where N
  a1 = first(a)
  @notimplementedif any(ai->size(ai)!=size(a1),a)
  r,c = cache
  for i in eachindex(a1.array)
    @notimplementedif any(ai->ai.touched[i]!=a1.touched[i],a)
    if a1.touched[i]
      ais = map(ai->ai.array[i],a)
      r.array[i] = evaluate!(c[i],k,ais...)
    end
  end
  r
end

function return_value(k::BroadcastingFieldOpMap,a::ArrayBlock...)
  evaluate(k,a...)
end

function return_cache(k::BroadcastingFieldOpMap,a::ArrayBlock...)
  @notimplemented
end

function evaluate!(cache,k::BroadcastingFieldOpMap,a::ArrayBlock...)
  @notimplemented
end

function return_value(
  k::BroadcastingFieldOpMap,a::Union{ArrayBlock,AbstractArray}...)
  evaluate(k,a...)
end

function return_cache(
  k::BroadcastingFieldOpMap,a::Union{ArrayBlock,AbstractArray}...)

  function _replace_nz_blocks(a::ArrayBlock,bi::AbstractArray)
    N = ndims(a.array)
    array = Array{typeof(bi),N}(undef,size(a))
    for i in eachindex(a.array)
      if a.touched[i]
        array[i] = bi
      end
    end
    ArrayBlock(array,a.touched)
  end

  function _replace_nz_blocks(a::ArrayBlock,bi::ArrayBlock)
    bi
  end

  inds = findall(ai->isa(ai,ArrayBlock),a)
  @notimplementedif length(inds) == 0
  a1 = a[inds[1]]
  b = map(ai->_replace_nz_blocks(a1,ai),a)
  c = return_cache(k,b...)
  c,b
end

function evaluate!(
  cache,k::BroadcastingFieldOpMap,a::Union{ArrayBlock,AbstractArray}...)

  function _replace_nz_blocks!(a::ArrayBlock,bi::AbstractArray)
    for i in eachindex(a.array)
      if a.touched[i]
        a.array[i] = bi
      end
    end
    a
  end

  function _replace_nz_blocks!(a::ArrayBlock,bi::ArrayBlock)
    bi
  end

  c, b = cache
  d = map((i,j)->_replace_nz_blocks!(i,j),b,a)
  evaluate!(c,k,d...)
end

function Base.:+(a::ArrayBlock,b::ArrayBlock)
  BroadcastingFieldOpMap(+)(a,b)
end

function Base.:-(a::ArrayBlock,b::ArrayBlock)
  BroadcastingFieldOpMap(-)(a,b)
end

function Base.:*(a::Number,b::ArrayBlock)
  bi = testitem(b)
  ci = a*bi
  array = Array{typeof(ci),ndims(b)}(undef,size(b))
  for i in eachindex(b.array)
    if b.touched[i]
      array[i] = a*b.array[i]
    end
  end
  ArrayBlock(array,b.touched)
end

function Base.:*(a::ArrayBlock,b::Number)
  b*a
end

function Base.:*(a::ArrayBlock{A,2},b::ArrayBlock{B,1}) where {A,B}
  @check size(a.array,2) == size(b.array,1)
  ai = testvalue(A)
  bi = testvalue(B)
  ri = ai*bi
  array = Vector{typeof(ri)}(undef,size(a.array,1))
  touched = fill(false,size(a.array,1))
  ni,nj = size(a.array)
  for i in 1:ni
    for j in 1:nj
      if a.touched[i,j] && b.touched[j]
        if !touched[i]
          array[i] = a.array[i,j]*b.array[j]
          touched[i] = true
        else
          array[i] = array[i] + a.array[i,j]*b.array[j]
        end
      end
    end
  end
  ArrayBlock(array,touched)
end

function Base.:*(a::ArrayBlock{A,2},b::ArrayBlock{B,2}) where {A,B}
  @check size(a.array,2) == size(b.array,1)
  ai = testvalue(A)
  bi = testvalue(B)
  ri = ai*bi
  ni = size(a.array,1)
  nj = size(b.array,2)
  nk = size(a.array,2)
  s = (ni,nj)
  array = Matrix{typeof(ri)}(undef,s)
  touched = fill(false,s)
  ni,nj = size(a.array)
  for i in 1:ni
    for j in 1:nj
      for k in 1:nk
        if a.touched[i,k] && b.touched[k,j]
          if !touched[i,j]
            array[i,j] = a.array[i,k]*b.array[k,j]
            touched[i,j] = true
          else
            array[i,j] = array[i,j] + a.array[i,k]*b.array[k,j]
          end
        end
      end
    end
  end
  ArrayBlock(array,touched)
end

function LinearAlgebra.rmul!(a::ArrayBlock,β)
  for i in eachindex(a.touched)
    if a.touched[i]
      rmul!(a.array[i],β)
    end
  end
end

function _zero_entries!(a::AbstractArray)
  fill!(a,zero(eltype(a)))
end

function _zero_entries!(a::ArrayBlock)
  for i in eachindex(a.touched)
    if a.touched[i]
      _zero_entries!(a.array[i])
    end
  end
end

function LinearAlgebra.mul!(c::ArrayBlock,a::ArrayBlock,b::ArrayBlock)
  _zero_entries!(c)
  mul!(c,a,b,1,0)
end

function LinearAlgebra.mul!(
  c::ArrayBlock{C,1} where C,
  a::ArrayBlock{A,2} where A,
  b::ArrayBlock{B,1} where B,
  α::Number,β::Number)

  ni, nj = size(a)
  @check nj == size(b.array,1)
  for i in 1:ni
    if β!=1 && c.touched[i]
      rmul!(c.array[i],β)
    end
    for j in 1:nj
      if a.touched[i,j] && b.touched[j]
        ci = c.array[i]
        aij = a.array[i,j]
        bj = b.array[j]
        mul!(ci,aij,bj,α,1)
      end
    end
  end
  c
end

function LinearAlgebra.mul!(
  c::ArrayBlock{C,2} where C,
  a::ArrayBlock{A,2} where A,
  b::ArrayBlock{B,2} where B,
  α::Number,β::Number)

  ni, nk = size(a.array)
  nj = size(b.array,2)
  @check nk == size(b.array,1)
  @check (ni,nj) == size(c.array)

  for i in 1:ni
    for j in 1:nj
      if β!=1 && c.touched[i,j]
        rmul!(c.array[i,j],β)
      end
      for k in 1:nk
        if a.touched[i,k] && b.touched[k,j]
          cij = c.array[i,j]
          aik = a.array[i,k]
          bkj = b.array[k,j]
          _mymul!(cij,aik,bkj,α)
          # hack to avoid allocations
          #mul!(cij,aik,bkj,α,1)
        end
      end
    end
  end
  c
end

function  _setsize_mul!(c,a::AbstractMatrix,b::AbstractVector)
  ni = size(a,1)
  @check length(b) == size(a,2)
  setsize!(c,(ni,))
  c
end

function  _setsize_mul!(c,a::AbstractMatrix,b::AbstractMatrix)
  ni = size(a,1)
  nj = size(b,2)
  @check size(a,2) == size(b,1)
  setsize!(c,(ni,nj))
  c
end

function  _setsize_mul!(c,a::MatrixBlock,b::VectorBlock)
  ni,nj = size(a)
  @check length(b) == nj
  @check length(c) == ni
  for i in 1:ni
    for j in 1:nj
      if a.touched[i,j] && b.touched[j]
        _setsize_mul!(c.array[i],a.array[i,j],b.array[j])
      end
    end
  end
  c
end

function  _setsize_mul!(c,a::MatrixBlock,b::MatrixBlock)
  ni,nk = size(a)
  nk2,nj = size(b)
  @check nk == nk2
  @check size(c) == (ni,nj)
  for i in 1:ni
    for j in 1:nj
      for k in 1:nk
        if a.touched[i,k] && b.touched[k,j]
          _setsize_mul!(c.array[i,j],a.array[i,k],b.array[k,j])
        end
      end
    end
  end
  c
end

function return_cache(::typeof(*),a::ArrayBlock,b::ArrayBlock)
  c1 = CachedArray(a*b)
  c2 = return_cache(unwrap_cached_array,c1)
  (c1,c2)
end

function evaluate!(cache,::typeof(*),a::ArrayBlock,b::ArrayBlock)
  c1,c2 = cache
  _setsize_mul!(c1,a,b)
  c = evaluate!(c2,unwrap_cached_array,c1)
  mul!(c,a,b)
  c
end

function  _setsize_as!(d,a::AbstractArray)
  setsize!(d,size(a))
end

function  _setsize_as!(d,a::ArrayBlock)
  @check size(d) == size(a)
  for i in eachindex(a.array)
    if a.touched[i]
      _setsize_as!(d.array[i],a.array[i])
    end
  end
  d
end

function return_cache(k::MulAddMap,a::ArrayBlock,b::ArrayBlock,c::ArrayBlock)
  c1 = CachedArray(a*b+c)
  c2 = return_cache(unwrap_cached_array,c1)
  (c1,c2)
end

function evaluate!(cache,k::MulAddMap,a::ArrayBlock,b::ArrayBlock,c::ArrayBlock)
  c1,c2 = cache
  _setsize_as!(c1,c)
  _setsize_mul!(c1,a,b)
  d = evaluate!(c2,unwrap_cached_array,c1)
  copyto!(d,c)
  mul!(d,a,b,k.α,k.β)
  d
end

function Base.copyto!(d::ArrayBlock,c::ArrayBlock)
  _zero_entries!(d)
  for i in eachindex(c.array)
    if c.touched[i]
      copyto!(d.array[i],c.array[i])
    end
  end
  d
end

# Autodiff related

function return_cache(k::Arrays.ConfigMap{typeof(ForwardDiff.gradient)},x::VectorBlock)
  xi = testitem(x)
  fi = return_cache(k,xi)
  array = Vector{typeof(fi)}(undef,length(x.array))
  for i in eachindex(x.array)
    if x.touched[i]
      array[i] = return_cache(k,x.array[i])
    end
  end
  ArrayBlock(array,x.touched)
end

function return_cache(k::Arrays.ConfigMap{typeof(ForwardDiff.jacobian)},x::VectorBlock)
  xi = testitem(x)
  fi = return_cache(k,xi)
  array = Vector{typeof(fi)}(undef,length(x.array))
  for i in eachindex(x.array)
    if x.touched[i]
      array[i] = return_cache(k,x.array[i])
    end
  end
  ArrayBlock(array,x.touched)
end

function return_cache(k::Arrays.DualizeMap,x::VectorBlock)
  cfg = return_cache(Arrays.ConfigMap(k.f),x)
  i = first(findall(x.touched))
  xi = x.array[i]
  cfgi = cfg.array[i]
  xidual = evaluate!(cfgi,k,xi)
  array = Vector{typeof(xidual)}(undef,length(x.array))
  cfg, ArrayBlock(array,x.touched)
end

function evaluate!(cache,k::Arrays.DualizeMap,x::VectorBlock)
  cfg, xdual = cache
  for i in eachindex(x.array)
    if x.touched[i]
      xdual.array[i] = evaluate!(cfg.array[i],k,x.array[i])
    end
  end
  xdual
end

function return_cache(k::Arrays.AutoDiffMap,ydual::VectorBlock,x,cfg::VectorBlock)
  i = first(findall(ydual.touched))
  yidual = ydual.array[i]
  xi = x.array[i]
  cfgi = cfg.array[i]
  ci = return_cache(k,yidual,xi,cfgi)
  ri = evaluate!(ci,k,yidual,xi,cfgi)
  c = Vector{typeof(ci)}(undef,length(ydual.array))
  array = Vector{typeof(ri)}(undef,length(ydual.array))
  for i in eachindex(ydual.array)
    if ydual.touched[i]
      c[i] = return_cache(k,ydual.array[i],x.array[i],cfg.array[i])
    end
  end
  ArrayBlock(array,ydual.touched), c
end

function evaluate!(cache,k::Arrays.AutoDiffMap,ydual::VectorBlock,x,cfg::VectorBlock)
  r,c = cache
  for i in eachindex(ydual.array)
    if ydual.touched[i]
      r.array[i] = evaluate!(c[i],k,ydual.array[i],x.array[i],cfg.array[i])
    end
  end
  r
end

function return_cache(k::Arrays.AutoDiffMap,ydual::MatrixBlock,x,cfg::VectorBlock)
  ci = first(findall(ydual.touched))
  i, j = Tuple(ci)
  yidual = ydual.array[ci]
  xi = x.array[j]
  cfgi = cfg.array[j]
  ci = return_cache(k,yidual,xi,cfgi)
  ri = evaluate!(ci,k,yidual,xi,cfgi)
  c = Matrix{typeof(ci)}(undef,size(ydual))
  array = Matrix{typeof(ri)}(undef,size(ydual))
  ni,nj = size(ydual)
  for j in 1:nj
    for i in 1:ni
      if ydual.touched[i,j]
        c[i,j] = return_cache(k,ydual.array[i,j],x.array[j],cfg.array[j])
      end
    end
  end
  ArrayBlock(array,ydual.touched), c
end

function evaluate!(cache,k::Arrays.AutoDiffMap,ydual::MatrixBlock,x,cfg::VectorBlock)
  r,c = cache
  ni,nj = size(ydual)
  for j in 1:nj
    for i in 1:ni
      if ydual.touched[i,j]
        r.array[i,j] = evaluate!(c[i,j],k,ydual.array[i,j],x.array[j],cfg.array[j])
      end
    end
  end
  r
end

# hack to avoid allocations
function _mymul!(
  c::ArrayBlock{C,2} where C,
  a::ArrayBlock{A,2} where A,
  b::ArrayBlock{B,2} where B,
  α::Number)
  mul!(c,a,b,α,1)
end

function _mymul!(cIJ,aIK,bKJ,α)
  @boundscheck begin
    @assert size(cIJ,1) == size(aIK,1)
    @assert size(cIJ,2) == size(bKJ,2)
    @assert size(aIK,2) == size(bKJ,1)
  end
  for i in 1:size(cIJ,1)
    for j in 1:size(bKJ,2)
      for k in 1:size(bKJ,1)
        @inbounds cIJ[i,j] += α*aIK[i,k]*bKJ[k,j]
      end
    end
  end
end

# Assembly related

for T in (:AddEntriesMap,:TouchEntriesMap)
  @eval begin

    function return_cache(
      k::$T,A,v::MatrixBlock,I::VectorBlock,J::VectorBlock)

      qs = findall(v.touched)
      i, j = Tuple(first(qs))
      cij = return_cache(k,A,v.array[i,j],I.array[i],J.array[j])
      ni,nj = size(v.touched)
      cache = Matrix{typeof(cij)}(undef,ni,nj)
      for j in 1:nj
        for i in 1:ni
          if v.touched[i,j]
            cache[i,j] = return_cache(k,A,v.array[i,j],I.array[i],J.array[j])
          end
        end
      end
      cache
    end

    function evaluate!(
      cache, k::$T,A,v::MatrixBlock,I::VectorBlock,J::VectorBlock)
      ni,nj = size(v.touched)
      for j in 1:nj
        for i in 1:ni
          if v.touched[i,j]
            evaluate!(cache[i,j],k,A,v.array[i,j],I.array[i],J.array[j])
          end
        end
      end
    end

    function return_cache(
      k::$T,A,v::VectorBlock,I::VectorBlock)

      qs = findall(v.touched)
      i = first(qs)
      ci = return_cache(k,A,v.array[i],I.array[i])
      ni = length(v.touched)
      cache = Vector{typeof(ci)}(undef,ni)
      for i in 1:ni
        if v.touched[i]
          cache[i] = return_cache(k,A,v.array[i],I.array[i])
        end
      end
      cache
    end

    function evaluate!(
      cache, k::$T,A,v::VectorBlock,I::VectorBlock)
      ni = length(v.touched)
      for i in 1:ni
        if v.touched[i]
          evaluate!(cache[i],k,A,v.array[i],I.array[i])
        end
      end
    end

  end
end


text:
numbers = parse.(Int, split(readline("day4/bingo.txt"), ','))

boards = open("day4/bingo.txt") do io
    readline(io)  # remove first line
    readlines(io)
end

boards = filter(!isempty, boards)
boards = permutedims(hcat(split.(boards)...))
boards = parse.(Int, boards)
boards = [boards[i:i+4, :] for i in 1:5:size(boards, 1)]

function play_bingo(numbers, boards)
    bingos = []
    marked_boards = [zeros(Bool, size(board)) for board in boards]
    is_bingo_board = zeros(Bool, length(boards))
    for number in numbers
        for i in eachindex(boards, marked_boards, is_bingo_board)
            is_bingo_board[i] && continue
            board = boards[i]
            marked_board = marked_boards[i]
            for j in eachindex(board, marked_board)
                if board[j] == number
                    marked_board[j] = true
                end
            end
            is_marked_col = sum(marked_board, dims = 1) .== 5
            is_marked_row = sum(marked_board, dims = 2) .== 5
            is_bingo_board[i] = any(is_marked_col) | any(is_marked_row)
            is_bingo_board[i] && push!(bingos, (number, board, marked_board))
        end
    end
    bingos
end

bingos = play_bingo(numbers, boards)

number, board, marked_board = bingos[1]
p1 = number * sum(marked_board[i] ? 0 : board[i] for i in eachindex(board))

@show p1

number, board, marked_board = bingos[end]
p2 = number * sum(marked_board[i] ? 0 : board[i] for i in eachindex(board))

@show p2


text:
"""
    BoltzmannNetworks

Implements neural sampling in second-order Boltzmann distributions as presented in

> Lars Buesing, Johannes Bill, Bernhard Nessler and Wolfgang Maass (2011).
> Neural Dynamics as Sampling: A Model for Stochastic Computation in Recurrent Networks of Spiking Neurons.
> PLoS Computational Biology 7(11), e1002211.

The goal is to approximate a distribution ``p(\\mathbf{z})`` over the binary variables
``\\mathbf{z} = z_1, z_2, \\ldots, z_k`` by a neural network ``\\mathcal{N}`` of neurons
``v_1, v_2, \\ldots, v_k``. The spiking activity of ``\\mathcal{N}`` is interpreted as a
Markov chain implementing Markov Chain Monte Carlo. Specifically, a spike of neuron ``v_i``
sets the variable ``z_i`` to `true` for ``\\tau`` timesteps / milliseconds / samples. The
passing of this period is kept track of via auxiliary variables ``\\zeta_1, \\zeta_2, \\ldots, \\zeta_k``.

Networks are represented by the type `BoltzmannNetwork`. This module implements the
following methods from the paper:

- sampling with absolute refractory period in discrete time via the function `sample_absolute!`
- sampling with relative refractory period in discrete time via the function `sample_relative!`
- Gibbs sampling (not using neural dynamics) as a baseline method

"""
module BoltzmannNetworks

using Random: shuffle
using ..Utils: σ
import ..Utils: spikes2chain
import Base: show, length

export BoltzmannNetwork
export randomnetwork, sample_absolute!, sample_relative!, sample_clamp!

###############################################################################
# Network type                                                                #
###############################################################################
"""
    BoltzmannNetwork

Represents a network of spiking neurons that performs MCMC sampling from the distribution

```math
p(\\mathbf{z}) = \\frac{1}{Z} \\exp \\left( \\sum_{i,j} \\frac{1}{2} W_{ij} z_i z_j + \\sum_i b_i z_i \\right)
```

where `Z` is the partition function, `W` are weights and `b` are biases.

# Fields
- `ζ::Vector{Int}` - vector of auxiliary variables `ζ`
- `z::Vector{Bool}` - vector of variables `z`
- `u::Vector{Float64}` - vector of neural membrane potentials
- `W::Matrix{Float64}` - weights relating the variables (symmetric, zeros on main diagonal)
- `b::Vector{Float64}` - biases governing neural excitability
- `τ::Int` - time constant encoding for how many steps `z` is `true` after a spike
- `g::G` - refractoriness function for sampling with relative refr. period
- `f::F` - activation function for sampling with relative refr. period
"""
struct BoltzmannNetwork{G, F}
    ζ::Vector{Int}
    z::Vector{Bool}
    u::Vector{Float64}
    W::Matrix{Float64}
    b::Vector{Float64}
    τ::Int
    g::G
    f::F

    # constructor to enforce constraints on W
    function BoltzmannNetwork(ζ, z, u, W, b, τ, g, f)
        for i in 1:size(W, 1)
            @assert W[i, i] == 0.0 "Main diagonal of weight matrix has to be all zeros."
            for j in i:size(W, 2)
                @assert W[i, j] == W[j, i] "Weight matrix has to be symmetric."
            end
        end
        return new{typeof(g), typeof(f)}(ζ, z, u, W, b, τ, g, f)
    end
end

"Generate a random `BoltzmannNetwork` of `n` neurons."
function randomnetwork(n::Int; τ=20, g=g) # τ as in Buesing et al. (2011)
    ζ = zeros(Int, n)
    z = zeros(Bool, n)
    u = zeros(Float64, n)
    W = randn(n, n) * 0.3 # scaling as in Buesing et al. (2011), Fig. 3, 5, 6, 7
    for i in 1:n, j in i:n
        if i==j
            W[i, j] = 0. # enforce that main axis is all zeros
        else
            W[i, j] = W[j, i] # enforce that weights are symmetric
        end
    end
    b = randn(n) * 0.5 .- 1.5 # scaling as in Buesing et al. (2011), Fig. 3, 5, 6, 7

    f = gettransferfun(g, τ=τ)
    return BoltzmannNetwork(ζ, z, u, W, b, τ, g, f)
end

###############################################################################
# Sampling with absolute refractory period.                                   #
###############################################################################
"""
    step_absolute!(net::BoltzmannNetwork)

Implements a single step of neural MCMC sampling with absolute refractory period
in discrete time.

# Details

Alters the state of the membrane potentials `net.u`, the variables `net.z` and
the auxiliary variables `net.ζ`. The neurons are updated sequentially (as
explained in the paper, this is necessary for correct sampling). For each
neuron, the membrane potential is calculated in accordance with the
**neural computability condition**:

```math
u_i = b_i + sum_j W_ij * z_i * z_j
```

If the neuron is refractory (``\\zeta \\ge 1``), ζ is decreased deterministically.
Otherwise the neurons spikes with probability ``\\sigma(u - \\log(\\tau))``. A spike
sets `z` to `true` and `ζ` to `τ`.
"""
function step_absolute!(net::BoltzmannNetwork)
    spikes = zeros(Bool, length(net.z))
    # update neurons
    # (in order, since transition operators have to be applied sequentially)
    for i in 1:length(spikes)
        # update membrane potential
        net.u[i] = net.b[i] + net.W[i, :]' * net.z
        # generate spike
        if (net.ζ[i] ≤ 1) && (rand() < σ(net.u[i] - log(net.τ)))
            spikes[i] = true
            net.ζ[i] = net.τ
            net.z[i] = 1
        else
            net.ζ[i] = max(net.ζ[i] - 1)
            if net.ζ[i] == 0
                net.z[i] = 0
            end
        end
    end
    return(spikes)
end

"""
    sample_absolute!(net, steps, burnin, dt)

Performs neural sampling with absolute refractory period in discrete time in
`BoltzmannNet` `net` for `steps` steps after a burnin time of `burnin` steps
with step size `dt`.

# Returns
- `ts::Vector{Float64}` - vector of spike times
- `inds::Vector{Int}` - vector of neuron indices of neurons that spiked
"""
function sample_absolute!(net, steps, burnin, dt)
    for _ in 1:burnin
        step_absolute!(net)
    end
    ts = zeros(0)
    inds = zeros(Int, 0)
    for t in 1:steps
        spikes = step_absolute!(net)
        for spike in findall(spikes)
            push!(ts, t * dt)
            push!(inds, spike)
        end
    end
    return ts, inds
end

###############################################################################
# Sampling with relative refractory period                                    #
###############################################################################
#
# In sampling with a relative refractory period, the neural transfer function
# is slightly more complicated. The refactoriness of a neuron is given by the
# function g(ζ): [0, τ] -> [0, 1]
#
# Depending on the exact function g chosen, the neural transfer function f
# has to fulfill some specific conditions. In general, a suitable f exists but
# cannot be found analytically. Therefore, it has to be approximated numerically.
#
# First, we implement the function g given in the paper.
"""Refractory function g(ζ).
Quantifies how excitable a neuron is based on variable ζ,
which encodes the time since last spike.

Assumes that ζ is in range [0, 1]"""
function g(ζ)
    x = 1 - ζ + sin(2π * ζ) / (2π) # according to Table 1 in Buesing et al. (2011)
    return min(1, max(0, x))
end

"Alternative refractory function"
function g2(ζ) # assumes that ζ is in range [0, 1]
    x = 4 * (1 - ζ) + sin(8π * ζ) / (2π) # according to Table 1 in Buesing et al. (2011)
    return min(1, max(0, x))
end

"Alternative refractory function"
function g3(ζ) # assumes that ζ is in range [0, 1]
    x = 1 - 2ζ + sin(4π * ζ) / (2π) # according to Table 1 in Buesing et al. (2011)
    return min(1, max(0, x))
end

# Next, we need to approximate f for a given g.
"""
    gettransferfun(g; τ=20)

Approximates the neural transfer function `f(u)` for a given function `g(ζ)`.

# Arguments
- `g::Function` - refractoriness; should accept values in range [0, 1].
"""
function gettransferfun(g; τ=20)
    # the following function calculates the error for a given u and f(u)
    ferror(u, fu) = abs(exp(u) - fu *
        sum(prod(1 - g(ζ/τ) * fu for ζ in η+1:τ) for η in 1:τ-1) /
        prod(1 - g(ζ/τ) * fu for ζ in 1:τ))
    # we assume that membrane potentials will stay in the range [-10, 10]
    us = -10:0.01:10
    fs = zeros(length(us))
    # we know that the output of f must lie in [0, 1]
    vals = 0:0.001:1
    # now, for each possible input value u, find the output value that minimizes the error
    for (i, u) in enumerate(us)
        fs[i] = vals[argmin(ferror.(u, vals))]
    end
    # now we can generate the function f:
    # for each input value, find the nearest known input value from us and
    # return the corresponding f(u)
    function f(u; us=us, fs=fs)
        i = argmin(abs.(u .- us))
        return fs[i]
    end
    return f
end

# Given the correct transfer function, we can implement the sampling
"""
    step_relative!(net::BoltzmannNetwork)

Implements a single step of neural MCMC sampling with relative refractory period
in discrete time.

# Details

Alters the state of the membrane potentials `net.u`, the variables `net.z` and
the auxiliary variables `net.ζ`. The neurons are updated sequentially (as
explained in the paper, this is necessary for correct sampling). For each
neuron, the membrane potential is calculated in accordance with the
**neural computability condition**:

```math
u_i = b_i + sum_j W_ij * z_i * z_j
```

The neural refactoriness is given by `g(ζ)`. At each point, the neuron fires with
probability `g(ζ) * f(u)`. A spike sets `z` to `true` and `ζ` to `τ`. If the neuron
does not spike, `ζ` decreases and `z` is set to `false` once `ζ` becomes `0`.
"""
function step_relative!(net::BoltzmannNetwork)
    spikes = zeros(Bool, length(net.z))
    # update neurons
    # (in order, since transition operators have to be applied sequentially)
    for i in 1:length(spikes)
        # update membrane potential
        net.u[i] = net.b[i] + net.W[i, :]' * net.z
        # generate spike with probability g(ζ) * f(u), otherwise decay ζ
        if (rand() < net.g(net.ζ[i] / net.τ) * net.f(net.u[i]))
            spikes[i] = true
            net.ζ[i] = net.τ
            net.z[i] = 1
        else
            net.ζ[i] = max(net.ζ[i] - 1)
            if net.ζ[i] == 0
                net.z[i] = 0
            end
        end
    end
    return(spikes)
end

"""
    sample_relative!(net, steps, burnin, dt)

Performs neural sampling with relative refractory period in discrete time in
`BoltzmannNet` `net` for `steps` steps after a burnin time of `burnin` steps
with step size `dt`.

# Returns
- `ts::Vector{Float64}` - vector of spike times
- `inds::Vector{Int}` - vector of neuron indices of neurons that spiked
"""
function sample_relative!(net, steps, burnin, dt)
    # run the Markov chain for burn-in period
    for _ in 1:burnin
        step_relative!(net)
    end

    # sample spikes
    ts = zeros(0)
    inds = zeros(Int, 0)
    for t in 1:steps
        spikes = step_relative!(net)
        for spike in findall(spikes)
            push!(ts, t * dt)
            push!(inds, spike)
        end
    end
    return ts, inds
end

###############################################################################
# Sampling with clamped values (inferring posterior)                          #
###############################################################################
"""
    step_clamp!(net::BoltzmannNetwork, clamp_on, clamp_off)

Performs one step of neural sampling with relative refractory period in `net`
while clamping the neurons at indices `clamp_on` to `true` and the neurons at
indices `clamp_off` to `false`.
"""
function step_clamp!(net::BoltzmannNetwork, clamp_on, clamp_off)
    spikes = zeros(Bool, length(net.z))
    # update neurons
    # (in order, since transition operators have to be applied sequentially)
    for i in shuffle(1:length(spikes))
        if i in clamp_on
            net.z[i] = true
            continue
        end
        if i in clamp_off
            net.z[i] = false
            continue
        end
        # update membrane potential
        net.u[i] = net.b[i] + net.W[i, :]' * net.z
        # generate spike with probability g(ζ) * f(u), otherwise decay ζ
        if (rand() < net.g(net.ζ[i] / net.τ) * net.f(net.u[i]))
            spikes[i] = true
            net.ζ[i] = net.τ
            net.z[i] = 1
        else
            net.ζ[i] = max(net.ζ[i] - 1)
            if net.ζ[i] == 0
                net.z[i] = 0
            end
        end
    end
    return(spikes)
end

"""
    sample_clamp!(net, steps, burnin, dt, clamp_on, clamp_off)

Performs neural sampling with relative refractory period in discrete time in
`BoltzmannNet` `net` for `steps` steps after a burnin time of `burnin` steps
with step size `dt` while clamping the neurons at indices `clamp_on` to `true`
and the neurons at indices `clamp_off` to `false`.

# Returns
- `ts::Vector{Float64}` - vector of spike times
- `inds::Vector{Int}` - vector of neuron indices of neurons that spiked
"""
function sample_clamp!(net, steps, burnin, dt, clamp_on, clamp_off)
    # run the Markov chain for burn-in period
    for _ in 1:burnin
        step_clamp!(net, clamp_on, clamp_off)
    end

    # sample spikes
    ts = zeros(0)
    inds = zeros(Int, 0)
    for t in 1:steps
        spikes = step_clamp!(net, clamp_on, clamp_off)
        for spike in findall(spikes)
            push!(ts, t * dt)
            push!(inds, spike)
        end
    end
    return ts, inds
end

###############################################################################
# Conversion of spike sequence into Markov chain                              #
###############################################################################
function spikes2chain(ts, inds, net, dt; start=0, stop=nothing)
    if stop == nothing
        stop = maximum(ts) + net.τ * dt
    end
    timesteps = round(Int, (stop - start) / dt)
    chain = zeros(Bool, length(net), timesteps)
    for (t, i) in zip(ts, inds)
        # each spike sets the corresponding variable to true for τ timesteps
        tᵢₙₜ = max(1, round(Int, (t - start) / dt))
        chain[i, tᵢₙₜ:min(timesteps, tᵢₙₜ+net.τ)] .= true
    end
    return chain
end

###############################################################################
# Overloads for pretty printing etc.                                          #
###############################################################################
function Base.show(io::IO, net::BoltzmannNetwork)
    print(io, "BoltzmannNetwork of $(length(net.z)) neurons.")
end

function Base.show(io::IO, ::MIME"text/plain", net::BoltzmannNetwork)
    print(io,
        "BoltzmannNetwork of $(length(net.z)) neurons.\n",
        "\t z = $(net.z) \n",
        "\t u = $(net.u) \n",
        "\t W = $(net.W) \n",
        "\t b = $(net.b) \n",
        "\t τ = $(net.τ)"
    )
end

function Base.length(net::BoltzmannNetwork)
    return length(net.z)
end

end # module BoltzmannNetworks