text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
Lara = Lar.Arrangement

#--------------------------------------------------------------------------------------------------------------------------------
#                               PIPELINE PART 1
#--------------------------------------------------------------------------------------------------------------------------------

"""
    frag_edge_channel(in_chan, out_chan, V, EV, bigPI)

Utility function for parallel edge fragmentation.

This function handles the edge fragmentation in the first part of arrangement's algorithmic pipeline
(see also [`Lar.Arrangement.planar_arrangement_1`](@ref)) during multiprocessing computation.
In order to do so it needs two `Distributed.RemoteChannel`s, one with the inputs and one for outputs.

See also: [`Lar.Arrangement.frag_edge`](@ref)
"""
function frag_edge_channel( in_chan::Distributed.RemoteChannel{Channel{Int64}},
                            out_chan::Distributed.RemoteChannel{Channel{Tuple}},
                            V::Lar.Points,
                            EV::Lar.ChainOp,
                            bigPI::Array{Array{Int64,1},1})
    run_loop = true
    while run_loop
        edgenum = take!(in_chan)
        if edgenum != -1
            put!(out_chan, (edgenum, Lar.Arrangement.frag_edge(V, EV, edgenum, bigPI)))
        else
            run_loop = false
        end
    end
end

"""
	frag_edge(V, EV, edge_idx, bigPI)

Splits the `edge_idx`-th edge in `EV`.

This method splits the `edge_idx`-th edge in `EV` into several parts by confronting it
with the others that intersect its bounding box (see also [`Lar.Arrangement.intersect_edges`](@ref)).

The method returns a set of the new vertices that the segment is made of (with redundancies) and
the associated cochain (with no redundancies).

See also:
 - [`Lar.Arrangement.planar_arrangement_1`](@ref)
 - [`Lar.Arrangement.frag_edge_channel`](@ref)

---

# Examples
```jldoctest
julia> V = [1.0 0.0; 0.0 1.0; 0.0 0.5; 0.5 1.0; 1.0 1.0]; # By Rows!

julia> EV = [[1, 2], [2, 5], [3, 4], [4, 5]];

julia> copEV = Lar.coboundary_0(EV::Lar.Cells);

julia> bigPI = Lar.spaceindex((convert(Lar.Points, V'), EV));

julia> Lara.frag_edge(V, copEV, 1, bigPI)[1]
5×2 Array{Float64,2}:
 1.0   0.0
 0.0   1.0
 1.0   0.0
 0.25  0.75
 0.0   1.0

julia> Lara.frag_edge(V, copEV, 1, bigPI)[2]
2×5 SparseMatrixCSC{Int8,Int64} with 4 stored entries:
  [1, 1]  =  1
  [2, 2]  =  1
  [1, 4]  =  1
  [2, 4]  =  1
```
"""
function frag_edge(V::Lar.Points, EV::Lar.ChainOp, edge_idx::Int, bigPI::Array{Array{Int64,1},1})::Tuple{Lar.Points, Lar.ChainOp}
    alphas = Dict{Float64, Int}()
    edge = EV[edge_idx, :]
    verts = V[edge.nzind, :]
    for i in bigPI[edge_idx]
        if i != edge_idx && edge_idx in bigPI[i]    #<------------------------- Could be usefull?
            intersection = Lar.Arrangement.intersect_edges(V, edge, EV[i, :])
            for (point, alpha) in intersection
                verts = [verts; point]
                alphas[alpha] = size(verts, 1)
            end
        end
    end
    alphas[0.0], alphas[1.0] = [1, 2]
    alphas_keys = sort(collect(keys(alphas)))
    edge_num = length(alphas_keys)-1
    verts_num = size(verts, 1)
    ev = SparseArrays.spzeros(Int8, edge_num, verts_num)
    for i in 1:edge_num
        ev[i, alphas[alphas_keys[i]]] = 1
        ev[i, alphas[alphas_keys[i+1]]] = 1
    end
    return verts, ev
end

"""
	intersect_edges(V, edge1, edge2)

Finds the intersection point (if there exists) between the two given edges.

This method compute the points where `edge2` intersect `edge1`.
If they are collinear only the vertices of `edge2` are considered (see the second example).

See also: [`Lar.Arrangement.frag_edge`](@ref)

---

# Examples
```jldoctest
# Cross
julia> V = [1.0 0.0; 0.0 1.0; 0.0 0.5; 0.5 1.0]; # By Rows!

julia> EV = [[1, 2], [3, 4]];

julia> copEV = Lar.coboundary_0(EV::Lar.Cells);

julia> Lara.intersect_edges(V, copEV[1, :], copEV[2, :])
1-element Array{Tuple{Array{T,2} where T,Float64},1}:
 ([0.25 0.75], 0.75)
```

```jldoctest
# Collinear
julia> V = [1.0 0.0; 0.0 1.0; 0.75 0.25; 0.5 0.5]; # By Rows!

julia> EV = [[1, 2], [3, 4]];

julia> copEV = Lar.coboundary_0(EV::Lar.Cells);

julia> Lar.Arrangement.intersect_edges(V, copEV[1, :], copEV[2, :])
2-element Array{Tuple{Array{T,2} where T,Float64},1}:
 ([0.75 0.25], 0.25)
 ([0.5 0.5], 0.5)

julia> Lar.Arrangement.intersect_edges(V, copEV[2, :], copEV[1, :])
0-element Array{Tuple{Array{T,2} where T,Float64},1}
```
"""
function intersect_edges(V::Lar.Points, edge1::Lar.Cell, edge2::Lar.Cell)::Array{Tuple{Lar.Points, Float64}, 1}
    err = 10e-8

    x1, y1, x2, y2 = vcat(map(c->V[c, :], edge1.nzind)...)
    x3, y3, x4, y4 = vcat(map(c->V[c, :], edge2.nzind)...)
    ret = Array{Tuple{Lar.Points, Float64}, 1}()

    v1 = [x2-x1, y2-y1];
    v2 = [x4-x3, y4-y3];
    v3 = [x3-x1, y3-y1];

    ang1 = Lar.dot(Lar.normalize(v1), Lar.normalize(v2))
    ang2 = Lar.dot(Lar.normalize(v1), Lar.normalize(v3))

    parallel = 1-err < abs(ang1) < 1+err
    colinear = parallel && (1-err < abs(ang2) < 1+err || -err < Lar.norm(v3) < err)

    if colinear
        o = [x1 y1]
        v = [x2 y2] - o
        alpha = 1/Lar.dot(v,v')
        ps = [x3 y3; x4 y4]
        for i in 1:2
            a = alpha*Lar.dot(v',(reshape(ps[i, :], 1, 2)-o))
            if 0 < a < 1
                push!(ret, (ps[i:i, :], a))
            end
        end

    elseif !parallel
        denom = (v2[2])*(v1[1]) - (v2[1])*(v1[2])
        a = ((v2[1])*(-v3[2]) - (v2[2])*(-v3[1])) / denom
        b = ((v1[1])*(-v3[2]) - (v1[2])*(-v3[1])) / denom

        if -err < a < 1+err && -err <= b <= 1+err
            p = [(x1 + a*(x2-x1))  (y1 + a*(y2-y1))]
            push!(ret, (p, a))
        end
    end

    return ret
end

"""
    merge_vertices!(V, EV[, edge_map[, err]])

Compact the vertices closer than `err` in a single one.

This method check one at time each vertex ``v`` in `V` and identifies each other vertex within `err` with ``v`` itself.
The cochain `EV` is coherently modified (multiple edges between two vertices are not allowed).
If an `edge_map` is given in input (this could be usefull during the planar arrangements), then also
the map is coherently modified and given back in output.

See also: [`Lar.Arrangement.planar_arrangement_1`](@ref)

---

# Examples
```jldoctest
# Collinear
julia> V = [0.5 0.5; 0.0 0.0; 0.5 0.5; 1.0 1.0; 0.5 0.5; 1.0 1.0]; # By Rows!

julia> EV = [[1, 4], [3, 2], [5, 6], [1, 6], [5, 3]];

julia> copEV = Lar.coboundary_0(EV::Lar.Cells);

julia> Lara.merge_vertices!(V, copEV)[1]
3×2 Array{Float64,2}:
 0.5  0.5
 0.0  0.0
 1.0  1.0

julia> Lara.merge_vertices!(V, copEV)[2]
2×3 SparseArrays.SparseMatrixCSC{Int8,Int64} with 4 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
  [2, 2]  =  1
  [1, 3]  =  1
```
"""
function merge_vertices!(V::Lar.Points, EV::Lar.ChainOp, edge_map::Array{Array{Int64,1},1}=[[-1]], err=1e-4)
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    kdtree = KDTree(permutedims(V))

    todelete = []

    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = Lar.inrange(kdtree, V[vi, :], err)

            newverts[nearvs] .= i

            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)

            i = i + 1
        end
    end

    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)

    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind

        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))

    end
    nedges = union(edges)
    nedges = filter(t->t[1]!=t[2], nedges)

    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))

    etuple2idx = Dict{Tuple{Int, Int}, Int}()

    for ei in 1:nedgenum
        nEV[ei, collect(nedges[ei])] .= 1
        etuple2idx[nedges[ei]] = ei
    end

    if edge_map != [[-1]]
        for i in 1:length(edge_map)
            row = edge_map[i]
            row = map(x->edges[x], row)
            row = filter(t->t[1]!=t[2], row)
            row = map(x->etuple2idx[x], row)
            edge_map[i] = row
        end
        return Lar.Points(nV), nEV, edge_map
    end

    return Lar.Points(nV), nEV
end

#--------------------------------------------------------------------------------------------------------------------------------
#                                       PIPELINE BICONNECTED COMPONENTS
#--------------------------------------------------------------------------------------------------------------------------------

"""
    biconnected_components(EV)

Compute sets of 2-cells on the same 3-cells biconnected components.

The method evaluate the 2-cells ``σ_i`` wich lies on the same 3-cells biconnected
component ``χ_j`` and gives back an array that contains an array of ``σ_i`` for each ``j``.
Do note that if the cochain EV contains more copies of the same 2-cell then
it will be considered like a 3-cell.

Intermediate part of Planar Arrangement's algorithmic pipeline.

See also: [`Lar.planar_arrangement`](@ref) for the complete pipeline.

---

# Examples
```jldoctest
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
    [1 1 0 0 0 0] #1 -> 1,2  |
    [1 0 1 0 0 0] #2 -> 1,3  |
    [1 0 0 1 0 0] #3 -> 1,4   |
    [1 0 0 0 1 0] #4 -> 1,5   |
    [1 0 0 0 0 1] #5 -> 1,6
    [0 1 1 0 0 0] #6 -> 2,3  |
    [0 0 0 1 1 0] #7 -> 4,5   |
    ]));
julia> Lara.biconnected_components(copEV)
2-element Array{Array{Int64,1},1}:
 [2, 6, 1]
 [4, 7, 3]
```

```jldoctest
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
    [1 1 0] #1 -> 1,2  |
    [1 1 0] #2 -> 1,2  |
    [1 0 1] #3 -> 1,2
    ]));
julia> Lara.biconnected_components(copEV)
1-element Array{Array{Int64,1},1}:
 [2, 1]
```
"""
function biconnected_components(EV::Lar.ChainOp)::Array{Array{Int, 1}, 1}
    ps = Array{Tuple{Int, Int, Int}, 1}()
    es = Array{Tuple{Int, Int}, 1}()
    todel = Array{Int, 1}()
    visited = Array{Int, 1}()
    bicon_comps = Array{Array{Int, 1}, 1}()
    hivtx = 1

    function an_edge(point)
        edges = setdiff(EV[:, point].nzind, todel)
        if length(edges) == 0
            edges = [false]
        end
        edges[1]
    end

    function get_head(edge, tail)
        setdiff(EV[edge, :].nzind, [tail])[1]
    end

    function v_to_vi(v)
        i = findfirst(t->t[1]==v, ps)
        # seems findfirst changed from 0 to Nothing
        if typeof(i) == Nothing
            return false
        elseif i == 0
            return false
        else
            return ps[i][2]
        end
    end

    push!(ps, (1,1,1))
    push!(visited, 1)
    exit = false
    while !exit
        edge = an_edge(ps[end][1])
        if edge != false
            tail = ps[end][2]
            head = get_head(edge, ps[end][1])
            hi = v_to_vi(head)
            if hi == false
                hivtx += 1
                push!(ps, (head, hivtx, ps[end][2]))
                push!(visited, head)
            else
                if hi < ps[end][3]
                    ps[end] = (ps[end][1], ps[end][2], hi)
                end
            end
            push!(es, (edge, tail))
            push!(todel, edge)
        else
            if length(ps) == 1
                found = false
                pop!(ps)
                for i in 1:size(EV,2)
                    if !(i in visited)
                        hivtx = 1
                        push!(ps, (i, hivtx, 1))
                        push!(visited, i)
                        found = true
                        break
                    end
                end
                if !found
                    exit = true
                end

            else
                if ps[end][3] == ps[end-1][2]
                    edges = Array{Int, 1}()
                    while true
                        edge, tail = pop!(es)
                        push!(edges, edge)
                        if tail == ps[end][3]
                            if length(edges) > 1
                                push!(bicon_comps, edges)
                            end
                            break
                        end
                    end

                else
                    if ps[end-1][3] > ps[end][3]
                        ps[end-1] = (ps[end-1][1], ps[end-1][2], ps[end][3])
                    end
                end
                pop!(ps)
            end
        end
    end
    bicon_comps = sort(bicon_comps, lt=(x,y)->length(x)>length(y))
    return bicon_comps
end

#--------------------------------------------------------------------------------------------------------------------------------
#                                       PIPELINE PART 2
#--------------------------------------------------------------------------------------------------------------------------------

"""
    get_external_cycle(V, EV, FE)

Evaluates the index of the external 3D-cell.

This method looks for and retrieve the external cycle of the model `(V, EV, FE)`.
If the cell does not exist then it return a `Nothing` element.

See also: [`Lar.Arrangement.componentgraph`](@ref).

---

# Examples
```jldoctest
# ``K^4`` graph
julia> V = [0.0 0.0; 4.0 0.0; 2.0 3.0; 2.0 1.5];
julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 0 0] #1 -> 1,2
		[0 1 1 0] #2 -> 2,3
		[1 0 1 0] #3 -> 3,1
		[1 0 0 1] #4 -> 1,2
		[0 1 0 1] #5 -> 2,3
		[0 0 1 1] #6 -> 3,1
	]));
julia> copFE = SparseArrays.sparse(Array{Int8, 2}([
		[1 0 0 1 1 0] #1 -> 1,4,5
		[0 1 0 0 1 1] #2 -> 2,5,6
		[0 0 1 1 0 1] #3 -> 3,4,6
		[1 1 1 0 0 0] #4 -> 1,2,3 External
	]));

julia> Lara.get_external_cycle(V, copEV, copFE)
4

julia> typeof(Lara.get_external_cycle(V, copEV, copFE[1:3,:]))
Nothing
```
"""
function get_external_cycle(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp)
    FV = abs.(FE)*EV
    vs = sparsevec(mapslices(sum, abs.(EV), dims=1)').nzind
    minv_x1 = maxv_x1 = minv_x2 = maxv_x2 = pop!(vs)
    for i in vs
        if V[i, 1] > V[maxv_x1, 1]
            maxv_x1 = i
        elseif V[i, 1] < V[minv_x1, 1]
            minv_x1 = i
        end
        if V[i, 2] > V[maxv_x2, 2]
            maxv_x2 = i
        elseif V[i, 2] < V[minv_x2, 2]
            minv_x2 = i
        end
    end
    cells = intersect(
        FV[:, minv_x1].nzind,
        FV[:, maxv_x1].nzind,
        FV[:, minv_x2].nzind,
        FV[:, maxv_x2].nzind
    )
    if length(cells) == 1
        return cells[1]
    else
        for c in cells
            if Lar.face_area(V, EV, FE[c, :]) < 0
                return c
            end
        end
    end
end

"""
    pre_containment_test(bboxes)

Generate the containment graph associated to `bboxes`.

The function evaluate a `SparseArrays{Int8, n, n}` PRE-containmnet graph in the sense
that it stores `1` when the bbox associated to the `i`-th component contains entirely
the bbox associated to the `j`-th one. (where `n = length(bboxes)`)

See also:
  - [`Lar.Arrangement.componentgraph`](@ref).
  - [`Lar.bbox_contains`](@ref).

---

# Examples
```jldoctest
# Planar Tiles
julia> bboxes = [
		([0.0 0.0], [1.0 1.0])
		([0.0 0.0], [0.4 0.4])
		([0.6 0.0], [1.0 0.4])
		([0.0 0.6], [0.4 1.0])
		([0.6 0.6], [1.0 1.0])
		([2.0 3.0], [2.0 3.0]) # no intersection
	];

julia> Lara.pre_containment_test(bboxes)
5×5 SparseMatrixCSC{Int8,Int64} with 4 stored entries:
  [2, 1]  =  1
  [3, 1]  =  1
  [4, 1]  =  1
  [5, 1]  =  1
```
"""
function pre_containment_test(bboxes)
    n = length(bboxes)
    containment_graph = SparseArrays.spzeros(Int8, n, n)

    for i in 1:n
        for j in 1:n
            if i != j && Lar.bbox_contains(bboxes[j], bboxes[i])
                containment_graph[i, j] = 1
            end
        end
    end

    return containment_graph
end

"""
    prune_containment_graph(n, V, EVs, shells, graph)

Prunes the containment `graph` from the non-included faces.

This method prunes the containment graph eliminating the included bouning boxes
that not corresponds to included faces.

Do note that this method expects to work on biconnectet clusters of faces.

See also:
  - [`Lar.Arrangement.componentgraph`](@ref).
  - [`Lar.point_in_face`](@ref).
```
"""
function prune_containment_graph(n, V, EVs, shells, graph)

    for i in 1:n
        an_edge = shells[i].nzind[1]
        origin_index = EVs[i][an_edge, :].nzind[1]
        origin = V[origin_index, :]

        for j in 1:n
            if i != j
                if graph[i, j] == 1
                    shell_edge_indexes = shells[j].nzind
                    ev = EVs[j][shell_edge_indexes, :]

                    if !Lar.point_in_face(origin, V, ev)
                        graph[i, j] = 0
                    end
                end
             end
         end

     end
     return graph
end

"""
    transitive_reduction!(graph)

Prune the redundancies in containment `graph`.

The ``m``-layer containment `graph` is reduced to the 1-layer containment `graph`
where only the first level nested components are considered.

See also: [`Lar.Arrangement.componentgraph`](@ref).

---

# Examples
```jldoctest
julia> containment = [
		0 0 0 0 0 0
		0 0 0 0 0 0
		1 0 0 0 0 0
		1 0 0 0 0 0
		1 0 1 0 0 0
		1 0 1 0 1 0
	];

julia> transitive_reduction!(containment)

julia> containment
6×6 Array{Int64,2}:
 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  1  0
```
"""
function transitive_reduction!(graph)
    n = size(graph, 1)
    for j in 1:n
        for i in 1:n
            if graph[i, j] > 0
                for k in 1:n
                    if graph[j, k] > 0
                        graph[i, k] = 0
                    end
                end
            end
        end
    end
end

"""
    cell_merging(n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)

Cells composing for the Topological Gift Wrapping algorithm.

This is the online part of the TGW algorithm.

See also: [`Lar.Arrangement.planar_arrangement_2`](@ref).
"""
function cell_merging(n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)
    function bboxes(V::Lar.Points, indexes::Lar.ChainOp)
        boxes = Array{Tuple{Any, Any}}(undef, indexes.n)
        for i in 1:indexes.n
            v_inds = indexes[:, i].nzind
            boxes[i] = Lar.bbox(V[v_inds, :])
        end
        boxes
    end
    # initialization
    sums = Array{Tuple{Int, Int, Int}}(undef, 0);
	# assembling child components with father components
    for father in 1:n
        if sum(containment_graph[:, father]) > 0
            father_bboxes = bboxes(V, abs.(EVs[father]')*abs.(boundaries[father]'))
            for child in 1:n
                if containment_graph[child, father] > 0
                    child_bbox = shell_bboxes[child]
                    for b in 1:length(father_bboxes)
                        if Lar.bbox_contains(father_bboxes[b], child_bbox)
                            push!(sums, (father, b, child))
                            break
                        end
                    end
                end
            end
        end
    end
    # offset assembly initialization
    EV = vcat(EVs...)
    edgenum = size(EV, 1)
    facenum = sum(map(x->size(x,1), boundaries))
    FE = spzeros(Int8, facenum, edgenum)
    shells2 = spzeros(Int8, length(shells), edgenum)
    r_offsets = [1]
    c_offset = 1
    # submatrices construction
    for i in 1:n
        min_row = r_offsets[end]
        max_row = r_offsets[end] + size(boundaries[i], 1) - 1
        min_col = c_offset
        max_col = c_offset + size(boundaries[i], 2) - 1
        FE[min_row:max_row, min_col:max_col] = boundaries[i]
        shells2[i, min_col:max_col] = shells[i]
        push!(r_offsets, max_row + 1)
        c_offset = max_col + 1
    end
    # offsetting assembly of component submatrices
    for (f, r, c) in sums
        FE[r_offsets[f]+r-1, :] += shells2[c, :]
    end

    return EV, FE
end

"""
    componentgraph(V, copEV, bicon_comps)

Topological Gift Wrapping algorithm on 2D skeletons.

This is the offline part of the TGW algorithm. It takes in input a model and its
biconnected components mapping and evaluates usefull informations:
 1. Number of biconnected components.
 2. Component Graph of the biconnected structure.
 3. The 1-cells structure (UNMODIFIED). <----------------------------- Could be removed?
 4. Association between non-dangling 2-cells and their orientation (for each component).
 5. Association between 3-cells and 2-cells (with orientation, for each component).
 6. Association between 3-cells and their orientation (for each component).
 7. Shell bounding boxes of the components.

See also: [`Lar.Arrangement.planar_arrangement_2`](@ref) for the TGW.

---

# Examples
```jldoctest
# Papillon
julia> V = [0.0 0.0; 0.0 3.0; 2.0 1.5; 4.0 0.0; 4.0 3.0];

julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 0 0 0] #1 -> 1,2
		[0 1 1 0 0] #2 -> 2,3
		[1 0 1 0 0] #3 -> 3,1
		[0 0 1 1 0] #4 -> 3,4
		[0 0 0 1 1] #5 -> 4,5
		[0 0 1 0 1] #6 -> 3,5
	]));

julia> copFE = SparseArrays.sparse(Array{Int8, 2}([
		[1 1 1 0 0 0] #1 -> 1,2,3
		[0 0 0 1 1 1] #2 -> 4,5,6
		[1 1 1 1 1 1] #3 -> 1,2,3,4,5,6 External
	]));

julia> bicon_comps = [[1, 2, 3], [4, 5, 6]];

julia> Lara.componentgraph(V, copEV, bicon_comps)[1]
2

julia> Lara.componentgraph(V, copEV, bicon_comps)[2]
2×2 SparseMatrixCSC{Int8,Int64} with 0 stored entries

julia> Lara.componentgraph(V, copEV, bicon_comps)[4]
2-element Array{SparseMatrixCSC{Int8,Int64},1}:

  [1, 1]  =  -1
  [3, 1]  =  -1
  [1, 2]  =  1
  [2, 2]  =  -1
  [2, 3]  =  1
  [3, 3]  =  1

  [1, 3]  =  -1
  [3, 3]  =  -1
  [1, 4]  =  1
  [2, 4]  =  -1
  [2, 5]  =  1
  [3, 5]  =  1

julia> Lara.componentgraph(V, copEV, bicon_comps)[5]
2-element Array{SparseMatrixCSC{Int8,Int64},1}:

  [1, 1]  =  -1
  [1, 2]  =  -1
  [1, 3]  =  1

  [1, 1]  =  1
  [1, 2]  =  1
  [1, 3]  =  -1

julia> Lara.componentgraph(V, copEV, bicon_comps)[6]
2-element Array{SparseVector{Int8,Int64},1}:
   [1]  =  1
  [2]  =  1
  [3]  =  -1
   [1]  =  -1
  [2]  =  -1
  [3]  =  1

julia> Lara.componentgraph(V, copEV, bicon_comps)[7]
2-element Array{Any,1}:
 ([0.0 0.0], [2.0 3.0])
 ([2.0 0.0], [4.0 3.0])
```
"""
function componentgraph(V, copEV, bicon_comps)

	# arrangement of isolated components
	n = size(bicon_comps, 1)
	shells = Array{Lar.Chain, 1}(undef, n)
	boundaries = Array{Lar.ChainOp, 1}(undef, n)
	EVs = Array{Lar.ChainOp, 1}(undef, n)
	# for each component
	for p=1:n
		ev = copEV[sort(bicon_comps[p]), :]
		# computation of 2-cells
		fe = Lar.Arrangement.minimal_2cycles(V, ev)
		# exterior cycle
		global shell_num = Lar.Arrangement.get_external_cycle(V, ev, fe)
		# decompose each fe (co-boundary local to component)
		EVs[p] = ev
		global tokeep = setdiff(1:fe.m, shell_num)
		boundaries[p] = fe[tokeep, :]
		shells[p] = fe[shell_num, :]
	end

#	@show shell_num;
#	@show [SparseArrays.findnz(EVs[k]) for k=1:length(EVs)];
#	@show tokeep;
#	@show [SparseArrays.findnz(boundaries[k]) for k=1:length(boundaries)];
#	@show [SparseArrays.findnz(shells[k]) for k=1:length(shells)];

	# computation of bounding boxes of isolated components
	shell_bboxes = []
	for i in 1:n
		vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
		@show vs_indexes
		push!(shell_bboxes, Lar.bbox(V[vs_indexes, :]))
		@show shell_bboxes
	end
	# computation and reduction of containment graph
	containment_graph = Lar.Arrangement.pre_containment_test(shell_bboxes)
	@show 1,Matrix(containment_graph)
	containment_graph = Lar.Arrangement.prune_containment_graph(n, V, EVs, shells, containment_graph)
	@show 2,Matrix(containment_graph)
	Lar.Arrangement.transitive_reduction!(containment_graph)
	@show 3,Matrix(containment_graph)
	return n, containment_graph, V, EVs, boundaries, shells, shell_bboxes
end

"""
    cleandecomposition(V, copEV, sigma)

!!! NOT WORKING

This function clears the model `(V, copEV)` from all edges outside ``σ``.

This function takes a model `(V, copEV)` and a `Lar.Chain` ``σ`` and
gives back the dropped model and and a vector `todel` of the 2-cells to drop.

---

# Examples
```jldoctest
# Nested Triangles
julia> V = [
		0.0 0.0; 2.0 0.0; 4.0 0.0;
		1.0 1.5; 3.0 1.5; 2.0 3.0;
		2.0 -3.; -2. 3.0; 6.0 3.0
	];

julia> copEV = SparseArrays.sparse(Array{Int8, 2}([
    [1 0 1 0 0 0 0 0 0] #1 -> 1,3
    [1 0 0 0 0 1 0 0 0] #2 -> 1,6
    [0 0 1 0 0 1 0 0 0] #3 -> 3,6
    [0 1 0 1 0 0 0 0 0] #4 -> 2,4
    [0 1 0 0 1 0 0 0 0] #5 -> 2,5
    [0 0 0 1 1 0 0 0 0] #6 -> 4,5
    [0 0 0 0 0 0 1 1 0] #7 -> 7,8
    [0 0 0 0 0 0 1 0 1] #8 -> 7,9
    [0 0 0 0 0 0 0 1 1] #9 -> 8,9
    ]));

julia> σ = SparseArrays.sparse([0; 0; 0; 1; 1; 1; 0; 0; 0]);

todel, V, copEV = Lara.cleandecomposition(V, copEV, convert(Lar.Chain, σ))

Plasm.view(convert(Lar.Points, V'), Lar.cop2lar(copEV));
```
"""
function cleandecomposition(V::Lar.Points, copEV::Lar.ChainOp, sigma::Lar.Chain)
    # Deletes edges outside sigma area
	todel = []
	new_edges = []
	map(i->new_edges=union(new_edges, edge_map[i]), sigma.nzind) # ???
	ev = copEV[new_edges, :]
	for e in 1:copEV.m
		if !(e in new_edges)  # ???  remove if ???

			vidxs = copEV[e, :].nzind
			v1, v2 = map(i->V[vidxs[i], :], [1,2])
			centroid = .5*(v1 + v2)

			if ! Lar.point_in_face(centroid, V, ev)
				push!(todel, e)
			end
		end
	end
	for i in reverse(todel)
		for row in edge_map

			filter!(x->x!=i, row)

			for j in 1:length(row)
				if row[j] > i
					row[j] -= 1
				end
			end
		end
	end
	V, copEV = Lar.delete_edges(todel, V, copEV)

    # biconnected components
    bicon_comps = Lar.Arrangement.biconnected_components(copEV) # -> arrays of edge indices
    @show bicon_comps,0
    if isempty(bicon_comps)
        println("No biconnected components found.")
        if (return_edge_map)
            return (nothing, nothing, nothing, nothing)
        else
            return (nothing, nothing, nothing)
        end
    end

	#remove dangling edges
    edges = sort(union(bicon_comps...))
    todel = sort(setdiff(collect(1:size(copEV,1)), edges))
    for i in reverse(todel)
        for row in edge_map

            filter!(x->x!=i, row)

            for j in 1:length(row)
                if row[j] > i
                    row[j] -= 1
                end
            end
        end
    end
    return todel, V, copEV
end

#--------------------------------------------------------------------------------------------------------------------------------
#								MAIN PIPELINE
#--------------------------------------------------------------------------------------------------------------------------------

"""
	function planar_arrangement_1(V, copEV[, sigma[, return_edge_map[, multiproc]]])

First part of arrangement's algorithmic pipeline.

This function computes the pairwise intersection between each edge of a given 2D cellular complex 1-skeleton.
The computation is speeded up via the evaluation of the Spatial Index. See [`Lar.spaceindex`](@ref).

See also: [`Lar.planar_arrangement`](@ref) for the complete pipeline.

---

# Examples
```jldoctest
julia> EV = [[1, 2], [3, 4], [1, 3], [2, 4], [5, 6], [7, 8], [5, 7], [6, 8]];

julia> V = [
           0.0 0.5 0.0 0.5 0.3 1.0 0.3 1.0;
           0.0 0.0 1.0 1.0 0.5 0.5 1.0 1.0
       ];

julia> W = convert(Lar.Points, V'); # Infering type for W = V'

julia> cop_EV = Lar.coboundary_0(EV::Lar.Cells);

julia> W1, copEV1 = Lar.planar_arrangement_1(W::Lar.Points, cop_EV::Lar.ChainOp)
([0.0 0.0; 0.5 0.0; … ; 1.0 0.5; 1.0 1.0],
  [1 ,  1]  =  1
  [4 ,  1]  =  1
  [1 ,  2]  =  1
  ⋮
  [11,  8]  =  1
  [9 ,  9]  =  1
  [11,  9]  =  1)

julia> EV1 = Lar.cop2lar(copEV1);

julia> V1 = convert(Lar.Points, W1');

julia> Plasm.view(Plasm.numbering(0.1)((V,[[[k] for k=1:size(V,2)], EV])));

julia> Plasm.view(Plasm.numbering(0.1)((V1,[[[k] for k=1:size(V1,2)], EV1])));
```
"""
function planar_arrangement_1(V::Lar.Points, copEV::Lar.ChainOp,
		sigma::Lar.Chain=SparseArrays.spzeros(Int8, 0),
		return_edge_map::Bool=false,
		multiproc::Bool=false)

	# data structures initialization
	edgenum = size(copEV, 1)
	edge_map = Array{Array{Int, 1}, 1}(undef,edgenum)
	rV = Lar.Points(zeros(0, 2))
	rEV = SparseArrays.spzeros(Int8, 0, 0)
	finalcells_num = 0

	# spaceindex computation
	model = (convert(Lar.Points,V'),Lar.cop2lar(copEV))
	bigPI = Lar.spaceindex(model::Lar.LAR)

	# multiprocessing of edge fragmentation
	if (multiproc == true)
		in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
		out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))
		ordered_dict = SortedDict{Int64,Tuple}()
		@async begin
			for i in 1:edgenum
				put!(in_chan,i)
			end
			for p in distributed.workers()
				put!(in_chan,-1)
			end
		end
		for p in distributed.workers()
			@async Base.remote_do(frag_edge_channel, p, in_chan, out_chan, V, copEV, bigPI)
		end
		for i in 1:edgenum
			frag_done_job = take!(out_chan)
			ordered_dict[frag_done_job[1]] = frag_done_job[2]
		end
		for (dkey, dval) in ordered_dict
			i = dkey
			v, ev = dval
			newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))
			edge_map[i] = newedges_nums
			finalcells_num += size(ev, 1)
			rV, rEV = Lar.skel_merge(rV, rEV, v, ev)
		end
	else
	# sequential (iterative) processing of edge fragmentation
		for i in 1:edgenum
			v, ev = Lar.Arrangement.frag_edge(V, copEV, i, bigPI)
			newedges_nums = map(x->x+finalcells_num, collect(1:size(ev, 1)))
			edge_map[i] = newedges_nums
			finalcells_num += size(ev, 1)
			rV, rEV = Lar.skel_merge(rV, rEV, v, ev) # block diagonal ...
		end
	end
	# merging of close vertices and edges (2D congruence)
	V, copEV = rV, rEV
    if return_edge_map
        V, copEV, edge_map = Lar.Arrangement.merge_vertices!(V, copEV, edge_map)
        return V, copEV, edge_map
    end
	V, copEV = Lar.Arrangement.merge_vertices!(V, copEV, edge_map)
	return V, copEV
end

"""
	function planar_arrangement_2(V, copEV, bicon_comps
		[, sigma[, return_edge_map[, multiproc]]]
	)

Second part of arrangement's algorithmic pipeline.

This function is the complete Topological Gift Wrapping (TGW) algorithm that is firstly
locally used in order to decompose the 2-cells and then globally to generate the 3-cells
of the arrangement of the ambient space ``E^3``.

During this process each dangling 2-cell is removed.
Do note that the isolated 1-cells are not removed by this procedure.

See also: [`Lar.planar_arrangement`](@ref) for the complete pipeline.

---

# Examples
```jldoctest
# Triforce
julia> V = [0.0 0.0; 2.0 0.0; 4.0 0.0; 1.0 1.5; 3.0 1.5; 2.0 3.0; 3.0 3.0];
julia> W = convert(Lar.Points, V');
julia> EV = SparseArrays.sparse(Array{Int8, 2}([
    [1 1 0 0 0 0 0] #1 -> 1,2
    [0 1 1 0 0 0 0] #2 -> 2,3
    [1 0 0 1 0 0 0] #3 -> 1,4
    [0 0 0 1 0 1 0] #4 -> 4,6
    [0 0 1 0 1 0 0] #5 -> 3,5
    [0 0 0 0 1 1 0] #6 -> 5,6
    [0 1 0 1 0 0 0] #7 -> 2,4
    [0 1 0 0 1 0 0] #8 -> 2,5
    [0 0 0 1 1 0 0] #9 -> 4,5
    [0 0 0 0 0 1 1]
    ]));

julia> Plasm.view(Plasm.numbering(0.5)(
		(W,[[[k] for k=1:size(W,2)], Lar.cop2lar(EV)])
	));

julia> bicon_comps = Lar.Arrangement.biconnected_components(EV);

julia> V, EV, FE = Lar.Arrangement.planar_arrangement_2(V, EV, bicon_comps);

julia> Plasm.view( Plasm.numbering1(0.5)((V, EV, FE)) );
```
"""
function planar_arrangement_2(V::Lar.Points, copEV::Lar.ChainOp, bicon_comps::Array{Array{Int, 1}, 1},
		sigma::Lar.Chain=spzeros(Int8, 0),
		return_edge_map::Bool=false,
		multiproc::Bool=false)

	# Topological Gift Wrapping
	n, containment_graph, V, EVs, boundaries, shells, shell_bboxes =
		Lar.Arrangement.componentgraph(V, copEV, bicon_comps)
	@show containment_graph
	# only in the context of 3D arrangement
	if sigma.n > 0 ## number of columns of sigma > 0
		todel, V, copEV = Lar.Arrangement.cleandecomposition(V, copEV, sigma)
		V, copEV = Lar.delete_edges(todel, V, copEV) # <----------------------non mi è chiaro perchè non è dentro la fun
	end
	# final shell poset aggregation and FE output
	copEV, FE = Lar.Arrangement.cell_merging(
		n, containment_graph, V, EVs, boundaries, shells, shell_bboxes)
	if (return_edge_map)
		return V, copEV, FE, edge_map
	end
	return V, copEV, FE
end

"""
	planar_arrangement(V::Points, copEV::ChainOp
		[, sigma::Chain[, return_edge_map::Bool[, multiproc::Bool]]]
	)

Compute the arrangement on the given cellular complex 1-skeleton in 2D.
Whole arrangement's algorithmic pipeline.

A cellular complex is arranged when the intersection of every possible pair of cell
of the complex is empty and the union of all the cells is the whole Euclidean space.
The method with no additional arguments specified gives back the full arranged complex
`V`, `EV` and `FE`.

See also:
  - [`Lar.Arrangement.planar_arrangement_1`](@ref)
  - [`Lar.Arrangement.biconnected_components`](@ref)
  - [`Lar.Arrangement.planar_arrangement_2`](@ref)

## Additional arguments:
  - `sigma::Chain`: if specified, `planar_arrangement` will delete from
    the output every edge and face outside this cell. Defaults to an empty cell.
  - `return_edge_map::Bool`: makes the function return also an `edge_map` which
    maps the edges of the imput to the one of the output. Defaults to `false`.
  - `multiproc::Bool`: Runs the computation in parallel mode. Defaults to `false`.

Many examples could be found in Lar Documentation and `examples` directory.
"""
function planar_arrangement( V::Lar.Points, copEV::Lar.ChainOp,
		sigma::Lar.Chain=spzeros(Int8, 0),
		return_edge_map::Bool=false,
		multiproc::Bool=false)

	# edge subdivision
	V, copEV = Lar.Arrangement.planar_arrangement_1(V::Lar.Points, copEV::Lar.ChainOp)
	# biconnected components
	bicon_comps = Lar.Arrangement.biconnected_components(copEV::Lar.ChainOp)
	# 2-complex and containment graph
	V, copEV, copFE = Lar.Arrangement.planar_arrangement_2(V::Lar.Points, copEV::Lar.ChainOp, bicon_comps)
	return V, copEV, copFE
end


text:
# Orthonormalize
ortho_qr(φk) = Matrix(qr(φk).Q)

ext: jl


text:
module test_julia_things

using Test

@test something(nothing, 1) == 1
@test something(nothing, Some(1)) == 1
@test_throws ArgumentError something(nothing)
@test sum([missing, 1]) === missing
@test sum(skipmissing([missing, 1])) == 1

@test nothing === Nothing()
@test missing === Missing()

if VERSION < v"1.3.0-rc"
    @test_throws ArgumentError print(nothing)
end

end # module test_julia_things


text:
"do nothing, this model does not have complex voltage constraints"
function constraint_voltage_dc(pm::_PM.AbstractPowerModel,  n::Int)
end

"""
```
sum(p_dcgrid[a] for a in bus_arcs_dcgrid) + sum(pconv_dc[c] for c in bus_convs_dc) == pd
```
"""
function constraint_power_balance_dc(pm::_PM.AbstractPowerModel, n::Int, i::Int, bus_arcs_dcgrid, bus_convs_dc, pd)
    p_dcgrid = _PM.var(pm, n, :p_dcgrid)
    pconv_dc = _PM.var(pm, n, :pconv_dc)

    JuMP.@constraint(pm.model, sum(p_dcgrid[a] for a in bus_arcs_dcgrid) + sum(pconv_dc[c] for c in bus_convs_dc) == (-pd))
end

"`pconv[i] == pconv`"
function constraint_active_conv_setpoint(pm::_PM.AbstractPowerModel, n::Int, i, pconv)
    pconv_var = _PM.var(pm, n, :pconv_tf_fr, i)

    JuMP.@constraint(pm.model, pconv_var == -pconv)
end

"`qconv[i] == qconv`"
function constraint_reactive_conv_setpoint(pm::_PM.AbstractPowerModel, n::Int, i, qconv)
    qconv_var = _PM.var(pm, n, :qconv_tf_fr, i)

    JuMP.@constraint(pm.model, qconv_var == -qconv)
end

###################### TNEP Constraints ############################
function constraint_voltage_dc_ne(pm::_PM.AbstractPowerModel,  n::Int)
end

function constraint_converter_limit_on_off(pm::_PM.AbstractDCPModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    pconv_tf_fr = _PM.var(pm, n, :pconv_tf_fr_ne)[i]
    pconv_tf_to = _PM.var(pm, n, :pconv_tf_to_ne)[i]
    pconv_pr_fr = _PM.var(pm, n, :pconv_pr_fr_ne)[i]
    z = _PM.var(pm, n, :conv_ne)[i]
    JuMP.@constraint(pm.model,  pconv_ac <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_ac >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_dc <= pmaxdc * z)
    JuMP.@constraint(pm.model,  pconv_dc >= pmindc * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_tf_to <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_to >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr >= pmin * z)

end

function constraint_converter_limit_on_off(pm::_PM.AbstractACPModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax)
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    pconv_tf_fr = _PM.var(pm, n, :pconv_tf_fr_ne)[i]
    pconv_tf_to = _PM.var(pm, n, :pconv_tf_to_ne)[i]
    pconv_pr_fr = _PM.var(pm, n, :pconv_pr_fr_ne)[i]

    qconv_ac = _PM.var(pm, n, :qconv_ac_ne)[i]
    qconv_tf_fr = _PM.var(pm, n, :qconv_tf_fr_ne)[i]
    qconv_tf_to = _PM.var(pm, n, :qconv_tf_to_ne)[i]
    qconv_pr_fr = _PM.var(pm, n, :qconv_pr_fr_ne)[i]
    iconv_ac = _PM.var(pm, n, :iconv_ac_ne)[i]
    vmc = _PM.var(pm, n, :vmc_ne, i)
    vmf = _PM.var(pm, n, :vmf_ne, i)

    z = _PM.var(pm, n, :conv_ne)[i]

    JuMP.@constraint(pm.model,  pconv_ac <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_ac >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_dc <= pmaxdc * z)
    JuMP.@constraint(pm.model,  pconv_dc >= pmindc * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_tf_to <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_to >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr >= pmin * z)
    #
    JuMP.@constraint(pm.model,  qconv_ac <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_ac >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_to <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_to >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr >= qmin * z)
    JuMP.@constraint(pm.model,  iconv_ac <= imax * z )

end

function constraint_converter_limit_on_off(pm::_PM.AbstractBFModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax)
    #converter
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne)[i]
    iconv_ac = _PM.var(pm, n, :iconv_ac_ne)[i]
    iconv_sq = _PM.var(pm, n, :iconv_ac_sq_ne)[i]
    z = _PM.var(pm, n, :conv_ne)[i]
    JuMP.@constraint(pm.model,  pconv_ac <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_ac >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_dc <= pmaxdc * z)
    JuMP.@constraint(pm.model,  pconv_dc >= pmindc * z)
    JuMP.@constraint(pm.model,  qconv_ac <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_ac >= qmin * z)
    JuMP.@constraint(pm.model,  iconv_ac <= imax * z)
    JuMP.@constraint(pm.model,  iconv_sq <= imax^2 * z)
    #transformer
    conv = PowerModels.ref(pm, n, :convdc_ne, i)
    busac_conv = PowerModels.ref(pm, n, :bus, conv["busac_i"])
    w_du = _PM.var(pm, n, :w_du, i)

    pconv_tf_fr = _PM.var(pm, n, :pconv_tf_fr_ne)[i]
    pconv_tf_to = _PM.var(pm, n, :pconv_tf_to_ne)[i]
    qconv_tf_fr = _PM.var(pm, n, :qconv_tf_fr_ne)[i]
    qconv_tf_to = _PM.var(pm, n, :qconv_tf_to_ne)[i]
    itf = _PM.var(pm, n, :itf_sq_ne, i)
    JuMP.@constraint(pm.model,  pconv_tf_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_tf_to <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_to >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_to <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_to >= qmin * z)
    bigM = 2;
    JuMP.@constraint(pm.model,  itf <= (bigM*imax)^2 * z)

    #filter

    #reactor
    pconv_pr_fr = _PM.var(pm, n, :pconv_pr_fr_ne)[i]
    qconv_pr_fr = _PM.var(pm, n, :qconv_pr_fr_ne)[i]
    ipr = _PM.var(pm, n, :irc_sq_ne, i)
    wc = _PM.var(pm, n , :wc_ac_ne,i)
    JuMP.@constraint(pm.model,  pconv_pr_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr >= qmin * z)
    bigM = 2;
    JuMP.@constraint(pm.model,  ipr <= (bigM * imax)^2 * z) #big M = 2
end

function constraint_converter_limit_on_off(pm::_PM.AbstractWModels, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax)
    #converter
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne)[i]
    iconv_ac = _PM.var(pm, n, :iconv_ac_ne)[i]
    iconv_sq = _PM.var(pm, n, :iconv_ac_sq_ne)[i]
    z = _PM.var(pm, n, :conv_ne)[i]
    JuMP.@constraint(pm.model,  pconv_ac <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_ac >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_dc <= pmaxdc * z)
    JuMP.@constraint(pm.model,  pconv_dc >= pmindc * z)
    JuMP.@constraint(pm.model,  qconv_ac <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_ac >= qmin * z)
    JuMP.@constraint(pm.model,  iconv_ac <= imax * z)
    JuMP.@constraint(pm.model,  iconv_sq <= imax^2 * z)

    # transformer
    pconv_tf_fr = _PM.var(pm, n, :pconv_tf_fr_ne)[i]
    pconv_tf_to = _PM.var(pm, n, :pconv_tf_to_ne)[i]
    qconv_tf_fr = _PM.var(pm, n, :qconv_tf_fr_ne)[i]
    qconv_tf_to = _PM.var(pm, n, :qconv_tf_to_ne)[i]
    JuMP.@constraint(pm.model,  pconv_tf_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_tf_to <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_to >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_to <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_to >= qmin * z)
    #filter
    #reactor?
    pconv_pr_fr = _PM.var(pm, n, :pconv_pr_fr_ne)[i]
    qconv_pr_fr = _PM.var(pm, n, :qconv_pr_fr_ne)[i]
    JuMP.@constraint(pm.model,  pconv_pr_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr >= qmin * z)
end

function constraint_converter_limit_on_off(pm::_PM.AbstractWRMModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax)
    #converter
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne)[i]
    iconv_ac = _PM.var(pm, n, :iconv_ac_ne)[i]
    iconv_sq = _PM.var(pm, n, :iconv_ac_sq_ne)[i]
    z = _PM.var(pm, n, :conv_ne)[i]
    JuMP.@constraint(pm.model,  pconv_ac <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_ac >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_dc <= pmaxdc * z)
    JuMP.@constraint(pm.model,  pconv_dc >= pmindc * z)
    JuMP.@constraint(pm.model,  qconv_ac <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_ac >= qmin * z)
    JuMP.@constraint(pm.model,  iconv_ac <= imax * z)
    JuMP.@constraint(pm.model,  iconv_sq <= imax^2 * z)

    # transformer
    pconv_tf_fr = _PM.var(pm, n, :pconv_tf_fr_ne)[i]
    pconv_tf_to = _PM.var(pm, n, :pconv_tf_to_ne)[i]
    qconv_tf_fr = _PM.var(pm, n, :qconv_tf_fr_ne)[i]
    qconv_tf_to = _PM.var(pm, n, :qconv_tf_to_ne)[i]
    JuMP.@constraint(pm.model,  pconv_tf_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_tf_to <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_to >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_to <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_to >= qmin * z)
    #filter
    #reactor?
    pconv_pr_fr = _PM.var(pm, n, :pconv_pr_fr_ne)[i]
    qconv_pr_fr = _PM.var(pm, n, :qconv_pr_fr_ne)[i]
    JuMP.@constraint(pm.model,  pconv_pr_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr >= qmin * z)
end

function constraint_converter_limit_on_off(pm::_PM.AbstractLPACModel, n::Int, i, pmax, pmin, qmax, qmin, pmaxdc, pmindc, imax)
    #converter
    pconv_ac = _PM.var(pm, n, :pconv_ac_ne)[i]
    pconv_dc = _PM.var(pm, n, :pconv_dc_ne)[i]
    qconv_ac = _PM.var(pm, n, :qconv_ac_ne)[i]
    iconv_ac = _PM.var(pm, n, :iconv_ac_ne)[i]
    z = _PM.var(pm, n, :conv_ne)[i]
    JuMP.@constraint(pm.model,  pconv_ac <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_ac >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_dc <= pmaxdc * z)
    JuMP.@constraint(pm.model,  pconv_dc >= pmindc * z)
    JuMP.@constraint(pm.model,  qconv_ac <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_ac >= qmin * z)
    JuMP.@constraint(pm.model,  iconv_ac <= imax * z) #used

    # transformer
    pconv_tf_fr = _PM.var(pm, n, :pconv_tf_fr_ne)[i]
    pconv_tf_to = _PM.var(pm, n, :pconv_tf_to_ne)[i]
    qconv_tf_fr = _PM.var(pm, n, :qconv_tf_fr_ne)[i]
    qconv_tf_to = _PM.var(pm, n, :qconv_tf_to_ne)[i]
    JuMP.@constraint(pm.model,  pconv_tf_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_fr >= pmin * z)
    JuMP.@constraint(pm.model,  pconv_tf_to <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_tf_to >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_fr >= qmin * z)
    JuMP.@constraint(pm.model,  qconv_tf_to <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_tf_to >= qmin * z)
    #filter
    #reactor?
    pconv_pr_fr = _PM.var(pm, n, :pconv_pr_fr_ne)[i]
    qconv_pr_fr = _PM.var(pm, n, :qconv_pr_fr_ne)[i]
    JuMP.@constraint(pm.model,  pconv_pr_fr <= pmax * z)
    JuMP.@constraint(pm.model,  pconv_pr_fr >= pmin * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr <= qmax * z)
    JuMP.@constraint(pm.model,  qconv_pr_fr >= qmin * z)
end

function constraint_branch_limit_on_off(pm::_PM.AbstractBFModel, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin)
    p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx]
    p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx]
    z = _PM.var(pm, n, :branchdc_ne)[i]
    ccm_dcgrid = _PM.var(pm, n, :ccm_dcgrid_ne, i)
    JuMP.@constraint(pm.model,  p_fr <= pmax * z)
    JuMP.@constraint(pm.model,  p_fr >= pmin * z)
    JuMP.@constraint(pm.model,  p_to <= pmax * z)
    JuMP.@constraint(pm.model,  p_to >= pmin * z)
    JuMP.@constraint(pm.model,  ccm_dcgrid <= imax^2 * z)
    JuMP.@constraint(pm.model,  ccm_dcgrid >= imin^2 * z)
end

function constraint_branch_limit_on_off(pm::_PM.AbstractWRModels, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin)
    p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx]
    p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx]
    z = _PM.var(pm, n, :branchdc_ne)[i]
    JuMP.@constraint(pm.model,  p_fr <= pmax * z)
    JuMP.@constraint(pm.model,  p_fr >= pmin * z)
    JuMP.@constraint(pm.model,  p_to <= pmax * z)
    JuMP.@constraint(pm.model,  p_to >= pmin * z)
end

function constraint_branch_limit_on_off(pm::_PM.AbstractPowerModel, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin)
    p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx]
    p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx]
    z = _PM.var(pm, n, :branchdc_ne)[i]

    JuMP.@constraint(pm.model,  p_fr <= pmax * z)
    JuMP.@constraint(pm.model,  p_fr >= pmin * z)
    JuMP.@constraint(pm.model,  p_to <= pmax * z)
    JuMP.@constraint(pm.model,  p_to >= pmin * z)
end

function constraint_branch_limit_on_off(pm::_PM.AbstractACPModel, n::Int, i, f_idx, t_idx, pmax, pmin, imax, imin)
    p_fr = _PM.var(pm, n, :p_dcgrid_ne)[f_idx]
    p_to = _PM.var(pm, n, :p_dcgrid_ne)[t_idx]
    z = _PM.var(pm, n, :branchdc_ne)[i]

    JuMP.@constraint(pm.model,  p_fr <= pmax * z)
    JuMP.@constraint(pm.model,  p_fr >= pmin * z)
    JuMP.@constraint(pm.model,  p_to <= pmax * z)
    JuMP.@constraint(pm.model,  p_to >= pmin * z)
end

function constraint_candidate_converters_mp(pm::_PM.AbstractPowerModel, n::Int, i::Int)
    z = _PM.var(pm, n, :conv_ne, i)
    z_1 = _PM.var(pm, n-1, :conv_ne, i)

    JuMP.@constraint(pm.model,  z == z_1)
end

function constraint_candidate_dcbranches_mp(pm::_PM.AbstractPowerModel, n::Int, i::Int)
    z = _PM.var(pm, n, :branchdc_ne, i)
    z_1 = _PM.var(pm, n-1, :branchdc_ne, i)

    JuMP.@constraint(pm.model,  z == z_1)
end

function constraint_candidate_acbranches_mp(pm::_PM.AbstractPowerModel, n::Int, i::Int)
    z = _PM.var(pm, n, :branch_ne, i)
    z_1 = _PM.var(pm, n-1, :branch_ne, i)

    JuMP.@constraint(pm.model,  z == z_1)
end

function constraint_power_balance_dc_dcne(pm::_PM.AbstractPowerModel, n::Int, i::Int, bus_arcs_dcgrid, bus_arcs_dcgrid_ne, bus_convs_dc, bus_convs_dc_ne, pd)
    p_dcgrid = _PM.var(pm, n, :p_dcgrid)
    p_dcgrid_ne = _PM.var(pm, n, :p_dcgrid_ne)
    pconv_dc = _PM.var(pm, n, :pconv_dc)
    pconv_dc_ne = _PM.var(pm, n, :pconv_dc_ne)

    JuMP.@constraint(pm.model, sum(p_dcgrid[a] for a in bus_arcs_dcgrid) + sum(p_dcgrid_ne[a] for a in bus_arcs_dcgrid_ne) + sum(pconv_dc[c] for c in bus_convs_dc) + sum(pconv_dc_ne[c] for c in bus_convs_dc_ne)  == (-pd))
end

function constraint_power_balance_dcne_dcne(pm::_PM.AbstractPowerModel, n::Int, i::Int, bus_arcs_dcgrid_ne, bus_ne_convs_dc_ne, pd_ne)
    p_dcgrid_ne = _PM.var(pm, n, :p_dcgrid_ne)
    pconv_dc_ne = _PM.var(pm, n, :pconv_dc_ne)
    xb = _PM.var(pm, n, :branchdc_ne)
    xc = _PM.var(pm, n, :conv_ne)
    JuMP.@constraint(pm.model, sum(p_dcgrid_ne[a] for a in bus_arcs_dcgrid_ne) + sum(pconv_dc_ne[c] for c in bus_ne_convs_dc_ne)  == (-pd_ne))
end


text:
module Lexicon

import Docile.Interface:

    parsedocs,
    macroname,
    name

import Base:

    start,
    next,
    done,
    length,
    push!,
    run,
    writemime,
    ==

using

    Base.Meta,
    Compat,
    Docile,
    Docile.Interface

export

    @query,
    query,
    save,
    update!,

    doctest,
    failed,
    passed,
    skipped,
    EachEntry,
    Config,
    Index

@document

# Conditional importing of the `Markdown` module.
if VERSION < v"0.4-dev+1488"
    include("../deps/Markdown/src/Markdown.jl")
    import .Markdown
end

include("compat.jl")
include("query.jl")
include("render.jl")
include("doctest.jl")
include("filtering.jl")

__init__() = setup_help() # Hook into the REPL's `?`.

end # module


text:
struct ConvergenceHistory{T, R}
    isconverged::Bool
    threshold::T
    mvps::Int
    residuals::R
end

##############################################################################
## LSMR
##
## Minimize ||Ax-b||^2 + λ^2 ||x||^2
##
## Adapted from the BSD-licensed Matlab implementation at
## http://web.stanford.edu/group/SOL/software/lsmr/
##
## A is a StridedVecOrMat or anything that implements
## mul!
## eltype(A)
## size(A)
## (this includes SparseMatrixCSC)
## x, v, h, hbar are AbstractVectors or anything that implements
## norm(x)
## copyto!(x1, x2)
## rmul!(x, α)
## axpy!(α, x1, x2)
## similar(x, T)
## length(x)
## b is an AbstractVector or anything that implements
## eltype(b)
## norm(b)
## copyto!(x1, x2)
## fill!(b, α)
## rmul!(b, α)
## similar(b, T)
## length(b)

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

## Arguments:
## x is initial guess for x0. Transformed in place to the solution.
## b is also transformed in place
## v, h, hbar are storage arrays of length size(A, 2)
function lsmr!(x, A, b, v, h, hbar;
    atol::Number = 1e-6, btol::Number = 1e-6, conlim::Number = 1e8,
    maxiter::Integer = max(size(A,1), size(A,2)), λ::Number = 0)

    # Sanity-checking
    m = size(A, 1)
    n = size(A, 2)
    length(x) == n || error("x has length $(length(x)) but should have length $n")
    length(v) == n || error("v has length $(length(v)) but should have length $n")
    length(h) == n || error("h has length $(length(h)) but should have length $n")
    length(hbar) == n || error("hbar has length $(length(hbar)) but should have length $n")
    length(b) == m || error("b has length $(length(b)) but should have length $m")

    T = Base.promote_op(/, eltype(b), eltype(A))
    Tr = real(T)
    normrs = Tr[]
    normArs = Tr[]
    conlim > 0 ? ctol = convert(Tr, inv(conlim)) : ctol = zero(Tr)
    # form the first vectors u and v (satisfy  β*u = b,  α*v = A'u)
    u = mul!(b, A, x, -1, 1)
    β = norm(u)
    β > 0 && rmul!(u, inv(β))
    mul!(v, A', u, 1, 0)
    α = norm(v)
    α > 0 && rmul!(v, inv(α))

    # Initialize variables for 1st iteration.
    ζbar = α * β
    αbar = α
    ρ = one(Tr)
    ρbar = one(Tr)
    cbar = one(Tr)
    sbar = zero(Tr)

    copyto!(h, v)
    fill!(hbar, zero(Tr))

    # Initialize variables for estimation of ||r||.
    βdd = β
    βd = zero(Tr)
    ρdold = one(Tr)
    τtildeold = zero(Tr)
    θtilde  = zero(Tr)
    ζ = zero(Tr)
    d = zero(Tr)

    # Initialize variables for estimation of ||A|| and cond(A).
    normA, condA, normx = -one(Tr), -one(Tr), -one(Tr)
    normA2 = abs2(α)
    maxrbar = zero(Tr)
    minrbar = 1e100

    # Items for use in stopping rules.
    normb = β
    istop = 0
    normr = β
    normAr = α * β
    tests = Tuple{Tr, Tr, Tr}[]
    iter = 0
    # Exit if b = 0 or A'b = 0.
    if normAr != 0
        while iter < maxiter
            iter += 1
            mul!(u, A, v, 1, -α)
            β = norm(u)
            if β > 0
                rmul!(u, inv(β))
                mul!(v, A', u, 1, -β)
                α = norm(v)
                α > 0 && rmul!(v, inv(α))
            end

            # Construct rotation Qhat_{k,2k+1}.
            αhat = sqrt(abs2(αbar) + abs2(λ))
            chat = αbar / αhat
            shat = λ / αhat

            # Use a plane rotation (Q_i) to turn B_i to R_i.
            ρold = ρ
            ρ = sqrt(abs2(αhat) + abs2(β))
            c = αhat / ρ
            s = β / ρ
            θnew = s * α
            αbar = c * α

            # Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar.
            ρbarold = ρbar
            ζold = ζ
            θbar = sbar * ρ
            ρtemp = cbar * ρ
            ρbar = sqrt(abs2(cbar * ρ) + abs2(θnew))
            cbar = cbar * ρ / ρbar
            sbar = θnew / ρbar
            ζ = cbar * ζbar
            ζbar = - sbar * ζbar

            # Update h, h_hat, x.
            rmul!(hbar, - θbar * ρ / (ρold * ρbarold))
            axpy!(1, h, hbar)
            axpy!(ζ / (ρ * ρbar), hbar, x)
            rmul!(h, - θnew / ρ)
            axpy!(1, v, h)

            ##############################################################################
            ##
            ## Estimate of ||r||
            ##
            ##############################################################################

            # Apply rotation Qhat_{k,2k+1}.
            βacute = chat * βdd
            βcheck = - shat * βdd

            # Apply rotation Q_{k,k+1}.
            βhat = c * βacute
            βdd = - s * βacute

            # Apply rotation Qtilde_{k-1}.
            θtildeold = θtilde
            ρtildeold = sqrt(abs2(ρdold) + abs2(θbar))
            ctildeold = ρdold / ρtildeold
            stildeold = θbar / ρtildeold
            θtilde = stildeold * ρbar
            ρdold = ctildeold * ρbar
            βd = - stildeold * βd + ctildeold * βhat

            τtildeold = (ζold - θtildeold * τtildeold) / ρtildeold
            τd = (ζ - θtilde * τtildeold) / ρdold
            d  = d + abs2(βcheck)
            normr = sqrt(d + abs2(βd - τd) + abs2(βdd))

            # Estimate ||A||.
            normA2 = normA2 + abs2(β)
            normA  = sqrt(normA2)
            normA2 = normA2 + abs2(α)

            # Estimate cond(A).
            maxrbar = max(maxrbar, ρbarold)
            if iter > 1
                minrbar = min(minrbar, ρbarold)
            end
            condA = max(maxrbar, ρtemp) / min(minrbar, ρtemp)

            ##############################################################################
            ##
            ## Test for convergence
            ##
            ##############################################################################

            # Compute norms for convergence testing.
            normAr  = abs(ζbar)
            normx = norm(x)

            # Now use these norms to estimate certain other quantities,
            # some of which will be small near a solution.
            test1 = normr / normb
            test2 = normAr / (normA * normr)
            test3 = inv(condA)
            push!(tests, (test1, test2, test3))

            t1 = test1 / (one(Tr) + normA * normx / normb)
            rtol = btol + atol * normA * normx / normb
            # The following tests guard against extremely small values of
            # atol, btol or ctol.  (The user may have set any or all of
            # the parameters atol, btol, conlim  to 0.)
            # The effect is equivalent to the normAl tests using
            # atol = eps,  btol = eps,  conlim = 1/eps.
            if iter >= maxiter istop = 7; break end
            if 1 + test3 <= 1 istop = 6; break end
            if 1 + test2 <= 1 istop = 5; break end
            if 1 + t1 <= 1 istop = 4; break end
            # Allow for tolerances set by the user.
            if test3 <= ctol istop = 3; break end
            if test2 <= atol istop = 2; break end
            if test1 <= rtol  istop = 1; break end
        end
    end
    converged = istop ∉ (3, 6, 7)
    tol = (atol, btol, ctol)
    ch = ConvergenceHistory(converged, tol, 2 * iter, tests)
    return x, ch
end


text:
using Pseudospectra, Test

@testset "Switch-method" begin
    # 400x400 real sparse
    A = Pseudospectra.convdiff_fd(10)

    # first direct, then ARPACK
    opts = Dict{Symbol,Any}(:direct => true,:npts=>24,:ax => [-10,100,-50,50])
    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
    ps_data_dir = deepcopy(ps_data)

    Pseudospectra.set_method!(ps_data,false)
    opts = Dict{Symbol,Any}(:direct => false, :ARPACK_plots => false,:npts=>40,
                            :ax => [30,100,-50,50])
    opts[:arpack_opts] = ArpackOptions{eltype(A)}(which=:SR,nev=20,ncv=60)

    driver!(ps_data,opts,gs,revise_method=true)
    @test iscomputed(ps_data)

    # first ARPACK, then direct
    opts = Dict{Symbol,Any}(:ARPACK_plots => false,:npts=>40,
                            :ax => [30,100,-50,50])
    opts[:arpack_opts] = ArpackOptions{eltype(A)}(which=:SR,nev=20,ncv=60)
    ps_data = new_matrix(A,opts)
    driver!(ps_data,opts,gs)
    @test iscomputed(ps_data)
    ps_data_arpack = deepcopy(ps_data)

    Pseudospectra.set_method!(ps_data,true)

    opts = Dict{Symbol,Any}(:direct => true,:npts=>24,:ax => [-10,100,-50,50])
    driver!(ps_data,opts,gs,revise_method=true)
    @test iscomputed(ps_data)

    # first ARPACK, then with new opts
    ps_data = deepcopy(ps_data_arpack)

    opts = Dict{Symbol,Any}(:ARPACK_plots => false,
                            :npts=>24,:ax => [-10,100,-50,50])
    opts[:arpack_opts] = ArpackOptions{eltype(A)}(which=:SR,nev=40,ncv=80)

    driver!(ps_data,opts,gs,revise_method=true)
    @test iscomputed(ps_data)
    # CHECKME: if 41 are computed, are all saved?
    @test length(ps_data.ps_dict[:ews]) >= 40
end


text:
"""
This module exports
PFparams, pfStateEstimator

`pfStateEstimator(meas, force, T_RB_TF, search_traj, T_TF_M, params::PFparams)`
Run state estimator based on particle filter\n
`Npart`: number of particles\n
`meas`: measurement from laser sensor\n
`T_RB_TF`: position of the robot tool flange from FK, in R4x4\n
`search_traj`: nominal trajectory in RTx3
"""
module PFstateEstimator
using Robotlib
# using Debug
using Plots
# using Devectorize
include("pf_funs.jl")
using StatsBase # logsumexp

export PFparams, PFtrace, PFstate, pfStateEstimator, init_pf

"""
σW::Vector{T}       # noise in state update\n
σV::Vector{T}       # Noise in sensor reading\n
σVFK::Vector{T}     # Noise in FK\n
Npart::Int              # Number of particles\n
NpartInit::Int          # Increasing number of particles with this factor during initialization
h::Float64              # Sample time\n
force_filter::Float64
"""
type PFparams{T}
    σW::Vector{T}       # noise in state update
    σV::Vector{T}       # Noise in sensor reading
    σVFK::Vector{T}     # Noise in FK
    Npart::Int              # Number of particles
    NpartInit::Int          # Increasing number of particles with this factor during initialization
    h::Float64              # Sample time
    force_filter::Float64   #
end

type PFtrace
    Neff::Float64
    resample::Bool
    x::Array{Float64,3}
    w::Array{Float64,1}
    reset::Bool
    resample_i::AbstractVector{Int}
    runMeasUpdate::Bool
end

type PFstate
    x::Array{Float64,3}
    w::Vector{Float64}
    expw::Vector{Float64}
    T_RB_TF_last::Matrix{Float64}
    meas_last::VecOrMat{Float64}
    force_last::Vector{Float64}
    force_filt::Vector{Float64}
    t::Int
end
# type PFtrace
# end

function init_pf(T_RB_TF, Npart, σW, meas, T_TF_M)
    # m_TF            = (T_RB_TF[1:3,1:3]*T_TF_M[1:3,1:3,1])*[meas[1:2,1];0]
    x_init          = T_RB_TF
    # x_init[1:3,4]  += 0.5*m_TF
    x,w,expw        = initParticles(x_init,Npart,σW)
    T_RB_TF_last = T_RB_TF
    meas_last       = zeros(meas)
    force_last      = zeros(6)
    force_filt      = zeros(6)
    return PFstate(x,w,expw,T_RB_TF_last,meas_last, force_last, force_filt,1)
end

# """
# `pfStateEstimator(meas, force, T_RB_TF, search_traj, T_TF_M, params::PFparams)`
# Run state estimator based on particle filter\n
# `Npart`: number of particles\n
# `meas`: measurement from laser sensor\n
# `T_RB_TF`: position of the robot tool flange from FK, in R4x4\n
# `search_traj`: nominal trajectory in RT×3
# """

function pfStateEstimator(s::PFstate, meas, force, T_RB_TF, traj_nom, search_traj, T_TF_M, params::PFparams; savetrace=false)

    x               = s.x
    w               = s.w
    expw            = s.expw
    T_RB_TF_last    = s.T_RB_TF_last
    meas_last       = s.meas_last
    force_last      = s.force_last
    force_filt      = s.force_filt
    σW              = params.σW
    σV              = params.σV
    σVFK            = params.σVFK
    Npart           = params.Npart
    NpartInit       = params.NpartInit
    h               = params.h
    force_filter    = params.force_filter
    #     m_TF = T_TF_M\[meas';0;1]

    force_filt = force_filter/h*abs(force-force_last) + (1-force_filter)*force_filt
    runMeasUpdate = any(meas .!= 0)

    # Resample--------------------------------------------------------------

    Neff = (1/sum(expw.^2))
    resample = Neff < Npart/2 || s.t == 2

    if resample
        j   = resample_systematic_exp(expw,Npart)
        x[:,:,1:Npart] = x[:,:,j]
        w[1:Npart]  = w[j]
        #x = x[:,:,j]
        #w = w[j]
        if s.t == 2
            x = x[:,:,1:Npart]
            w = w[1:Npart]
            expw = expw[1:Npart]
            s.x = x
            s.w = w
            s.expw = expw
        end
    else
        j = 1:size(x,3)
    end

    # Time update-----------------------------------------------------------
    increment   = T_RB_TF*trinv(T_RB_TF_last)
    # σWt     = (increment[1:3,4] + √(h)).*σW[1] + √(h)*σW[3]*force_filt[1:3] # TODO: This line is ad-hoc
    σWt         = σW[1].*ones(3)
    timeUpdate!(x,increment,σWt,σW[2],h)

    # Meas update-----------------------------------------------------------
    if runMeasUpdate
        for i = 1:size(T_TF_M,3)
            x_μ2, mx_μ, index, γ   = distToSeam(meas[1:2,i], search_traj[i], x, T_TF_M[:,:,i], tries = 1)
            wU      = measUpdate(mx_μ, meas[3,i], x, traj_nom, σV, index, γ, T_TF_M[:,:,i])
            # wU[wU .< -12] = -8 # Robustify against outliers
            # wU[(meas[1,i] != 0) & (x_μ2 .== 0)] -= 1 # If a measurement is obtained, but the seam not found for a particle, punish that particle
            # do NOT add opposite kind of error, no sensor will not work then
            w    += wU
        end
    end

    w    += measFKUpdate(x,T_RB_TF,force,σVFK)
    w,expw[:]  = normalizeWeights!(w)

    # # Reset the particle filter---------------------------------------------
    reset = any(isnan(w)) || any(isinf(w)) #|| (1/sum(expw.^2)) < 2
    # if reset
    #     m_TF            = (T_RB_TF[1:3,1:3]*T_TF_M[1:3,1:3,1])'*[meas[1:2,1];0]
    #     x_init          = T_RB_TF
    #     x_init[1:3,4]  -= m_TF
    #     x, w, expw      = initParticles(x_init,Npart,σW)
    # end

    # Calculate state estimate----------------------------------------------
    x̂ = zeros(3,4)
    for i = 1:size(x,3)
        x̂[1:3,1:4] += expw[i].*x[1:3,1:4,i] # calculate mean particle
    end
    toOrthoNormal!(x̂)

    # x̂ = x[:,:,findmax(w)[2]]
    x̂ = [x̂; 0 0 0 1]

    # Calculate position uncertainty----------------------------------------
    se      = zeros(3)
    for i = 1:size(x,3)
        se += expw[i].*(x[1:3,4,i]-x̂[1:3,4]).^2
    end
    se              = √(se)

    # Update filter states--------------------------------------------------
    T_RB_TF_last[:,:] = T_RB_TF
    meas_last[:]      = meas
    s.t              += 1
    savetrace && return x̂, se, s, PFtrace(Neff,resample, deepcopy(x), deepcopy(expw), reset, j, runMeasUpdate)
    return x̂, se, s
end

function normalizeWeights!(w)
    w -= logsumexp(w)
    expw = exp(w)
    return w, expw
end

function initParticles(x_init,Npart,σW)
    x = repeat(x_init,inner=[1,1],outer=[1,1,Npart])
    x[1:3,4,:] = squeeze(x[1:3,4,:],2) + 5*σW[1]*randn(3,Npart)
    I = eye(3)
    for i = 1:Npart
        k = randn(3)
        k /= norm(k) # Axis of rotation
        theta = 5*σW[2]*randn()
        kx = skew(k)
        R = I + sin(theta)*kx + (1-cos(theta))*(k*k'-I)
        x[1:3,1:3,i] = R*x[1:3,1:3,i]
    end
    expw = 1/Npart*ones(Npart)
    w = log(1/Npart)*ones(Npart)
    return x,w,expw
end

function orthoNormalize(x,j)
    x[:,:,1] = toOrthoNormal(x[:,:,1])
    for i = 2:length(j)
        if j[i-1] == j[i]
            continue
        end
        x[:,:,j[i]] = toOrthoNormal(x[:,:,j[i]])
    end
    return x
end
end


text:
function prove_periodic_orbit(filename::AbstractString)
    println("Starting proof for ", filename, ":")
    f = load(filename)
    return _prove_periodic_orbit(f["x"], f["C₁"], f["C₂"], f["ν"], f["R"])
end

function _prove_periodic_orbit(x, C₁::Interval{Float64}, C₂::Interval{Float64}, ν::Interval{Float64}, R::Float64)
    n = order(space(space(x)[1]))

    # ensures the numerical approximation satisfies the complex conjugacy symmetry

    x₁, x₂ = component(x, 1), component(x, 2)
    for x₁ⱼ ∈ eachcomponent(x₁)
        x₁ⱼ .= (x₁ⱼ .+ conj.(view(x₁ⱼ, n:-1:-n))) ./ 2
        x₁ⱼ[0] = real(x₁ⱼ[0])
    end
    x₂ .= real.(x₂)

    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(x₁)
    γ = x₂[1]
    s₁, s₂, s₃, s₄, s₅ = sum(Σ₊), sum(Σ₋), sum(N₁), sum(N₂), sum(N₃)

    #

    print("    (i)\n        • R = ", R, " and ν = "); showfull(ν); print("\n")

    space_F = Fourier(2n, 1.0)^5 × ParameterSpace()^3
    F = Sequence(space_F, Vector{Complex{Interval{Float64}}}(undef, dimension(space_F)))
    F!(F, x, C₁, C₂)

    domain_DF = Fourier(3n, 1.0)^5 × ParameterSpace()^3
    DF = LinearOperator(domain_DF, space_F, Matrix{Complex{Interval{Float64}}}(undef, dimension(space_F), dimension(domain_DF)))
    DF!(DF, x, C₁, C₂)

    A = inv(mid.(project(DF, space_F, space_F)))

    X_Fourier = ℓ¹(GeometricWeight(ν))
    X_cartesian_Fourier = NormedCartesianSpace(X_Fourier, ℓ∞())
    X = NormedCartesianSpace((X_cartesian_Fourier, NormedCartesianSpace(ℓ∞(), ℓ∞())), ℓ∞())

    opnorm_A = opnorm(Interval.(A), X)
    bound_tail_A = inv(Interval(2n+1))

    Y = norm(A * F, X)
    print("        • Y = "); showfull(Y); print("\n")

    Z₁ = opnorm(A * DF - I, X) +
        bound_tail_A * opnorm_Df(x, X_Fourier) +
        opnorm_A * inv(ν^(3n+1)) * max(
            1,
            2(abs(s₁) + abs(s₂) + 3(abs(s₃) + abs(s₄) + abs(s₅))),
            abs(s₄)*abs(s₅) + abs(s₃)*abs(s₅) + abs(s₃)*abs(s₄)
            )
    print("        • Z₁ = "); showfull(Z₁); print("\n")

    Z₂ = (opnorm_A + bound_tail_A) *
        (opnorm_∂Σ₊DF(x, R, X_Fourier) + opnorm_∂Σ₋DF(x, R, X_Fourier) +
        opnorm_∂N₁DF(x, R, X_Fourier) + opnorm_∂N₂DF(x, R, X_Fourier) + opnorm_∂N₃DF(x, R, X_Fourier) +
        opnorm_∂γDF(x, R, X_Fourier) + 1)
    print("        • Z₂ = "); showfull(Z₂); print("\n")

    ie = interval_of_existence(Y, Z₁, Z₂, R, C²Condition())
    print("        • interval of existence = "); showfull(ie); print("\n")

    #

    if isempty(ie)
        println("    => (i) failed")
    elseif _prove_periodic_orbit_ii(γ, Σ₊, Σ₋, N₁, N₃, inf(ie), R, X_Fourier)
        println("    => (i) and (ii) succeeded")
    else
        println("    => (i) succeeded and (ii) failed")
    end

    return ie
end

function _prove_periodic_orbit(x, C₁::Interval{BigFloat}, C₂::Interval{BigFloat}, ν::Interval{BigFloat}, R::BigFloat)
    n = order(space(space(x)[1]))

    # ensures the numerical approximation satisfies the complex conjugacy symmetry

    x₁, x₂ = component(x, 1), component(x, 2)
    for x₁ⱼ ∈ eachcomponent(x₁)
        x₁ⱼ .= (x₁ⱼ .+ conj.(view(x₁ⱼ, n:-1:-n))) ./ 2
        x₁ⱼ[0] = real(x₁ⱼ[0])
    end
    x₂ .= real.(x₂)

    # set some variables in interval with lesser precision

    v = coefficients(x)
    x_₆₄ = Sequence(
        space(x),
        complex.(
            Interval.(Float64.(inf.(real.(v)), RoundDown), Float64.(sup.(real.(v)), RoundUp)),
            Interval.(Float64.(inf.(imag.(v)), RoundDown), Float64.(sup.(imag.(v)), RoundUp))
            )
        )
    C₁_₆₄ = Interval(Float64(inf(C₁), RoundDown), Float64(sup(C₁), RoundUp))
    C₂_₆₄ = Interval(Float64(inf(C₂), RoundDown), Float64(sup(C₂), RoundUp))
    ν_₆₄ = Interval(Float64(inf(ν), RoundDown), Float64(sup(ν), RoundUp))
    R_₆₄ = Interval(Float64(R, RoundDown), Float64(R, RoundUp))
    Σ₊_₆₄, Σ₋_₆₄, N₁_₆₄, N₂_₆₄, N₃_₆₄ = eachcomponent(component(x_₆₄, 1))
    γ_₆₄ = component(x_₆₄, 2)[1]
    s₁_₆₄, s₂_₆₄, s₃_₆₄, s₄_₆₄, s₅_₆₄ = sum(Σ₊_₆₄), sum(Σ₋_₆₄), sum(N₁_₆₄), sum(N₂_₆₄), sum(N₃_₆₄)

    #

    print("    (i)\n        • R = ", R, " and ν = "); showfull(ν); print("\n")

    space_F = Fourier(2n, 1.0)^5 × ParameterSpace()^3
    F = Sequence(space_F, Vector{Complex{Interval{BigFloat}}}(undef, dimension(space_F)))
    F!(F, x, C₁, C₂)

    domain_DF = Fourier(3n, 1.0)^5 × ParameterSpace()^3
    DF_₆₄ = LinearOperator(domain_DF, space_F, Matrix{Complex{Interval{Float64}}}(undef, dimension(space_F), dimension(domain_DF)))
    DF!(DF_₆₄, x_₆₄, C₁_₆₄, C₂_₆₄)

    A_₆₄ = inv(mid.(project(DF_₆₄, space_F, space_F)))

    X_Fourier_₆₄ = ℓ¹(GeometricWeight(ν_₆₄))
    X_₆₄ = NormedCartesianSpace((NormedCartesianSpace(X_Fourier_₆₄, ℓ∞()), NormedCartesianSpace(ℓ∞(), ℓ∞())), ℓ∞())
    X_Fourier = ℓ¹(GeometricWeight(ν))
    X = NormedCartesianSpace((NormedCartesianSpace(X_Fourier, ℓ∞()), NormedCartesianSpace(ℓ∞(), ℓ∞())), ℓ∞())

    opnorm_A_₆₄ = opnorm(Interval.(A_₆₄), X_₆₄)
    bound_tail_A_₆₄ = inv(Interval(2n+1))

    Y = norm(A_₆₄ * F, X)
    print("        • Y = "); showfull(Y); print("\n")

    Z₁_₆₄ = opnorm(A_₆₄ * DF_₆₄ - I, X_₆₄) +
        bound_tail_A_₆₄ * opnorm_Df(x_₆₄, X_Fourier_₆₄) +
        opnorm_A_₆₄ * inv(ν_₆₄^(3n+1)) * max(
            1,
            2(abs(s₁_₆₄) + abs(s₂_₆₄) + 3(abs(s₃_₆₄) + abs(s₄_₆₄) + abs(s₅_₆₄))),
            abs(s₄_₆₄)*abs(s₅_₆₄) + abs(s₃_₆₄)*abs(s₅_₆₄) + abs(s₃_₆₄)*abs(s₄_₆₄)
            )
    Z₁ = Interval(BigFloat(inf(Z₁_₆₄), RoundDown), BigFloat(sup(Z₁_₆₄), RoundUp))
    print("        • Z₁ = "); showfull(Z₁); print("\n")

    Z₂_₆₄ = (opnorm_A_₆₄ + bound_tail_A_₆₄) *
        (opnorm_∂Σ₊DF(x_₆₄, R_₆₄, X_Fourier_₆₄) + opnorm_∂Σ₋DF(x_₆₄, R_₆₄, X_Fourier_₆₄) +
        opnorm_∂N₁DF(x_₆₄, R_₆₄, X_Fourier_₆₄) + opnorm_∂N₂DF(x_₆₄, R_₆₄, X_Fourier_₆₄) + opnorm_∂N₃DF(x_₆₄, R_₆₄, X_Fourier_₆₄) +
        opnorm_∂γDF(x_₆₄, R_₆₄, X_Fourier_₆₄) + 1)
    Z₂ = Interval(BigFloat(inf(Z₂_₆₄), RoundDown), BigFloat(sup(Z₂_₆₄), RoundUp))
    print("        • Z₂ = "); showfull(Z₂); print("\n")

    ie = interval_of_existence(Y, Z₁, Z₂, R, C²Condition())
    print("        • interval of existence = "); showfull(ie); print("\n")

    #

    if isempty(ie)
        println("    => (i) failed")
    elseif _prove_periodic_orbit_ii(γ_₆₄, Σ₊_₆₄, Σ₋_₆₄, N₁_₆₄, N₃_₆₄, Float64(inf(ie), RoundUp), R_₆₄, X_Fourier_₆₄)
        println("    => (i) and (ii) succeeded")
    else
        println("    => (i) succeeded and (ii) failed")
    end

    return ie
end

function _prove_periodic_orbit_ii(γ, Σ₊, Σ₋, N₁, N₃, inf_r₀, R, X_Fourier)
    println("    (ii)")
    check_1 = iszero(imag(γ)) & (real(γ) > R)
    println("        • check 1: γ ∈ (0, +∞) -> ", check_1)
    check_2 = check_not_identically_zero(Σ₊, inf_r₀)
    println("        • check 2: Σ₊ is not identically zero in [0, 2π) -> ", check_2)
    check_3 = check_not_identically_zero(N₁, inf_r₀)
    println("        • check 3: N₁ is not identically zero in [0, 2π) -> ", check_3)
    check_4 = check_not_identically_zero(N₃, inf_r₀)
    println("        • check 4: N₃ is not identically zero in [0, 2π) -> ", check_4)
    bound_DN₁ = 4(abs(γ) + inf_r₀) * (norm(Σ₊, X_Fourier) + inf_r₀) * (norm(N₁, X_Fourier) + inf_r₀)
    check_5 = check_no_zeros(N₁, bound_DN₁, inf_r₀, "N₁", "5")
    println("        • check 5: N₁ has no zeros in [0, 2π) -> ", check_5)
    bound_DN₃ = 2(abs(γ) + inf_r₀) * (norm(N₃, X_Fourier) + inf_r₀) * (norm(Σ₊, X_Fourier) + inf_r₀ + sqrt(Interval(3.0))*(norm(Σ₋, X_Fourier) + inf_r₀))
    check_6 = check_no_zeros(N₃, bound_DN₃, inf_r₀, "N₃", "6")
    println("        • check 6: N₃ has no zeros in [0, 2π) -> ", check_6)
    return check_1 & check_2 & check_3 & check_4 & check_5 & check_6
end

#

function F!(F, x, C₁, C₂)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)
    s₁, s₂, s₃, s₄, s₅ = sum(Σ₊), sum(Σ₋), sum(N₁), sum(N₂), sum(N₃)

    F₁, F₂, F₃, F₄, F₅ = eachcomponent(component(F, 1))
    F₆ = component(F, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))
    𝒟 = Derivative(1)

    project!(F₁, 2γ * ((N₂ - N₃) ^ 2 - N₁ * (2N₁ - N₂ - N₃)) + η₁ * Σ₊ - 𝒟(Σ₊))
    project!(F₂, 2γ*sqrt3 * (N₃ - N₂) * (N₁ - N₂ - N₃) - 𝒟(Σ₋))
    project!(F₃, 4γ * (Σ₊ * N₁) - 𝒟(N₁))
    project!(F₄, -2γ * N₂ * (Σ₊ + sqrt3 * Σ₋) + η₂ - 𝒟(N₂))
    project!(F₅, -2γ * N₃ * (Σ₊ - sqrt3 * Σ₋) - 𝒟(N₃))

    F₆[1] = s₁
    F₆[2] = s₁^2 + s₂^2 + s₃^2 + s₄^2 + s₅^2 - 2s₃*s₄ - 2s₃*s₅ - 2s₄*s₅ - C₁
    F₆[3] = s₃*s₄*s₅ - C₂

    return F
end

function DF!(DF, x, C₁, C₂)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)
    s₁, s₂, s₃, s₄, s₅ = sum(Σ₊), sum(Σ₋), sum(N₁), sum(N₂), sum(N₃)

    DF .= zero(eltype(DF))
    M₁₁ = component(DF, 1, 1)
    ∂Σ₊F₁, ∂N₁F₁, ∂N₂F₁, ∂N₃F₁ = component(M₁₁, 1, 1), component(M₁₁, 1, 3), component(M₁₁, 1, 4), component(M₁₁, 1, 5)
    ∂Σ₋F₂, ∂N₁F₂, ∂N₂F₂, ∂N₃F₂ = component(M₁₁, 2, 2), component(M₁₁, 2, 3), component(M₁₁, 2, 4), component(M₁₁, 2, 5)
    ∂Σ₊F₃, ∂N₁F₃ = component(M₁₁, 3, 1), component(M₁₁, 3, 3)
    ∂Σ₊F₄, ∂Σ₋F₄, ∂N₂F₄ = component(M₁₁, 4, 1), component(M₁₁, 4, 2), component(M₁₁, 4, 4)
    ∂Σ₊F₅, ∂Σ₋F₅, ∂N₃F₅ = component(M₁₁, 5, 1), component(M₁₁, 5, 2), component(M₁₁, 5, 5)
    M₁₂ = component(DF, 1, 2)
    DγF₁, DγF₂, DγF₃, DγF₄, DγF₅ = component(M₁₂, 1, 1), component(M₁₂, 2, 1), component(M₁₂, 3, 1), component(M₁₂, 4, 1), component(M₁₂, 5, 1)
    Dη₁F₁ = component(M₁₂, 1, 2)
    Dη₂F₄ = component(M₁₂, 4, 3)
    M₂₁ = component(DF, 2, 1)
    ∂Σ₊F₆, ∂Σ₋F₆, ∂N₁F₆, ∂N₂F₆, ∂N₃F₆ = component(M₂₁, :, 1), component(M₂₁, :, 2), component(M₂₁, :, 3), component(M₂₁, :, 4), component(M₂₁, :, 5)

    sqrt3 = sqrt(3one(real(eltype(x))))
    𝒟 = Derivative(1)

    # F₁ = 2γ * ((N₂ - N₃) ^ 2 - N₁ * (2N₁ - N₂ - N₃)) + η₁ * Σ₊ - 𝒟(Σ₊)
    lsub!(η₁*I, project!(∂Σ₊F₁, 𝒟))
    project!(∂N₁F₁, Multiplication(-2γ * (4N₁ - N₂ - N₃)))
    project!(∂N₂F₁, Multiplication(2γ * (N₁ + 2(N₂ - N₃))))
    project!(∂N₃F₁, Multiplication(2γ * (N₁ - 2(N₂ - N₃))))
    project!(Sequence(codomain(DγF₁), vec(coefficients(DγF₁))), 2((N₂ - N₃) ^ 2 - N₁ * (2N₁ - N₂ - N₃)))
    project!(Sequence(codomain(Dη₁F₁), vec(coefficients(Dη₁F₁))), Σ₊)

    # F₂ = 2γ * sqrt(3) * (N₃ - N₂) * (N₁ - N₂ - N₃) - 𝒟(Σ₋)
    project!(∂Σ₋F₂, 𝒟) .*= -1
    project!(∂N₁F₂, Multiplication(2γ*sqrt3 * (N₃ - N₂)))
    project!(∂N₂F₂, Multiplication(2γ*sqrt3 * (2N₂ - N₁)))
    project!(∂N₃F₂, Multiplication(2γ*sqrt3 * (N₁ - 2N₃)))
    project!(Sequence(codomain(DγF₂), vec(coefficients(DγF₂))), 2sqrt3 * (N₃ - N₂) * (N₁ - N₂ - N₃))

    # F₃ = 4γ * Σ₊ * N₁ - 𝒟(N₁)
    project!(∂Σ₊F₃, Multiplication(4γ * N₁))
    project!(∂N₁F₃, Multiplication(4γ * Σ₊)) .-= project(𝒟, domain(∂N₁F₃), codomain(∂N₁F₃), eltype(DF))
    project!(Sequence(codomain(DγF₃), vec(coefficients(DγF₃))), 4(Σ₊ * N₁))

    # F₄ = -2γ * N₂ * (Σ₊ + sqrt(3) * Σ₋) + η₂ - 𝒟(N₂)
    project!(∂Σ₊F₄, Multiplication(-2γ * N₂))
    project!(∂Σ₋F₄, Multiplication(-2γ*sqrt3 * N₂))
    project!(∂N₂F₄, Multiplication(-2γ * (Σ₊ + sqrt3 * Σ₋))) .-= project(𝒟, domain(∂N₂F₄), codomain(∂N₂F₄), eltype(DF))
    project!(Sequence(codomain(DγF₄), vec(coefficients(DγF₄))), -2N₂ * (Σ₊ + sqrt3 * Σ₋))
    Dη₂F₄[0,1] = one(eltype(DF))

    # F₅ = -2γ * N₃ * (Σ₊ - sqrt(3) * Σ₋) - 𝒟(N₃)
    project!(∂Σ₊F₅, Multiplication(-2γ * N₃))
    project!(∂Σ₋F₅, Multiplication(2γ*sqrt3 * N₃))
    project!(∂N₃F₅, Multiplication(-2γ * (Σ₊ - sqrt3 * Σ₋))) .-= project(𝒟, domain(∂N₃F₅), codomain(∂N₃F₅), eltype(DF))
    project!(Sequence(codomain(DγF₅), vec(coefficients(DγF₅))), -2N₃ * (Σ₊ - sqrt3 * Σ₋))

    # F₆[1] = sum(Σ₊)
    ∂Σ₊F₆[1,:] .= one(eltype(DF))
    # F₆[2] = sum(Σ₊)^2 + sum(Σ₋)^2 + sum(N₁)^2 + sum(N₂)^2 + sum(N₃)^2 - 2sum(N₁)*sum(N₂) - 2sum(N₁)*sum(N₃) - 2sum(N₂)*sum(N₃) - C₁
    ∂Σ₊F₆[2,:] .= 2s₁
    ∂Σ₋F₆[2,:] .= 2s₂
    ∂N₁F₆[2,:] .= 2(s₃ - s₄ - s₅)
    ∂N₂F₆[2,:] .= 2(s₄ - s₃ - s₅)
    ∂N₃F₆[2,:] .= 2(s₅ - s₃ - s₄)
    # F₆[3] = sum(N₁)*sum(N₂)*sum(N₃) - C₂
    ∂N₁F₆[3,:] .= s₄*s₅
    ∂N₂F₆[3,:] .= s₃*s₅
    ∂N₃F₆[3,:] .= s₃*s₄

    return DF
end

#

function opnorm_Df(x, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))

    # f₁ = 2γ * ((N₂ - N₃) ^ 2 - N₁ * (2N₁ - N₂ - N₃)) + η₁ * Σ₊
    v₁ = abs(η₁) + # ∂Σ₊f₁
        2abs(γ) * norm(4N₁ - N₂ - N₃, X_Fourier) + # ∂N₁f₁
        2abs(γ) * norm(N₁ + 2(N₂ - N₃), X_Fourier) + # ∂N₂f₁
        2abs(γ) * norm(N₁ - 2(N₂ - N₃), X_Fourier) # ∂N₃f₁

    # f₂ = 2γ * sqrt(3) * (N₃ - N₂) * (N₁ - N₂ - N₃)
    v₂ = 2abs(γ)*sqrt3 * norm(N₃ - N₂, X_Fourier) + # ∂N₁f₂
        2abs(γ)*sqrt3 * norm(2N₂ - N₁, X_Fourier) + # ∂N₂f₂
        2abs(γ)*sqrt3 * norm(N₁ - 2N₃, X_Fourier) # ∂N₃f₂

    # f₃ = 4γ * Σ₊ * N₁
    v₃ = 4abs(γ) * norm(N₁, X_Fourier) + # ∂Σ₊f₃
        4abs(γ) * norm(Σ₊, X_Fourier) # ∂N₁f₃

    # f₄ = -2γ * N₂ * (Σ₊ + sqrt(3) * Σ₋) + η₂
    v₄ = 2abs(γ) * norm(N₂, X_Fourier) + # ∂Σ₊f₄
        2abs(γ)*sqrt3 * norm(N₂, X_Fourier) + # ∂Σ₋f₄
        2abs(γ) * norm(Σ₊ + sqrt3 * Σ₋, X_Fourier) # ∂N₂f₄

    # f₅ = -2γ * N₃ * (Σ₊ - sqrt(3) * Σ₋)
    v₅ = 2abs(γ) * norm(N₃, X_Fourier) + # ∂Σ₊f₅
        2abs(γ)*sqrt3 * norm(N₃, X_Fourier) + # ∂Σ₋f₅
        2abs(γ) * norm(Σ₊ - sqrt3 * Σ₋, X_Fourier) # ∂N₃f₅

    return max(v₁, v₂, v₃, v₄, v₅)
end

function opnorm_∂Σ₊DF(x, r, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    n_N₁, n_N₂, n_N₃ = norm(N₁, X_Fourier), norm(N₂, X_Fourier), norm(N₃, X_Fourier)

    M₁₁ = 4(abs(γ) + r)
    M₁₂ = max(1, 4(n_N₁ + r), 2(n_N₂ + r), 2(n_N₃ + r))
    M₂₁ = 2

    return max(M₁₁ + M₁₂, M₂₁)
end

function opnorm_∂Σ₋DF(x, r, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))
    n_N₂, n_N₃ = norm(N₂, X_Fourier), norm(N₃, X_Fourier)

    M₁₁ = 2(abs(γ) + r)*sqrt3
    M₁₂ = max(2sqrt3 * (n_N₂ + r), 2sqrt3 * (n_N₃ + r))
    M₂₁ = 2

    return max(M₁₁ + M₁₂, M₂₁)
end

function opnorm_∂N₁DF(x, r, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))
    n_Σ₊ = norm(Σ₊, X_Fourier)
    n_N₁, n_N₂, n_N₃ = norm(N₁, X_Fourier), norm(N₂, X_Fourier), norm(N₃, X_Fourier)

    M₁₁ = 12(abs(γ) + r)
    M₁₂ = max(2(4n_N₁ + n_N₂ + n_N₃ + 6r), 2sqrt3 * (n_N₂ + n_N₃ + 2r), 4(n_Σ₊ + r))
    M₂₁ = max(6, n_N₂ + n_N₃ + 2r)

    return max(M₁₁ + M₁₂, M₂₁)
end

function opnorm_∂N₂DF(x, r, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))
    n_Σ₊, n_Σ₋ = norm(Σ₊, X_Fourier), norm(Σ₋, X_Fourier)
    n_N₁, n_N₂, n_N₃ = norm(N₁, X_Fourier), norm(N₂, X_Fourier), norm(N₃, X_Fourier)

    M₁₁ = 6(abs(γ) + r)*sqrt3
    M₁₂ = max(2(n_N₁ + 2n_N₂ + 2n_N₃ + 5r), 2sqrt3 * (n_N₁ + 2n_N₂ + 3r), 2(n_Σ₊ + r + sqrt3 * (n_Σ₋ + r)))
    M₂₁ = max(6, n_N₁ + n_N₃ + 2r)

    return max(M₁₁ + M₁₂, M₂₁)
end

function opnorm_∂N₃DF(x, r, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))
    n_Σ₊, n_Σ₋ = norm(Σ₊, X_Fourier), norm(Σ₋, X_Fourier)
    n_N₁, n_N₂, n_N₃ = norm(N₁, X_Fourier), norm(N₂, X_Fourier), norm(N₃, X_Fourier)

    M₁₁ = 6(abs(γ) + r)*sqrt3
    M₁₂ = max(2(n_N₁ + 2n_N₂ + 2n_N₃ + 5r), 2sqrt3 * (n_N₁ + 2n_N₃ + 3r), 2(n_Σ₊ + r + sqrt3 * (n_Σ₋ + r)))
    M₂₁ = max(6, n_N₁ + n_N₂ + 2r)

    return max(M₁₁ + M₁₂, M₂₁)
end

function opnorm_∂γDF(x, r, X_Fourier)
    Σ₊, Σ₋, N₁, N₂, N₃ = eachcomponent(component(x, 1))
    γ, η₁, η₂ = component(x, 2)

    sqrt3 = sqrt(3one(real(eltype(x))))
    n_Σ₊, n_Σ₋ = norm(Σ₊, X_Fourier), norm(Σ₋, X_Fourier)
    n_N₁, n_N₂, n_N₃ = norm(N₁, X_Fourier), norm(N₂, X_Fourier), norm(N₃, X_Fourier)

    return max(
        2(4n_N₁ + n_N₂ + n_N₃ + 6r) + 4(n_N₁ + 2n_N₂ + 2n_N₃+ 5r),
        2sqrt3 * (2n_N₁ + 3n_N₂ + 3n_N₃ + 8r),
        4(n_Σ₊ + n_N₁ + 2r),
        2(n_N₂ + r) + 2sqrt3 * (n_N₂ + r) + 2(n_Σ₊ + r + sqrt3 * (n_Σ₋ + r)),
        2(n_N₃ + r) + 2sqrt3 * (n_N₃ + r) + 2(n_Σ₊ + r + sqrt3 * (n_Σ₋ + r))
        )
end

#

function check_not_identically_zero(f, r)
    tmp = f(Interval(0.0))
    for tⱼ ∈ 0.0:1e-3:2π
        evaluate!(tmp, f, Interval(tⱼ))
        real_f = real(tmp[0])
        if (inf(real_f - r) > 0) || (sup(real_f + r) < 0)
            return true
        end
    end
    return false
end

function check_no_zeros(f, bound_Df, r, name, checknumber)
    π2 = sup(2convert(Interval{Float64}, π))
    tⱼ = Interval(0.0)
    tmp = f(tⱼ)
    count = 0
    while tⱼ < π2
        evaluate!(tmp, f, tⱼ)
        real_f = real(tmp[0])
        if (inf(real_f - r) ≤ 0) && (0 ≤ sup(real_f + r))
            return false
        else
            tⱼ = Interval(inf(tⱼ + prevfloat(inf((abs(real_f) - r)/bound_Df))))
            count += 1
        end
        if count % 10_000 == 1
            print("        • check ", checknumber, ": ", name, " has no zeros in [0, ", round(inf(tⱼ); digits = 6), ")    \u001b[1000D")
        end
    end
    return true
end


text:
using VT100

@nodegen mutable struct VT100Node <: AbstractNode
  screen::ScreenEmulator
end
VT100Node(w, h) = VT100Node(ScreenEmulator(w, h))


text:
mutable struct Simulation{S <: DiseaseStateSequence, M <: ILM}
  time::Float64
  iterations::Int64
  population::Population
  risk_functions::RiskFunctions{S}
  risk_parameters::RiskParameters{S}
  substitution_model::Union{Nothing, NucleicAcidSubstitutionModel}
  disease_states::DiseaseStates
  transmission_rates::TransmissionRates
  event_rates::EventRates{S}
  events::Events{S}
  transmission_network::TransmissionNetwork

  function Simulation(
    pop::Population,
    states::DiseaseStates,
    rf::RiskFunctions{S},
    rp::RiskParameters{S};
    time::Float64 = 0.0,
    skip_checks::Bool=false) where {
      S <: DiseaseStateSequence}
    @debug "Initializing $S TNILM Simulation with the following starting states:" states = convert(DiseaseStates, S) counts = [sum(states .== Ref(s)) for s in convert(DiseaseStates, S)]
    if !skip_checks
      if length(states) != individuals(pop)
        error("Length of initial disease state vector must match number of individuals")
      elseif !all(in.(states, Ref(convert(DiseaseStates, S))))
        error("All states in initial disease state vector must be valid within specified epidemic model")
      end
    end
    tr = initialize(TransmissionRates, states, pop, rf, rp)
    rates = initialize(EventRates, tr, states, pop, rf, rp)
    events = Events{S}(states)
    net = TransmissionNetwork(states)
    return new{S, TNILM}(time, 0, pop, rf, rp, nothing, copy(states), tr, rates, events, net)
  end

  function Simulation(
    pop::Population,
    states::DiseaseStates,
    rf::RiskFunctions{S},
    rp::RiskParameters{S},
    sm::NASM;
    time::Float64 = 0.0,
    skip_checks::Bool=false) where {
      S <: DiseaseStateSequence}
    @debug "Initializing $S PhyloILM Simulation with the following starting states:" states
    if !skip_checks
      if length(states) != individuals(pop)
        error("Length of initial disease state vector must match number of individuals")
      elseif !all(in.(states, Ref(convert(DiseaseStates, S))))
        error("All states in initial disease state vector must be valid within specified epidemic model")
      end
    end
    tr = initialize(TransmissionRates, states, pop, rf, rp)
    rates = initialize(EventRates, tr, states, pop, rf, rp)
    events = Events{S}(states)
    net = TransmissionNetwork(states)
    return new{S, PhyloILM}(time, 0, pop, rf, rp, sm, copy(states), tr, rates, events, net)
  end
end

function Simulation(
  pop::Population,
  rf::RiskFunctions{S},
  rp::RiskParameters{S},
  sm::NASM;
  time::Float64 = 0.0) where {
    S <: DiseaseStateSequence}
  return Simulation(pop, fill(State_S, individuals(pop)), rf, rp, sm, time=time, skip_checks=true)
end

function Simulation(
  pop::Population,
  rf::RiskFunctions{S},
  rp::RiskParameters{S};
  time::Float64 = 0.0) where {
    S <: DiseaseStateSequence}
  return Simulation(pop, fill(State_S, individuals(pop)), rf, rp, time=time, skip_checks=true)
end

function Base.show(
  io::IO,
  x::Simulation{S, M}) where {
    S <: DiseaseStateSequence,
    M <: ILM}
  y = "$S $M epidemic simulation @ time = $(round(x.time, digits=2))\n"
  for s in convert(DiseaseStates, S)
    y *= "\n$s = $(sum(x.disease_states .== Ref(s)))"
  end
  return print(io, y)
end