text:
struct Giωn
    ωn::Real
    val::Complex
end

struct MasubaraGF
    number::Int64
    GF::Vector{Giωn}
end

function make_input(path::String, num::Int64)
    data = readdlm(path)
    fnum = min(num, length(data[:,1])) # number of frenquencies used
    res = [Giωn(data[i,1],data[i,2]+1.0im*data[i,3]) for i = 1:fnum]
    return MasubaraGF(fnum, res |> reverse)
end

function make_input(path::String)
    data = readdlm(path)
    fnum = length(data[:,1]) # number of frenquencies used
    res = [Giωn(data[i,1],data[i,2]+1.0im*data[i,3]) for i = 1:fnum]
    return MasubaraGF(fnum, res |> reverse)
end

function pick_matrix(G::MasubaraGF)
    # input :freq: Masubara frenquencies; val: G
    N = G.number; g = G.GF
    pick = zeros(ComplexF64,N, N)
    for i = 1:N, j = 1:N
        up = 1 - λ(g[i])*λ(g[j])'
        down = 1 - h(1.0im*g[i].ωn)*h(1.0im*g[j].ωn)'
        pick[i,j] = up/down
    end
    return pick
end

function isNevanlinnasolvable(G::MasubaraGF; err::Float64 = 1.e-6)
    # check if the pick_matrix is semi positive semidefinite
    pick = pick_matrix(G)
    evals = eigvals(pick)
    return all((evals .+ err) .>= 0), minimum(evals)
end

function λ(G::Giωn)
    return h(-G.val)
end

function θk(fac::Matrix, θp::Number)
    # θp = θ(k-1)
    up = -fac[2,2] * θp + fac[1,2]
    down = fac[2,1] * θp - fac[1,1]
    return up/down
end

function core(G::MasubaraGF; etype::DataType=Float64)
    dtype = Complex{etype}
    M = G.number; g = G.GF
    phis = zeros(dtype, M) #store λk(iω_(k+1))
    phis[1] = λ(g[1])
    abcds = Vector{Matrix{dtype}}(undef,M)
    for i = 1:M abcds[i] = eye(dtype,2) end
    for j = 1:(M-1)
        for k=j:M
            prod = zeros(dtype,2,2)
            prod[1,1] = h1(1.0im*g[k].ωn,1.0im*g[j].ωn)
            prod[1,2] = phis[j]
            prod[2,1] = phis[j]' * h1(1.0im*g[k].ωn,1.0im*g[j].ωn)
            prod[2,2] = 1. + 0.0im
            abcds[k] *= prod
        end
        phis[j+1] = θk(abcds[j+1], λ(g[j+1]))
    end
    return phis
end

function evaluation(z::Number, G::MasubaraGF;
    optim::Symbol = :none, etype::DataType=Float64)
    dtype = Complex{etype}
    M = G.number
    g = G.GF
    phis = core(G, etype=etype)
    abcd = eye(dtype, 2)
    for j = 1:M
        prod = zeros(dtype, 2,2)
        prod[1,1] = h1(z,1.0im*g[j].ωn)
        prod[1,2]= phis[j]
        prod[2,1] = phis[j]' * h1(z,1.0im*g[j].ωn)
        prod[2,2] = 1. + 0.0im
        abcd *= prod
    end

    if optim == :none
        θm(z::Number) = 0.0 + 0.0im
    end

    up = abcd[1,1] * θm(z) + abcd[1,2]
    down = abcd[2,1] * θm(z) + abcd[2,2]
    θ = up/down
    return invh(θ)
end

function spectrum_density(ω::Real, η::Real, G::MasubaraGF;
    etype::DataType=Float64)
    z = ω + 1.0im * η
    return 2*(evaluation(z, G, etype=etype) |> imag)
end


text:
export stnKD, PointNet

stnKD(K::Int) = Chain(
    Conv((1,), K=>64, relu),
    BatchNorm(64),
    Conv((1,), 64=>128, relu),
    BatchNorm(128),
    Conv((1,), 128=>1024, relu),
    BatchNorm(1024),
    x -> maximum(x, dims=1),
    x -> reshape(x, :, size(x, 3)),
    Dense(1024, 512, relu),
#     BatchNorm(512),
    Dense(512, 256, relu),
    BatchNorm(256),
    Dense(256,K*K),
    x -> reshape(x, K, K, size(x, 2)),
    x -> x .+ reshape(Matrix{Float64}(I, K, K), K, K, 1),
    x -> PermutedDimsArray(x, (2,1,3)),
)

# PointNet(k::Int) = Chain(
#     x -> batched_mul(x, stn(x)),
#     Conv((1,), 3=>64),
#     BatchNorm(64),
#     x -> batched_mul(x, fstn(x)),
#     Conv((1,), 64=>128, relu),
#     BatchNorm(128),
#     Conv((1,), 128=>1024),
#     BatchNorm(1024),
#     x -> maximum(x, dims=1),
#     x -> reshape(x, 1024, :),
#     Dense(1024, 512, relu),
#     BatchNorm(512),
#     Dense(512, 256, relu),
#     Dropout(0.4),
#     BatchNorm(256),
#     Dense(256, k, relu),
#     x -> softmax(x, dims=1),
# )

struct PointNet
    stn
    fstn
    conv_block1
    feat
    cls
end

function PointNet(num_classes::Int=10, K::Int=64)
    stn = stnKD(3)
    fstn = stnKD(K)
    conv_block1 = conv_bn_block(3, 64, (1,))
    feat = Chain(
        Conv((1,), 64=>128, relu),
        BatchNorm(128),
        Conv((1,), 128=>1024),
        BatchNorm(1024),
        x -> maximum(x, dims=1),
        x -> reshape(x, 1024, :),
        Dense(1024, 512, relu),
        BatchNorm(512),
        Dense(512, 256, relu),
        Dropout(0.4),
        BatchNorm(256),
        )
    cls = Dense(256, num_classes, relu)
    PointNet(stn,fstn, conv_block1, feat, cls)
end

function (m::PointNet)(X)

        # X: [N, 3, B]

        X = batched_mul(X, m.stn(X))
        # X: [3, 3, B]

        X = m.conv_block1(X)
        # X: [3, 64, B]

        X = batched_mul(X, m.fstn(X))
        # X: [3, 64, B]

        X = m.feat(X)
        # X: [256, B]

        X = m.cls(X)
        # X: [num_classes, B]

        return softmax(X, dims=1)
end

@functor PointNet

ext: jl


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

alpha_init = 10. *pi/180
alphadot_init = 0.
h_init = 0.
hdot_init = 0.
u = 0.467
udot = 0
kinem = UNSflow.KinemPar2DOF(alpha_init, h_init, alphadot_init, hdot_init, u)

x_alpha = 0.05
r_alpha = 0.5
kappa = 0.05
w_alpha = 1.
w_h = 1.
w_alphadot = 0.
w_hdot = 0.
cubic_h_1 = 1.
cubic_h_3 = 0.
cubic_alpha_1 = 1.
cubic_alpha_3 = 0.
strpar = UNSflow.TwoDOFPar(x_alpha, r_alpha, kappa, w_alpha, w_h, w_alphadot, w_hdot, cubic_h_1, cubic_h_3, cubic_alpha_1, cubic_alpha_3)

#Dummy kinematic definitions to initialise surface
alphadef = UNSflow.ConstDef(alpha_init)
hdef = UNSflow.ConstDef(h_init)
udef = UNSflow.ConstDef(1.) #This is relative to uref
startkinem = UNSflow.KinemDef(alphadef, hdef, udef)

lespcrit = [0.11;]

pvt = 0.35

c = 1.

surf = UNSflow.TwoDSurf("FlatPlate", pvt, startkinem, lespcrit, c=c, uref=u)

curfield = UNSflow.TwoDFlowField()

dtstar = 0.015

nsteps = 50000

t_tot = nsteps * dtstar / u

startflag = 0

writeflag = 0
writeInterval = dtstar * nsteps/20.

writeInterval = t_tot/20.

delvort = UNSflow.delSpalart(500, 12, 1e-5)

mat, surf, curfield = UNSflow.ldvm2DOF(surf, curfield, strpar, kinem, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

UNSflow.makeForcePlots2D()

#cleanWrite()


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

struct NDArrayStyle{N} <: Broadcast.AbstractArrayStyle{N} end
NDArrayStyle(::Val{N}) where N        = NDArrayStyle{N}()
NDArrayStyle{M}(::Val{N}) where {N,M} = NDArrayStyle{N}()

# Determin the output type
Base.BroadcastStyle(::Type{<:NDArray{T,N}}) where {T,N} = NDArrayStyle{N}()

Base.broadcastable(x::NDArray) = x

# Make it non-lazy
broadcasted(f, x::NDArray, args...)    = f(x, args...)
broadcasted(f, y, x::NDArray, args...) = f(y, x, args...)
broadcasted(f, x::NDArray{T,N}, y::NDArray{T,N}, args...) where {T,N} =
  f(x, y, args...)


text:
@testset "nash_conv" begin
    env = OpenSpielEnv(
        "kuhn_poker";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = TabularRandomPolicy()

    @test RLZoo.nash_conv(p, env) == 11 / 12

    p = get_optimal_kuhn_policy()
    @test RLZoo.nash_conv(p, env) == 0.0

    env = OpenSpielEnv(
        "leduc_poker";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 4.747222222222222

    env = OpenSpielEnv(
        "kuhn_poker(players=3)";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 2.0624999999999996

    env = OpenSpielEnv(
        "kuhn_poker(players=4)";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = TabularRandomPolicy()
    @test RLZoo.nash_conv(p, env) ≈ 3.4760416666666663
end


text:
using Accessors

"""
    keysort(nt::NamedTuple)

Recursively sort the keys of a NamedTuple.

EXAMPLE:
```
julia> x = randnt(3,2)
(p = (z = :z, e = :e, s = :s), w = (g = :g, o = :o), g = (c = :c, f = :f, k = :k))

julia> keysort(x)
(g = (c = :c, f = :f, k = :k), p = (e = :e, s = :s, z = :z), w = (g = :g, o = :o))
```
"""
function keysort end

_keysort(nt::NamedTuple{(), Tuple{}}) = nt

function _keysort(nt::NamedTuple{K,V}) where {K,V}
    π = sortperm(collect(K))
    k = K[π]
    v = @inbounds (_keysort.(values(nt)))[π]
    return namedtuple(k)(v)
end

_keysort(t::Tuple) = _keysort.(t)
_keysort(x) = x

export keysort

struct Lenses{T}
    ls::T
end

export call
call(f,args...; kwargs...) = f(args...; kwargs...)

(ℓ::Lenses)(nt::NamedTuple) = Base.Fix2(call, nt).(ℓ.ls)

@generated function keysort(nt::NamedTuple{K,V}) where {K,V}
    s = _keysort(schema(nt))
    ℓ = Lenses(lenses(s))
    return :(leaf_setter($s)($ℓ(nt)...))
end

keysort(t::T) where {T<:Tuple} = keysort.(t)
keysort(nt::NamedTuple{(), Tuple{}}) = nt
keysort(x) = x

function keysort(lm::LazyMerge)
    x = getfield(lm, :x)
    y = getfield(lm, :y)
    return lazymerge(keysort(x), keysort(y))
end


text:
include("BEMModel.jl")
using Base.Test

################################################################################
# Connectivity dictionary key functions
#
"""
    inverse_key(key::AbstractArray)

Given the key for the mesh connectivity d -> d', return the key for the
connectivity d' -> d
"""
function inverse_key(key::AbstractString)
    if key == "V2V"
        return "V2V"
    elseif key == "V2E"
        return "E2V"
    elseif key == "V2F"
        return "F2V"
    elseif key == "E2V"
        return "V2E"
    elseif key == "E2E"
        return "E2E"
    elseif key == "E2F"
        return "F2E"
    elseif key == "F2V"
        return "V2F"
    elseif key == "F2E"
        return "E2F"
    elseif key == "F2F"
        return "F2F"
    elseif key == "E2N"
        return "N2E"
    elseif key == "N2E"
        return "E2N"
    else
        throw(KeyError(key))
    end
end
################################################################################

################################################################################
# Connectivity access functions
#
"""
    get_indices(connectivity::MeshConnectivity, i::Int64)

Get the indices of the desired mesh entity. The argument `connectivity`
specifies a mesh connectivity d -> d'. The argument `i` specifies the mesh
entitiy (d, i). This function returns an array of `indices` that specifies
which entities are connected to (d, i). That is, it allows one to determine the
set of entities connected to (d, i) defined by { (d',j) : j ∈ `indices` }.
"""
function get_indices(connectivity::MeshConnectivity, i::Int64)
    start_index = connectivity.offsets[i]
    end_index = connectivity.offsets[i + 1] - 1
    indices = connectivity.indices[start_index:end_index]
    return indices
end

"""
    get_vertices(connectivity::MeshConnectivity, geometry::MeshGeometry,
                 i::Int64)

Get the vertices connected to the desired mesh entity.
"""
function get_vertices(connectivity::MeshConnectivity, geometry::MeshGeometry,
                      i::Int64)
    indices = get_indices(connectivity, i)
    vertices = geometry.vertices[indices]
    return vertices
end

"""
    update_connectivity!(connectivity::MeshConnectivity, index::Int64,
                         rep_number::Int64)

Update the mesh connectivity fields between specified mesh entities. If the mesh
connectivity we wish to update is d -> d', the `index` argument specifies the
mesh entitiy (d, index) while the argument `rep_number` specifies the mesh
entitiy (d', rep_number).
"""
function update_connectivity!(connectivity::MeshConnectivity, index::Int64,
                              rep_number::Int64)
    # The index in the indices array to add the replacement number
    rep_index = connectivity.offsets[index]
    if rep_index > length(connectivity.indices)
        # Add the replacement number to the end
        push!(connectivity.indices, rep_number)
    else
        # Insert the replacement number at replacement index. This looks strange
        # but this is how it's done in Julia...
        splice!(connectivity.indices, rep_index, [rep_number, connectivity.indices[rep_index]])
    end
    # Offset all the entity offsets after index to take into account the fact
    # we just added an index.
    for i in index + 1:length(connectivity.offsets)
        connectivity.offsets[i] += 1
    end
end

################################################################################
# Computational mesh functions and algorithms based on A. Logg paper
#
"""
    connectivity_transpose(connectivity::MeshConnectivity)

Given a mesh connectivity d -> d', compute the connectivity d' -> d for d' < d.
"""
function connectivity_transpose!(topology::MeshTopology, key::AbstractString)
    # Get the transpose key
    transpose_key = inverse_key(key)

    # Starting Connectivity
    connectivity = topology.connectivity[key]

    # Number of dimension d' entities
    ndprime = maximum(connectivity.indices)
    # Number of dimension d entities
    nd = length(connectivity.offsets) - 1

    # Initialize connectivity indices
    new_indices = Int64[] # A priori unknown index array length

    # Initialize the offset values with 1. This will be used to place the
    # cell indices in the correct position of the array. The algorithm
    # increments the appropriate offset elements after placing an index
    new_offsets = ones(Int64, ndprime + 1) # Known number of d' entities

    topology.connectivity[transpose_key] = MeshConnectivity(new_indices,
                                                            new_offsets)
    # The idea here is to use the indices of d' connected to the jth
    # d-entity as an index into the current offset array. This value in the
    # current offset array gives the location to place the d-entity number in
    # the indices array. After placing the d-entity index in the indices array,
    # increment the elements in the offsets array at positions greater than the
    # d'-entity number by 1. After looping through all the cells, this will give
    # the correct connectivity information.
    for j in 1:nd
        # Get the d' indices for the jth cell
        indices = get_indices(connectivity, j)
        for index in indices
            update_connectivity!(topology.connectivity[transpose_key], index, j)
        end
    end
end

"""
    connectivity_build!(topology::MeshTopology)
"""
# TO DO: Implament this from the A. Logg. Paper. This may me needed for 3D
# Mesh computation
function connectivity_build!(topology::MeshTopology)
end

"""
    connectivity_intersection(topology::MeshTopology, key_start::AbstractString,
                              key_end::AbstractString)

Compute the connectivity d -> d' from d -> d" and d" -> d' for d >= d'.
key_start is the key for d -> d" and key_end is the key for d" -> d'.
"""
# TO DO: Implament this from the A. Logg. Paper. This may me needed for 3D
# Mesh computation
function connectivity_intersection!(topology::MeshTopology,
                                    key_start::AbstractString,
                                    key_end::AbstractString)
end

"""
    generate_connectivity!(mesh::Mesh, key::AbstractString)

Build the mesh connectivity corresponding to the key parameter.
"""
# TO DO: Might have to make this general based on the A. Logg. algorithm.
# For now, only special cases for key are implemented. Currently only
# to find E2V for 2D meshes.
function generate_connectivity!(mesh::Mesh{Dimension2}, key::AbstractString)
end
################################################################################

################################################################################
# Miscellaneous mesh functions
#
"""
    generate_normals!(mesh::Mesh)

Compute and store the unit normal vectors of the mesh. Each mesh cell has an
associated unit normal. This normal is used during integral evaluations of the
normal derivative of the free space Green's function as well as in determining
if discontinuous elements are required for flux on the boundary.

Currently only implemented for 2D linear elements.
"""
# TO DO: It would be nice if we didn't have to count on the
# vertices being given in order. This is minor, so fix up other things first.
function generate_normals!(mesh::Mesh{Dimension2})
    # Number of mesh cells
    nd = length(mesh.topology.connectivity["E2V"].offsets) - 1

    # Treat each mesh segment as a directed line segment. So each segment
    # can be written as [dx, dy] = [pt2.x - pt1.x, pt2.y - pt1.y]. The
    # outward pointing normal is given by [dy, -dx]. This requires that the
    # vertices of the 2D mesh are given in counterclockwise order.
    for i in 1:nd
        vertices = get_vertices(mesh.topology.connectivity["E2V"],
                                mesh.geometry, i)
        dx = vertices[2][1] - vertices[1][1]
        dy = vertices[2][2] - vertices[1][2]
        push!(mesh.geometry.normals, Point(dy, -1.0*dx))
    end
end

# TO DO: 3D version. Some linear algebra stuff to look up (cross products, etc.)
function generate_normals!(mesh::Mesh{Dimension3})
end
################################################################################

################################################################################
# Polynomial interpolation node functions
#
"""
    node_stepping(vertices::Array{Point}, include_verts::Array{Bool},
                  order:Int64)

Compute shape function interpolation nodes on boundary edge.
# Arguments
* `vertices`: the vertices defining the endpoints of the edge
* `include_verts`: an array of booleans indicating if a vertex
should be included as a node. include_verts[i] = true means that vertices[i]
will be included as a vertex
* `order`: the order of the polynomial on the edge. The convention is
different from mathematical convention: 1 = constant, 2 = linear, 3 = quadratic,
etc. This means the number of nodes is equal to order.
"""
function node_stepping(vertices::Array{Point}, include_verts::Array{Bool},
                       order::Int64)
    @test length(vertices) == 2 # The number of vertices for an edge is always 2
    @test length(include_verts) == 2 # Must match the number of vertices.

    nodes = Array{Point}(order)
    delta = Point(vertices[2] - vertices[1])/(order - 1) # step size
    for i in 1:order
        nodes[i] = vertices[1] + delta*(i - 1)
    end

    if include_verts[1] == false # remove the first node
        shift!(nodes)
    end
    if include_verts[2] == false # remove the last node
        pop!(nodes)
    end
    return nodes
end

"""
    function make_edge_nodes!(mesh::Mesh, index::Int64, order::Int64)

Build the node geometry and E2N connectivity for edge (2, `index`). The number
of nodes added is determined by the order of the polynomial on the edge as
indicated by `order`.
"""
##### TO DO: Can this be tidied up. A lots going on here with some repetition.
#####        Can there be some more functions?
function make_edge_nodes!(mesh::Mesh, index::Int64, order::Int64)
    # Get the vertices of the current edge
    vertices = get_vertices(mesh.topology.connectivity["E2V"], mesh.geometry,
                            index)
    # Get the indices of these vertices
    vertex_indices = get_indices(mesh.topology.connectivity["E2V"], index)
    @test length(vertices) == 2 # Edges should always have 2 vertices
    @test length(vertex_indices) == 2

    if order == 1 # Special discontinuous case for constant elements
        push!(mesh.geometry.nodes, mean(vertices)) # Node in the centre of edge
        global_number = length(mesh.geometry.nodes) # Total number of mesh nodes
        # Edge 'index' is connected to the node at index 'global_nodes'
        push!(mesh.topology.connectivity["E2N"].indices, global_number)
        # Increment all other edge offset indices by 1
        mesh.topology.connectivity["E2N"].offsets[index + 1:end] += 1
    else
        # Current number of nodes on edge
        local_nodes = mesh.topology.connectivity["E2N"].offsets[index + 1] -
                      mesh.topology.connectivity["E2N"].offsets[index]
        if local_nodes == 0
            # Make all nodes including nodes coinciding with vertices
            include_verts = [true, true] # Include both vertices as nodes
            nodes = node_stepping(vertices, include_verts, order)
            for i in 1:length(nodes)
                # Add nodes for the current edge (given by index)
                push!(mesh.geometry.nodes, nodes[i]) # Add node to geometry
                # Index of recently added node
                global_number = length(mesh.geometry.nodes)
                update_connectivity!(mesh.topology.connectivity["E2N"], index,
                                     global_number)
                # Since local_nodes == 0, the edges connected to vertices[1] and
                # vertices[2] share nodes[1] and nodes[end] respectively. So
                # we add more copies of global_number to the indices array
                # for these connected edges and modify offsets appropriately
                if i == 1
                    edge_indices = get_indices(mesh.topology.connectivity["V2E"],
                                               vertex_indices[1])
                    for edge in edge_indices
                        if edge == index
                            # This is the current edge we're operating on
                            # so we skip over this. We're only concerned with
                            # the connected edges here.
                        else
                            update_connectivity!(mesh.topology.connectivity["E2N"],
                                                 edge, global_number)
                        end
                    end
                end
                if i == length(nodes)
                    edge_indices = get_indices(mesh.topology.connectivity["V2E"],
                                               vertex_indices[end])
                    for edge in edge_indices
                        if edge == index
                            # This is the current edge we're operating on
                            # so we skip over this. We're only concerned with
                            # the connected edges here.
                        else
                            update_connectivity!(mesh.topology.connectivity["E2N"],
                                                 edge, global_number)
                        end
                    end
                end
            end
        elseif local_nodes == 1 # One of the edge vertices already has a node.
            # To determine this vertex, we first find the edges connected to
            # each vertex.
            i = 1
            include_verts = [true, true] # initialize flags
            for vertex_index in vertex_indices
                edge_indices = get_indices(mesh.topology.connectivity["V2E"],
                                           vertex_index)
                # Since local_nodes == 1, at least one of the edges connected
                # to one of the vertices must already have nodes. At the same
                # time, it MUST be the case that the other vertex CANNOT be
                # connected to edges with nodes (otherwise local_nodes == 2).
                # So the strategy here is to inspect the edges connected to each
                # vertex to find which of the vertices is connected to an edge
                # with nodes.
                for edge in edge_indices
                    if edge == index
                        # This is the current edge we're operating on
                        # so we skip over this. We're only concerned with
                        # the connected edges here.
                    else
                        edge_nodes = mesh.topology.connectivity["E2N"].offsets[edge + 1] -
                                     mesh.topology.connectivity["E2N"].offsets[edge]
                        if edge_nodes == order # This means the edge has nodes
                            include_verts[i] = false
                        end
                    end
                end
                i += 1
            end
            # As stated earlier, only ONE of the vertices can have a connected
            # edge with nodes. We test that here.
            @test (include_verts[1] == true) $ (include_verts[2] == true)

            nodes = node_stepping(vertices, include_verts, order)
            for i in 1:length(nodes)
                # Add nodes for the current edge
                push!(mesh.geometry.nodes, nodes[i]) # Add node to geometry
                # Index of recently added node
                global_number = length(mesh.geometry.nodes)
                update_connectivity!(mesh.topology.connectivity["E2N"], index,
                                     global_number)

                # Update the connected edges to the vertex without the initial node
                if (i == 1) && (include_verts[1] == true)
                    edge_indices = get_indices(mesh.topology.connectivity["V2E"],
                                               vertex_indices[1])
                    for edge in edge_indices
                        if edge == index
                            # This is the current edge we're operating on
                            # so we skip over this. We're only concerned with
                            # the connected edges here.
                        else
                            update_connectivity!(mesh.topology.connectivity["E2N"],
                                                 edge, global_number)
                        end
                    end
                end
                if (i == length(nodes)) && (include_verts[2] == true)
                    edge_indices = get_indices(mesh.topology.connectivity["V2E"],
                                               vertex_indices[end])
                    for edge in edge_indices
                        if edge == index
                        else
                            update_connectivity!(mesh.topology.connectivity["E2N"],
                                                 edge, global_number)

                        end
                    end
                end
            end
        elseif local_nodes == 2
            # Only make the internal nodes. No need to worry about updating
            # connected edges in this case
            include_verts = [false, false]
            nodes = node_stepping(vertices, include_verts, order)
            for i in 1:length(nodes)
                push!(mesh.geometry.nodes, nodes[i])
                global_number = length(mesh.geometry.nodes)
                update_connectivity!(mesh.topology.connectivity["E2N"], index,
                                     global_number)
            end
        else
            throw(ErrorException("Too many nodes on edge"))
        end
    end
end

"""
    generate_nodes!(mesh::Mesh, order::Int64)

Add the polynomial interpolation nodes to the mesh. This function is called once
during the matrix assembly procedure to calculate the geometry and topology of
the nodes with respect to the input mesh.
"""
function generate_nodes!(mesh::Mesh, order::Int64)
    # For this BEM implementation, the interpolation nodes will always be on
    # the mesh edges, for both 2D and 3D meshes.

    # Check if the edge topology is known:
    if haskey(mesh.topology.connectivity, "E2V") == false # Edges unknown
        # Generate the edge -> vertex connectivity. Note that we start with
        # E2V for 2d meshes (since the edge is of maximal dimension D for
        # 2d meshes). So this will only be called for 3d meshes.
        # TO DO: generate_connectivity is not implemented yet.
        generate_connectivity!(mesh, "E2V") # Does not exist yet
    end
    # We need to know which edges are connected to a given vertex. So we need
    # to find the vertex -> edge connectivity. This is simply the transpose of
    # E2V
    if haskey(mesh.topology.connectivity, "V2E") == false
        connectivity_transpose!(mesh.topology, "E2V")
    end

    nedges = length(mesh.topology.connectivity["E2V"].offsets) - 1 # Number of edges

    # The nodes have both a topological and geometric aspect. The geometric
    # aspect is simply their cartesian coordinates, similar to the mesh
    # vertices. Their topological aspect is how they are connected to the mesh.
    # In particular, we are interested in how they are connected to the edges
    # and cells of the mesh.

    # Since the number of nodes for edge k = offsets[k+1] - offsets[k], we give
    # a starting representation of zero nodes for each edge
    offsets = ones(Int64, nedges + 1) # Offsets for edge to node connectivity
    indices = Int64[] # Indices for edge to node connectivity
    mesh.topology.connectivity["E2N"] = MeshConnectivity(indices, offsets)

    # Each edge must have a number of nodes equal to order. Note that edges
    # will share nodes whenever they meet at a vertex.
    for i in 1:nedges
        # Build nodes for each edge
        make_edge_nodes!(mesh, i, order)
    end
    ###### TO DO: This only works for dimension 2!!!! Dimension 3 requires a
    ######        general connectivity generator. We can use transpose here
    ######        because the edges are the cells in dimension 2. So to get
    ######        the node to cell connectivity, we simply take the transpose of
    ######        the E2N
    # Make the node -> cell connectivity

    # TO DO: calculating connectivity transpose after making the nodes takes
    # some time. There's a potential performance increase by calculating the
    # node to cell connectivity while the nodes are being generated in
    # make_edge_nodes!. Performance is already fairly good, so this is a minor
    # issue.
    connectivity_transpose!(mesh.topology, "E2N") # Only for dim 2
end
################################################################################
"""
    discretize_mesh!(mesh::Mesh, order::Int64)

Compute the Lagrange interpolation node points (and the associated connectivity)
and the mesh cell outward unit normals
"""
function discretize_mesh!(mesh::Mesh, order::Int64)
    # Assemble global node system. Skip if the nodes are already in place.
    if length(mesh.geometry.nodes) == 0
        generate_nodes!(mesh, order)
    end
    # Compute mesh normals
    if length(mesh.geometry.normals) == 0
        generate_normals!(mesh)
    end
end

"""
    initialize_matrix(mesh::Mesh, mesh_bcs::MeshBoundaryCondition)

Initialize the system matrix. Returns an NxN matrix with dummy values where N
equals the number of nodes plus the number of split collocation points.
"""
function initialize_matrix(mesh::Mesh{Dimension2},
                           mesh_bcs::MeshBoundaryCondition)
    # The minimum size of the matrix is equal to the number of nodes in the mesh
    N = length(mesh.geometry.nodes)

    # We add a collocation equation whenever there are two degrees of freedom at
    # a mesh node. This happens at corner vertices whose connected edges both
    # have a Dirichlet boundary condition
    index = 1
    for vertex in mesh.geometry.vertices
        edge_indices = get_indices(mesh.topology.connectivity["V2E"], index)
        @test length(edge_indices) == 2 # Must be true for 2D mesh
        if isa(mesh_bcs.bcs[mesh_bcs.markers[edge_indices[1]]], Dirichlet) &
           isa(mesh_bcs.bcs[mesh_bcs.markers[edge_indices[2]]], Dirichlet)
           # Both boundary connected edges have a Dirichlet condition
           N += 1
       end
       index += 1
   end
   A = Array{Float64}(N,N) # Initialize the matrix
   return A
end

"""
    collocate!(A::Array{Float64}, index::Int64, coll_point::Point,
               geometry::MeshGeometry, bcs::MeshBoundaryCondition)

Modify row `index` of matrix `A` with a collocation equation corresponding to
the impulse source located at `coll_point`.
"""
function collocate!(A::Array{Float64}, index::Int64, coll_point::Point,
                    geometry::MeshGeometry, bcs::MeshBoundaryCondition)
end

"""
    assemble_system(mesh::Mesh, mesh_bcs::MeshBoundaryCondition)

Assemble the BEM system Ax = b. This is based on the direct BEM where the system
of equations is generated by selecting collocation points at every element
interpolation node. Extra degrees of freedom are treated with the multipoint
method or with discontinuous elements by moving collocation points inside the
mesh cells whenever required.
"""
function assemble_system(mesh::Mesh{Dimension2},
                         mesh_bcs::MeshBoundaryCondition, order::Int64)
    # Discretize the mesh using Lagrangian elements specified by order. For now,
    # we only use geometrically linear boundary elements although our function
    # can be any order. Future implementation will allow for geometrically
    # curved boundary elements with B-splines (up to cubic). With a B-spline
    # boundary, we can also consider using a functional basis that matches the
    # B-spline basis. Although I don't see why the Lagrangian basis wouldn't
    # work on B-splines as well. This should be tested once both are
    # implemented.
    discretize_mesh!(mesh, order)

    A = initialize_matrix(mesh, mesh_bcs)

    # Use the collocation method to generate the system of equations. The
    # collocation points are chosen to be at the nodes of the mesh. We
    # introduce mulitnodes wherever there is a geometric discontinuity (corner).
    # If the mulitnode method results in redundant equations (as in the case of
    # Dirichlet boundary conditions specified on adjacent cells), the
    # collocation points are shifted inside each boundary cell to ensure
    # unique equations.
    index = 1
    for node in mesh.geometry.nodes
        # Get the edges connected to the current node.
        edge_indices = get_indices(mesh.topology.connectivity["N2E"], index)
        # Check if there are multiple Dirichlet conditions on the connected
        # edges. If so, we need to split and move the collocation point.
        if length(edge_indices) == 1 # This is an internal node
            collocate!(A, index, node, mesh.geometry, mesh_bcs)
        else
            # If there is not one edge attached to a node, there must be two
            # when working with a 2d mesh.
            @test length(edge_indices) == 2
            for edge in edge_indices
                # For the current node, find the nearest adjacent node on the
                # current edge
                node_indices = get_indices(mesh.topology.connectivity["E2N"], edge)
                for i in 1:length(node_indices)

                    if node_indices[i] == node
                        # Skip the current node. We want the closest node
                    else
                        if i == 1
                            # Take this as the minimum distance automatically
                            delta = mesh.geometry.nodes[i] - node
                        else
                            # The minimum distance, delta, must have been
                            # initialized by this point
                            @test isdefined(:delta)
                            test_delta = mesh.geometry.nodes[i] - node
                            if norm(test_delta) < norm(delta)
                                delta = test_delta
                            end
                        end
                    end
                end
                # Place a collocation point some fractional distance d between
                # the current node and the nearest internal node on the current
                # edge.
                d = 0.25 # Test with 0.25 for now
                coll_point = node + d*delta
                collocate!(A, index, coll_point, mesh.geometry, mesh_bcs)
            end
        end
        index += 1
    end
    return A
end


text:
using DickeModel
using Test

@testset "DickeModel.jl" begin
    # Write your tests here.
end


text:
using MuttsVersionedTrees
using Test
using Documenter

@testset "MuttsVersionedTrees" begin

    @testset "jldoctests" begin
    # Set up trees used in jldoctests
    Documenter.DocMeta.setdocmeta!(
        MuttsVersionedTrees,
        :DocTestSetup,
        quote
            using MuttsVersionedTrees
            using MuttsInterface
            head = VTree{Int}();

            head = insert!(head, 5);
            head = insert!(head, 1);
            head = insert!(head, 10);
            head = insert!(head, 2);

            println(head)

            mark_immutable!(head)

            head = insert!(head, 7);
            head = delete!(head, 2);

            mark_immutable!(head)
            head2 = branch!(head)

            head2 = insert!(head2, 2);
            head2 = insert!(head2, 8);
        end;
        recursive=true)

    Documenter.doctest(nothing, [MuttsVersionedTrees,])
    end
end


text:
@testset "asynchronous.jl" begin
    x, A = demodata()
    Ã = A .- mean(A, dims=1)

    Ψ = asynchronous(Ã)
    @plottest begin
        heatmap(x, x, Ψ)
        plot!(xflip=true, yflip=true)
        xlabel!("Peak position")
        ylabel!("Peak position")
        title!("Asynchronous")
    end "asynchronous.png"

    # Can this be prettier?
    @test .≈(diag(Ψ), 0, atol=2e-19) |> all
end


text:
struct Hubbard{A,B,N,O<:AbstractOperator{Float64}} <: Hamiltonian{A,Float64}
    address::A
    basis::B
    u::SMatrix{N,N,Float64}
    t::SVector{N,Float64}
    terms::O
end

function Hubbard(
    address;
    u::Union{AbstractMatrix{<:Real},Real}=1.0,
    t::Union{AbstractVector{<:Real},Real}=1.0,
    momentum=true,
    dispersion=hubbard_dispersion,
)
    basis = momentum ? MomentumSpace() : RealSpace()

    C = num_components(address)
    M = num_modes(address)
    # TODO here
    if u isa AbstractMatrix
        # Check that matrix is well-formed
        size(u) == (C, C) || throw(ArgumentError("`u` must be a $C×$C matrix or a scalar"))
        issymmetric(u) || throw(ArgumentError("`u` must be symmetric"))
        Rimu.Hamiltonians.warn_fermi_interaction(address, u)
        us = SMatrix{C,C,Float64}(u)
    else
        us = SMatrix{C,C,Float64}(fill(u, (C, C)))
    end
    # Create appropriate terms
    if momentum
        kinetic_term = ParticleCountTerm(KineticEnergyFunction(address, t, dispersion))
        interaction_term = MomentumTwoBodyTerm(InteractionMatrixFunction(us./M); fold=true)
    else
        if dispersion ≢ hubbard_dispersion
            throw(ArgumentError("setting dispersion is not supported in real space"))
        end
        kinetic_term = NeighbourOneBodyTerm(ParameterColumnFunction(address, t))
        interaction_term = OnsiteInteractionTerm(InteractionMatrixFunction(us))
    end
    terms = kinetic_term + interaction_term

    return Hubbard(address, basis, us, to_parameter_vector(address, t), terms)
end

function Base.show(io::IO, h::Hubbard)
    print(io, "Hubbard(")
    print(IOContext(io, :compact => true), h.address)
    print(io, ", u=$(h.u), t=$(h.t)")
    if basis(h) ≡ MomentumSpace()
        print(io, ", momentum=true)")
    else
        print(io, ", momentum=false)")
    end
end

starting_address(h::Hubbard) = h.address
LOStructure(::Hubbard) = IsHermitian()
terms(h::Hubbard) = h.terms
basis(h::Hubbard) = h.basis


text:
abstract type AbstractLeafNode{T, N, p} <: AbstractNode{T, N, p} end

Base.size(x::AbstractLeafNode) = size(typeof(x))
Base.IndexStyle(::Type{<: AbstractLeafNode}) = IndexLinear()

@inline function Base.getindex(x::AbstractLeafNode, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.data[i]
end

@inline function Base.setindex!(x::AbstractLeafNode, v, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds x.data[i] = v
    x
end

@generated function construct(::Type{T}) where {T}
    exps = [:(zero($T)) for T in fieldtypes(T)]
    quote
        @_inline_meta
        T($(exps...))
    end
end
@generated function allocate!(x::AbstractLeafNode{T}, i::Int) where {T}
    if T.mutable
        quote
            @_inline_meta
            @boundscheck checkbounds(x, i)
            @inbounds begin
                isactive(x, i) && return x[i]
                if isassigned(x.data, i)
                    leaf = unsafe_getindex(x.data, i)
                else
                    leaf = construct(T)
                end
                x[i] = leaf
            end
            Activated(x, i)
        end
    else
        quote
            @_inline_meta
            @boundscheck checkbounds(x, i)
            @inbounds getmask(x)[i] = true
            Activated(x, i)
        end
    end
end

function allocate!(x::AbstractLeafNode, mask::AbstractArray{Bool})
    checkbounds(x, CartesianIndices(mask))
    @simd for i in eachindex(mask)
        @inbounds mask[i] && allocate!(x, i)
    end
end

@inline function deactivate!(x::AbstractLeafNode)
    isnull(x) || fill!(getmask(x), false)
    x
end

isallocated(x::AbstractLeafNode, i::Int) = isassigned(x.data, i)
cleanup!(x::AbstractLeafNode) = x
nleaves(x::AbstractLeafNode) = count(getmask(x))

struct LeafNode{T, N, p} <: AbstractLeafNode{T, N, p}
    data::MaskedDenseArray{T, N}
    LeafNode{T, N, p}(data) where {T, N, p} = new(data)
    LeafNode{T, N, p}(::UndefInitializer) where {T, N, p} = new()
end

function LeafNode{T, N, p}() where {T, N, p}
    dims = size(LeafNode{T, N, p})
    data = MaskedDenseArray{T}(undef, dims)
    LeafNode{T, N, p}(data)
end

struct StructLeafNode{T, N, p, Ttuple} <: AbstractLeafNode{T, N, p}
    data::MaskedStructArray{T, N, Ttuple}
    StructLeafNode{T, N, p, Ttuple}(data) where {T, N, p, Ttuple} = new(data)
    StructLeafNode{T, N, p, Ttuple}(::UndefInitializer) where {T, N, p, Ttuple} = new()
end

function StructLeafNode{T, N, p, Ttuple}() where {T, N, p, Ttuple}
    dims = size(StructLeafNode{T, N, p})
    data = MaskedStructArray{T}(undef, dims)
    StructLeafNode{T, N, p, Ttuple}(data)
end

Base.propertynames(x::StructLeafNode) = (:data, propertynames(x.data)...)
@inline function Base.getproperty(x::StructLeafNode{T, N, p}, name::Symbol) where {T, N, p}
    name == :data && return getfield(x, :data)
    LeafNode{fieldtype(T, name), N, p}(getproperty(getfield(x, :data), name))
end

macro StructLeafNode(ex)
    @assert Meta.isexpr(ex, :braces)
    @assert length(ex.args) == 3
    T, N, p = ex.args
    esc(quote
        Ttuple = $StructArrays.map_params(t -> Array{t, $N}, $StructArrays.staticschema($T))
        StructLeafNode{$T, $N, $p, Ttuple}
    end)
end


text:
"""
    Epi(k,m)

compute ``\\sqrt{(x^2+m)}``
"""
function Epi(q, m2)
    sqrt(q^2 + m2)
end

"""
    star1fun(qp, qm, ps, m, T)

compute ``-\\frac{1}{\\pi}\\mathcal{F}_{1}\\left(q_{+}, q_{-}, p, \\bar{m}_{\\pi, k}^{2}\\right)``
"""
function star1fun(qp, qm, ps, m, T)
    (2 * π)^-2 *
    (4 * ps)^-1 *
    (
        Epi(qp, m) - Epi(qm, m) +
        2 * T * log((1 - exp(-Epi(qp, m) / T)) / (1 - exp(-Epi(qm, m) / T)))
    )
end

@doc raw"""
    star1funpm(qp, qm, ps, m, T)

compute $-\frac{1}{\pi}\mathcal{F}_{1}^{\prime}\left(q_{+}, q_{-}, p, \bar{m}_{\pi, k}^{2}\right)$

Compared with $\mathcal{F}$, the $\mathcal{F'}$ has subtracted vacuum contributions. it will be used in $\mathrm{Im} I_2$
"""
function star1funpm(qp, qm, ps, m, T)
    (2 * π)^-2 *
    (4 * ps)^-1 *
    (2 * T * log((1 - exp(-Epi(qp, m) / T)) / (1 - exp(-Epi(qm, m) / T))))
end

@doc raw"""
    star2fun(qp, qm, ps, k, m, T)

compute $-\frac{1}{\pi}\mathcal{F}_{2}\left(q_{+}, q_{-}, k, p, \bar{m}_{\pi, k}^{2}\right)$
"""
function star2fun(qp, qm, ps, k, m, T)
    (2 * π)^-2 *
    (4 * ps)^-1 *
    (
        Epi(k, m)^-1 * coth(Epi(k, m) / (2 * T)) * (1 / 2) * (k^2 - qm^2) + (
            Epi(qp, m) - Epi(k, m) +
            2 * T * log((1 - exp(-Epi(qp, m) / T)) / (1 - exp(-Epi(k, m) / T)))
        )
    )
end

@doc raw"""
    star2funpm(qp, qm, ps, k, m, T)

compute $-\frac{1}{\pi}\mathcal{F}^{\prime}_{2}\left(q_{+}, q_{-}, k, p, \bar{m}_{\pi, k}^{2}\right)$
"""
function star2funpm(qp, qm, ps, k, m, T)
    (
        (k - qm) * (k + qm) * (-1 + coth(Epi(k, m) / (2 * T))) +
        4 *
        T *
        Epi(k, m) *
        log((-1 + exp(-(Epi(qp, m) / T))) / (-1 + exp(-(Epi(k, m) / T))))
    ) / (32 * pi^2 * ps * Epi(k, m))
end

@doc raw"""
    star3fun(qp, ps, k, m, T)

compute $-\frac{1}{\pi}\mathcal{F}_{3}\left(q_{+},k, p, \bar{m}_{\pi, k}^{2}\right)$
"""
function star3fun(qp, ps, k, m, T)
    (2 * π)^-2 *
    (4 * Epi(k, m))^-1 *
    qp *
    coth(Epi(qp, m) / (2 * T)) *
    (1 - (qp^2 + ps^2 - k^2) / (2 * qp * ps))
end

@doc raw"""
    star3funpm(qp, ps, k, m, T)

compute $-\frac{1}{\pi}\mathcal{F}^{\prime}_{3}\left(q_{+}, k, p, \bar{m}_{\pi, k}^{2}\right)$
"""
function star3funpm(qp, ps, k, m, T)
    ((k + ps - qp) * (k - ps + qp) * (-1 + coth(Epi(qp, m) / (2 * T)))) /
    (32 * pi^2 * ps * Epi(k, m))
end

@doc raw"""
    star4fun(p0, p, k, m, T)

compute $-\frac{1}{\pi}\mathcal{F}_{4}\left(p_{0}, k, p, \bar{m}_{\pi, k}^{2}\right)$
"""
function star4fun(p0, p, k, m, T, dϵ = 0.02)
    if k > p / 2
        return ((-2 * k + p)^2 * (4 * k + p) * coth(Epi(k, m) / (2 * T))) /
               (384 * pi^2 * Epi(k, m)^2) * deltafun(p0 - 2 * Epi(k, m), dϵ)
    else
        return 0.0
    end
end

@doc raw"""
    loopfunppfix(p0, ps, k, m, T)

compute $-\frac{1}{\pi}\Im I_{1, k}(p)$
"""
function loopfunppfix(p0, ps, k, m, T, dϵ = 0.02)
    star4fun(p0, ps, k, m, T, dϵ) + loopfunpp(p0, ps, k, m, T)
end

@doc raw"""
    loopfunpm(p0, ps, k, m, T)

compute $-\frac{1}{\pi}\Im I_{2, k}(p)$
"""
function loopfunpm(p0, ps, k, m, T)
    if m >= 0
        if k > ps
            if p0 > ps
                return 0.0
            elseif Epi(k + ps, m) - Epi(k, m) < p0 <= ps
                qm =
                    -ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1funpm(qp, qm, ps, m, T)
            elseif p0 <= Epi(k + ps, m) - Epi(k, m)
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            end
        elseif ps / 2 < k <= ps
            if p0 > ps
                return 0.0
            elseif Epi(k + ps, m) - Epi(k, m) < p0 <= ps
                qm =
                    -ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1funpm(qp, qm, ps, m, T)
            elseif Epi(ps, m) - Epi(k, m) < p0 <= Epi(k + ps, m) - Epi(k, m)
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            elseif p0 <= Epi(ps, m) - Epi(k, m)
                qm = ps - sqrt((p0 + Epi(k, m))^2 - m)
                qp = sqrt((p0 + Epi(k, m))^2 - m)
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            end
        elseif k <= ps / 2
            if p0 > ps
                return 0.0
            elseif Epi(k + ps, m) - Epi(k, m) < p0 <= ps
                qm =
                    -ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1funpm(qp, qm, ps, m, T)
            elseif Epi(ps, m) - Epi(k, m) < p0 <= Epi(k + ps, m) - Epi(k, m)
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            elseif Epi(k - ps, m) - Epi(k, m) < p0 <= Epi(ps, m) - Epi(k, m)
                qm = ps - sqrt((p0 + Epi(k, m))^2 - m)
                qp = sqrt((p0 + Epi(k, m))^2 - m)
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            elseif p0 <= Epi(k - ps, m) - Epi(k, m)
                qm =
                    ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1funpm(qp, qm, ps, m, T)
            end
        end
    elseif m < 0
        if k > ps
            if p0 > Epi(k + ps, m) - Epi(k, m)
                return 0.0
            elseif ps < p0 <= Epi(k + ps, m) - Epi(k, m)
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                qmp =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (-4 * m + p0^2 - ps^2)) /
                    (2 * (p0^2 - ps^2))
                qpp = qmp + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T) -
                       star1funpm(qpp, qmp, ps, m, T)
            elseif p0 <= ps
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            end
        elseif ps / 2 < k <= ps
            if p0 > Epi(k + ps, m) - Epi(k, m)
                return 0.0
            elseif ps < p0 <= Epi(k + ps, m) - Epi(k, m)
                qm3 = sqrt(-m + (sqrt(k^2 + m) + p0)^2) - ps
                qm3p =
                    (
                        -ps +
                        sqrt(p0^2 * (-p0^2 + ps^2) * (4 * m - p0^2 + ps^2)) /
                        (p0^2 - ps^2)
                    ) / 2
                qp3 = sqrt(-m + (sqrt(k^2 + m) + p0)^2)
                qp3p =
                    (
                        ps +
                        sqrt(p0^2 * (-p0^2 + ps^2) * (4 * m - p0^2 + ps^2)) /
                        (p0^2 - ps^2)
                    ) / 2
                return -star1funpm(qp3p, qm3p, ps, m, T) +
                       star2funpm(qp3, qm3, ps, k, m, T) -
                       star3funpm(qp3, ps, k, m, T)
            elseif Epi(ps, m) - Epi(k, m) < p0 <= ps
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            elseif p0 <= Epi(ps, m) - Epi(k, m)
                qm = ps - sqrt((p0 + Epi(k, m))^2 - m)
                qp = sqrt((p0 + Epi(k, m))^2 - m)
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            end
        elseif k <= ps / 2
            if p0 > Epi(k + ps, m) - Epi(k, m)
                return 0.0
            elseif ps < p0 <= Epi(k + ps, m) - Epi(k, m)
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                qmp =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (-4 * m + p0^2 - ps^2)) /
                    (2 * (p0^2 - ps^2))
                qpp = qmp + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T) -
                       star1funpm(qpp, qmp, ps, m, T)
            elseif Epi(ps, m) - Epi(k, m) < p0 <= ps
                qm = -ps + sqrt((p0 + Epi(k, m))^2 - m)
                qp = qm + ps
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            elseif Epi(k - ps, m) - Epi(k, m) < p0 <= Epi(ps, m) - Epi(k, m)
                qm = ps - sqrt((p0 + Epi(k, m))^2 - m)
                qp = sqrt((p0 + Epi(k, m))^2 - m)
                return star2funpm(qp, qm, ps, k, m, T) -
                       star3funpm(qp, ps, k, m, T)
            elseif p0 <= Epi(k - ps, m) - Epi(k, m)
                qm =
                    ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1funpm(qp, qm, ps, m, T)
            end
        end
    end
end

function compensate_delta(ps, k, m, T)
    ((2 * k - ps) * coth(Epi(k, m) / (2 * T))) / (16 * pi^2 * Epi(k, m))
end

@doc raw"""
    loopfunpp(p0, ps, k, m, T)

compute $-\frac{1}{\pi}\Im I_{2, k}(p)$, but without $\mathcal{F}_4$
contribution
"""
function loopfunpp(p0, ps, k, m, T)
    if m >= 0
        if k > ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif 2 * Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif p0 < 2 * Epi(k, m)
                return 0.0
            end
        elseif ps / 2 < k <= ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif 2 * Epi(k, m) <= p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif p0 < 2 * Epi(k, m)
                return 0.0
            end
        elseif k <= ps / 2
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif Epi(k, m) + Epi(ps - k, m) <= p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif 2 * Epi(ps / 2, m) <= p0 < Epi(k, m) + Epi(ps - k, m)
                qm =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1fun(qp, qm, ps, m, T)
            elseif p0 < 2 * Epi(ps / 2, m)
                return 0.0
            end
        end
    elseif m < 0
        if k > ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif 2 * Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif p0 < 2 * Epi(k, m)
                return 0.0
            end
        elseif ps / 2 < k <= ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif 2 * Epi(k, m) <= p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif p0 < 2 * Epi(k, m)
                return 0.0
            end
        elseif k <= ps / 2
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif 2 * Epi(ps / 2, m) <= p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif Epi(k, m) + Epi(ps - k, m) <= p0 < 2 * Epi(ps / 2, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                qpp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qmp =
                    ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star2fun(qp, qm, ps, k, m, T) +
                       star3fun(qp, ps, k, m, T) - star1fun(qpp, qmp, ps, m, T)
            elseif p0 < Epi(k, m) + Epi(ps - k, m)
                return 0.0
            end
        end
    end
end

function loopfunppSmooth(p0, ps, k, m, T,slop=5)
    if m >= 0
        if k > ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif 1.7*Epi(k, m)<p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return 0.5*(star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T))*(tanh(slop*(p0-2 * Epi(k, m)))+1)
            elseif p0<=1.7*Epi(k, m)
                return 0.0
            end
        elseif ps / 2 < k <= ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif Epi(k, m)<p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return 0.5*(star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T))*(tanh(slop*(p0-2 * Epi(k, m)))+1)
            elseif p0<=Epi(k, m)
                return 0.0
            end
        elseif k <= ps / 2
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif Epi(k, m) + Epi(ps - k, m) <= p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif 2 * Epi(ps / 2, m) <= p0 < Epi(k, m) + Epi(ps - k, m)
                qm =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp =
                    ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star1fun(qp, qm, ps, m, T)
            elseif p0 < 2 * Epi(ps / 2, m)
                return 0.0
            end
        end
    elseif m < 0
        if k > ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return 0.5*(star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T))*(tanh(slop*(p0-2 * Epi(k, m)))+1)
            end
        elseif ps / 2 < k <= ps
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return 0.5*(star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T))*(tanh(slop*(p0-2 * Epi(k, m)))+1)
            end
        elseif k <= ps / 2
            if p0 >= Epi(k + ps, m) + Epi(k, m)
                qm =
                    -ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qp = qm + ps
                return star1fun(qp, qm, ps, m, T)
            elseif Epi(ps, m) + Epi(k, m) <= p0 < Epi(k + ps, m) + Epi(k, m)
                qm = sqrt((p0 - Epi(k, m))^2 - m) - ps
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif 2 * Epi(ps / 2, m) <= p0 < Epi(ps, m) + Epi(k, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                return star2fun(qp, qm, ps, k, m, T) + star3fun(qp, ps, k, m, T)
            elseif Epi(k, m) + Epi(ps - k, m) <= p0 < 2 * Epi(ps / 2, m)
                qm = ps - sqrt((p0 - Epi(k, m))^2 - m)
                qp = sqrt((p0 - Epi(k, m))^2 - m)
                qpp =
                    ps / 2 -
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                qmp =
                    ps / 2 +
                    sqrt(p0^2 * (p0^2 - ps^2) * (p0^2 - ps^2 - 4 * m)) /
                    (2 * (p0^2 - ps^2))
                return star2fun(qp, qm, ps, k, m, T) +
                       star3fun(qp, ps, k, m, T) - star1fun(qpp, qmp, ps, m, T)
            elseif p0 < Epi(k, m) + Epi(ps - k, m)
                return 0.0
            end
        end
    end
end


text:
using FactCheck
using SimpleHttpCommon
using StreamReader
const S = SimpleHttpCommon

if !isdefined(:TEST_DIR)
    include("utils.jl")
end

facts("Test Response.jl") do
    context("Basic") do
        context("Initial") do
            r = mkres()
            @fact r.headers_sent --> false
            @fact r.writer == S.startwriter --> true
            @fact r.content_length --> S.EMPTY_CL
            S.start(r)
            @fact_throws S.start(r)
            @fact r.headers_sent --> true
            @fact r.content_length --> 0

            expected_headers = ["Content-Length",
                                "Content-Type",
                                "Date",
                                "Server"]
            for k in keys(r.headers)
                @fact k in expected_headers --> true "Key '$k' not found"
            end

            @fact r.headers["Content-Length"] --> "0"
        end

        context("Test contentlength") do
            r = mkres()
            @fact r.content_length --> -1
            r.content_length = 10
            @fact r.content_length --> 10
        end

        context("Test sendheaders") do
            context("basic") do
                r = mkres()
                r.headers = S.headers()
                @fact r.content_length --> S.EMPTY_CL
                r.content_length = 10
                @fact r.content_length --> 10
                S.sendheaders(r)
                edata = "HTTP/1.1 200 OK\r\nContent-Length: 10\r\n\r\n"
                @fact str(r.io.io.data) --> edata
                @fact r.headers_size --> sizeof(edata)
                @fact r.headers["Content-Length"] --> "10"
            end

            context("response status not accepts body") do
                r = mkres()
                r.content_length = 10
                r.status = r.protocol.status_wo_msg[1]

                fre(;key="status_not_accept_msg_body",
                    header="Content-Length",
                    status=r.protocol.status_wo_msg[1]) do
                    S.sendheaders(r)
                end
            end

            context("chunked") do
                context("manual") do
                    r = mkres()
                    r.headers = ["Transfer-Encoding" => "chunked"]
                    S.sendheaders(r)
                    @fact r.headers_sent --> true
                    @fact r.content_length --> -1

                    r = mkres()
                    r.headers = ["Transfer-Encoding" => "chunked", "Content-Length" => "10"]

                    fre(;key="invalid_header",
                        header="Content-Length") do
                        S.sendheaders(r)
                    end
                end

                context("chunk()") do
                    r = mkres()
                    chunk(r)
                    S.sendheaders(r)
                    @fact r.headers_sent --> true
                    @fact r.content_length --> -1

                    r = mkres()
                    r.headers = ["Transfer-Encoding" => "chunked", "Content-Length" => "10"]

                    fre(;key="invalid_header",
                        header="Content-Length") do
                        S.sendheaders(r)
                    end
                end
            end
        end
    end

    context("Basic Response writer") do
        context("writer functions") do
            function t(fn::Function, title)
                l = {"God is love!", "Jesus is my best Friend!!"}
                l = {(l[1],l[1]), (l[2].data,l[2])}
                for (data, data_r) in l
                    context(@sprintf "The %s for write '%s'" title data) do
                        r = mkres()
                        cl = sizeof(data)
                        fn(r, data)
                        @fact r.headers_sent --> true
                        @fact sentcl(r) --> cl
                        @fact r.content_length --> cl
                        @fact r.headers["Content-Length"] --> "$cl"
                        @fact str(r.io.io.data)[end-(cl-1):end] --> data_r
                    end
                end
            end

            t("Base.write") do r, data
                write(r, data)
            end
            t("'data >> response'") do r, data
                data >> r
            end
            t("'response << data'") do r, data
                r << data
            end
        end

        context("Simple writer") do
            context("unique data block") do
                r = mkres()
                data = "God is love!"
                write(r, data)
                @fact r.headers_sent --> true
                cl = 12
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact r.headers["Content-Length"] --> "$cl"
                @fact str(r.io.io.data)[end-(cl-1):end] --> data

                extrapolates = "Jesus Christ is God in human form 1!"

                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)) do
                    write(r, extrapolates)
                end
            end

            context("multiples data block") do
                data = ["God is love!", "Jesus is my King!"]
                size = [sizeof(_) for _ in data]
                cl = sum(size)
                r = mkres()
                r.content_length = cl
                write(r, data[1])
                @fact r.headers_sent --> true
                @fact sentcl(r) --> size[1]
                @fact r.content_length --> cl
                @fact r.headers["Content-Length"] --> "$cl"
                @fact str(r.io.io.data)[end-(size[1]-1):end] --> data[1]
                write(r, data[2])
                @fact str(r.io.io.data)[end-(cl-1):end] --> join(data, "")
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact r.headers["Content-Length"] --> "$cl"

                extrapolates = "Jesus Christ is God in human form 2!"
                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)) do
                    write(r, extrapolates)
                end
            end

            context("multiples data block with partial extrapolation") do
                data = ["God is love!", "Jesus is my King!"]
                size = [sizeof(_) for _ in data]
                cl = sum(size) + 3
                r = mkres()
                r.content_length = cl

                write(r, data[1])
                @fact r.headers_sent --> true
                @fact sentcl(r) --> size[1]
                @fact r.content_length --> cl
                @fact r.headers["Content-Length"] --> "$cl"
                @fact str(r.io.io.data)[end-(size[1]-1):end] --> data[1]

                write(r, data[2])
                @fact str(r.io.io.data)[end-(cl-4):end] --> join(data, "")
                @fact sentcl(r) --> (cl-3)
                @fact r.content_length --> cl
                @fact r.headers["Content-Length"] --> "$cl"

                extrapolates = "Jesus!!"

                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)-3) do
                    write(r, extrapolates)
                end

                @fact str(r.io.io.data)[end-(cl-1):end] --> string(join(data, ""), "\0\0\0")
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact r.headers["Content-Length"] --> "$cl"
            end
        end

        context("Buffered writer") do
            context("unique block") do
                r = mkres()
                @fact r.writer --> exactly(S.startwriter)
                buffered(r)
                @fact r.writer --> exactly(S.startwriter)

                data = "God is love!"
                cl = sizeof(data)

                r << data
                @fact r.writer --> exactly(S.bufferedwriter)
                @fact r.headers_sent --> false # is in buffer
                @fact sentcl(r) --> 0
                @fact r.content_length --> S.EMPTY_CL

                flush(r)

                @fact r.writer --> exactly(S.bufferedwriter)
                @fact r.headers_sent --> true
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact str(r.io.io.data)[end-(cl-1):end] --> data

                extrapolates = "Jesus Christ is God in human form 3!"

                # write in buffer
                r << extrapolates

                @fact str(r.buffer.data) --> extrapolates

                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)) do
                    flush(r)
                end
            end
            context("unique block with set contentlength") do
                r = mkres()
                buffered(r)

                data = "God is love!"
                cl = sizeof(data)
                r.content_length = cl

                r << data

                @fact r.headers_sent --> false # is in buffer
                @fact sentcl(r) --> 0
                @fact r.content_length --> cl

                flush(r)

                @fact r.headers_sent --> true
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact str(r.io.io.data)[end-(cl-1):end] --> data

                extrapolates = "Jesus Christ is God in human form 4!"

                # write in buffer
                r << extrapolates

                @fact str(r.buffer.data) --> extrapolates

                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)) do
                    flush(r)
                end
            end
            context("multiples blocks") do
                r = mkres()
                buffered(r)

                data = ["God is", " love!"]
                size = [sizeof(_) for _ in data]
                cl = sum(size)

                r << data[1]

                @fact r.headers_sent --> false # is in buffer
                @fact sentcl(r) --> 0
                @fact r.content_length --> S.EMPTY_CL
                @fact str(r.buffer.data) --> data[1]
                @fact sizeof(r.io.io.data) --> 0

                r << data[2]

                @fact r.headers_sent --> false # is in buffer
                @fact sentcl(r) --> 0
                @fact r.content_length --> S.EMPTY_CL
                @fact str(r.buffer.data) --> join(data, "")
                @fact sizeof(r.io.io.data) --> 0

                flush(r)

                @fact r.headers_sent --> true
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact str(r.io.io.data)[end-(cl-1):end] --> join(data, "")

                extrapolates = "Jesus Christ is God in human form 5!"

                # write in buffer
                r << extrapolates

                @fact str(r.buffer.data) --> extrapolates

                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)) do
                    flush(r)
                end
            end
            context("multiples blocks with set contentlength") do
                r = mkres()
                buffered(r)

                data = ["God is", " love!"]
                size = [sizeof(_) for _ in data]
                cl = sum(size)

                r.content_length = cl

                write(r, data[1])

                @fact r.headers_sent --> false # is in buffer
                @fact sentcl(r) --> 0
                @fact r.content_length --> cl
                @fact str(r.buffer.data) --> data[1]
                @fact sizeof(r.io.io.data) --> 0

                flush(r)

                @fact r.headers_sent --> true
                @fact sentcl(r) --> size[1]
                @fact r.content_length --> cl
                @fact sizeof(r.buffer.data) --> 0
                @fact str(r.io.io.data)[end-(size[1]-1):end] --> data[1]
                @fact sizeof(r.io.io.data) --> (r.headers_size + size[1])

                write(r, data[2])

                @fact r.headers_sent --> true
                @fact sentcl(r) --> size[1]
                @fact r.content_length --> cl
                @fact str(r.buffer.data) --> data[2]
                @fact str(r.io.io.data)[end-(size[1]-1):end] --> data[1]
                @fact (sizeof(r.io.io.data)-r.headers_size) --> size[1]

                flush(r)

                @fact r.headers_sent --> true
                @fact sentcl(r) --> cl
                @fact r.content_length --> cl
                @fact str(r.buffer.data) --> ""
                @fact str(r.io.io.data)[end-(cl-1):end] --> join(data, "")
                @fact (sizeof(r.io.io.data)-r.headers_size) --> cl

                extrapolates = "Jesus Christ is God in human form 6!"

                # write in buffer
                write(r, extrapolates)

                @fact str(r.buffer.data) --> extrapolates

                fre(;key="exceeded_content_data",
                    content_length=cl,
                    data_size=sizeof(extrapolates),
                    exceeded=sizeof(extrapolates)) do
                    flush(r)
                end
            end
        end

        context("chunked writer") do
            r = mkres()
            r.headers = S.headers()
            chunk(r)

            blocks = String["For God so loved the world, that he gave his ",
                    "only begotten Son, that whosoever ",
                    "believeth in him should not perish, but ",
                    "have everlasting life.",
                    "(John 3:16)"]

            for b in blocks
                b >> r
            end

            S.done(r)

            io = r.io
            seekstart(io.io)

            edata = "HTTP/1.1 200 OK\r\nTransfer-Encoding: chunked\r\n\r\n"
            @fact UTF8String(readbytes(io.io, r.headers_size)) --> edata
            @fact String[UTF8String(_) for _ in collect(S.Chunk.iterator(io))] --> blocks
        end

        context("chunked writer with buffer") do
            r = mkres()
            r.headers = S.headers()

            chunk(r)
            buffered(r)

            blocks = String["For God so loved the world, that he gave his ",
                    "only begotten Son, that whosoever ",
                    "believeth in him should not perish, but ",
                    "have everlasting life.",
                    "(John 3:16)"]

            writeblock = (i, b) -> begin
                s = sizeof(b)
                total_parts, parts_size, last_part_size = calculate(s, 7)
                prev_pn = 1
                start_ = 0

                for pn in 1:(total_parts-1)
                    pd = b[(start_+1):start_ + parts_size]
                    @fact sizeof(pd) --> parts_size
                    r << pd
                    prev_pn = pn
                    start_ += parts_size
                end

                pd = b[(start_+1):end]
                @fact sizeof(pd) --> last_part_size
                r << pd

                flush(r)
            end

            for (i, b) in enumerate(blocks)
                writeblock(i, b)
            end

            S.done(r)

            io = r.io
            seekstart(io.io)

            edata = "HTTP/1.1 200 OK\r\nTransfer-Encoding: chunked\r\n\r\n"
            @fact UTF8String(readbytes(io.io, r.headers_size)) --> edata
            @fact String[UTF8String(_) for _ in collect(S.Chunk.iterator(io))] --> blocks
        end

        context("chunked on done") do
            r = mkres()
            r.headers = S.headers()
            chunk(r)

            blocks = String["For God so loved the world, that he gave his ",
                    "only begotten Son, that whosoever ",
                    "believeth in him should not perish, but ",
                    "have everlasting life.",
                    "(John 3:16)"]

            r.data = blocks
            S.done(r)

            io = r.io
            seekstart(io.io)

            edata = "HTTP/1.1 200 OK\r\nTransfer-Encoding: chunked\r\n\r\n"
            @fact UTF8String(readbytes(io.io, r.headers_size)) --> edata
            @fact String[UTF8String(_) for _ in collect(S.Chunk.iterator(io))] --> blocks
        end
    end
end


text:
using GOESRLevel0
using Documenter

makedocs(;
    sitename="GOESRLevel0.jl",
    authors="Ben Elsen <mail@benelsen.com>",
    modules=[GOESRLevel0],
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://benelsen.github.io/GOESRLevel0.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Guide" => "guide.md",
        "Data" => "data.md",
        "Library" => "library.md",
    ],
    repo="https://github.com/benelsen/GOESRLevel0.jl/blob/{commit}{path}#L{line}",
)

deploydocs(;
    repo="github.com/benelsen/GOESRLevel0.jl",
)


text:
function solution()
  sum = 0

  for x = 1:999
    if x % 3 == 0 || x % 5 == 0
      sum += x
    end
  end

  sum
end

# XXX Shorter, but 100%+ slower. std_dev is way higher as well
#=function solution()=#
  #=sum([x % 3 == 0 || x % 5 == 0 ? x : 0 for x = 0:999])=#
#=end=#

if length(ARGS) == 1 && ARGS[1] == "-a"
  println(solution())
  exit(0)
end

for line in eachline(STDIN)
  iters = int(line)

  start = time_ns()
  for i in 1:iters
    solution()
  end
  end_ = time_ns()

  println(end_ - start)
end


text:
"""
    isdegenerate(N::AbstractEcologicalNetwork)

Networks are called degenerate if some species have no interactions, either at
all, or with any species other than themselves. This is particularly useful to
decide the networks to keep when generating samples for null models.
"""
function isdegenerate(N::AbstractEcologicalNetwork)
    return minimum(values(degree(nodiagonal(N)))) == 0.0
end

"""
    simplify{T<:AbstractBipartiteNetwork}(N::T)

Returns a new network in which species with no interactions have been removed.
"""
function simplify(N::T) where {T<:AbstractBipartiteNetwork}
    Y = copy(N)
    simplify!(Y)
    return Y
end

"""
    simplify!{T<:AbstractBipartiteNetwork}(N::T)

Returns a new network in which species with no interactions have been removed.
"""
function simplify!(N::T) where {T<:AbstractBipartiteNetwork}
    d1 = degree(N; dims=1)
    d2 = degree(N; dims=2)
    s1 = species(N; dims=1)
    s2 = species(N; dims=2)
    p1 = filter(i -> d1[s1[i]] > zero(eltype(N)[1]), 1:richness(N; dims=1))
    p2 = filter(i -> d2[s2[i]] > zero(eltype(N)[1]), 1:richness(N; dims=2))
    N.T = N.T[p1]
    N.B = N.B[p2]
    N.A = N.A[p1,p2]
end

"""
    simplify(N::AbstractUnipartiteNetwork)

Returns a new network in which species with no interactions have been removed.
"""
function simplify(N::T) where {T <: AbstractUnipartiteNetwork}
    Y = copy(N)
    simplify!(Y)
    return Y
end

"""
    simplify!(N::AbstractUnipartiteNetwork)

Modifies the network to drop all species without an interaction.
"""
function simplify!(N::T) where {T <: AbstractUnipartiteNetwork}
    d = degree(N)
    s = species(N)
    positions = filter(i -> d[s[i]] > zero(eltype(N)[1]), 1:richness(N))
    N.S = N.S[positions]
    N.A = N.A[positions, positions]
end


text:
using Base: IEEEFloat

#=

  - error-free transformations (`two_sum`, `two_diff`, `two_square`, `two_prod`,
                                `two_hilo_sum`, `two_lohi_sum`, `two_hilo_diff`, `two_lohi_diff`)
  - least-error transformations (`two_sqrt`, `two_inv`, `two_div`)

 nomenclature:
 - `two_<op>`: the "two" refers to the number of values returned
    two and three argument versions of `two_[sum,diff,prod]` are given
 - `two_hilo_<op>` the argments are ordered by non-increasing absolute value
    - this is `fast_two_sum` in the literature
 - `two_lohi_<op>` the argments are ordered by non-decreasing absolute value
    - this is not available in the literature

=#

"""
    two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_sum(a::T, b::T) where {T}
    hi = a + b
    v  = hi - a
    lo = (a - (hi - v)) + (b - v)
    return hi, lo
end

"""
   two_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
- Unchecked Precondition: !(isinf(a) | isinf(b) | isinf(c))
"""
function two_sum(a::T, b::T, c::T) where {T}
    a, b, c = magnitude_mintomax(a, b, c)
    md, lo = two_sum(b, c)
    hi, md = two_sum(a, md)
    hi, lo = two_hilo_sum(hi, md+lo)
    return hi, lo
end

"""
    two_diff(a, b)

Computes `s = fl(a-b)` and `e = err(a-b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_diff(a::T, b::T) where {T}
    hi = a - b
    v  = hi - a
    lo = (a - (hi - v)) - (b + v)
    return hi, lo
end

"""
    two_diff(a, b, c)

Computes `s = fl(a-b-c)` and `e1 = err(a-b-c), e2 = err(e1)`.
- Unchecked Precondition: !(isinf(a) | isinf(b) | isinf(c))
"""
function two_diff(a::T,b::T,c::T) where {T}
    a, b, c = magnitude_maxtomin(a, b, c)
    s, t = two_diff(-b, c)
    hi, u = two_sum(a, s)
    hi, lo = two_hilo_sum(hi, u+t)
    return hi, lo
end

"""
    two_hilo_sum(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_hilo_sum(a::T, b::T) where {T}
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

"""
    two_lohi_sum(a, b)

*unchecked* requirement `|b| ≥ |a|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_lohi_sum(a::T, b::T) where {T}
    hi = b + a
    lo = a - (hi - b)
    return hi, lo
end

"""
    two_hilo_diff(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a-b)` and `lo = err(a-b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_hilo_diff(a::T, b::T) where {T}
    hi = a - b
    lo = (a - hi) - b
    hi, lo
end

"""
    two_lohi_diff(a, b)

*unchecked* requirement `|b| ≥ |a|`
Computes `hi = fl(a-b)` and `lo = err(a-b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_lohi_diff(a::T, b::T) where {T}
    hi = b - a
    lo = (b - hi) - a
    hi, lo
end

"""
    two_square(a)

Computes `hi = fl(a*a)` and `lo = fl(err(a*a))`.
- Unchecked Precondition: !(isinf(a))
"""
@inline function two_square(a::T) where {T}
    hi = a * a
    lo = fma(a, a, -hi)
    hi, lo
end

"""
    two_prod(a, b)

Computes `hi = fl(a*b)` and `lo = fl(err(a*b))`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_prod(a::T, b::T) where {T}
    hi = a * b
    lo = fma(a, b, -hi)
    hi, lo
end

"""
    two_prod(a, b, c)

Computes `hi = fl(a*b*c)` and `lo = err(a*b*c)`.
- Unchecked Precondition: !(isinf(a) | isinf(b) | isinf(c))
"""
function two_prod(a::T, b::T, c::T) where {T}
    abhi, ablo = two_prod(a, b)
    hi, abhiclo = two_prod(abhi, c)
    ablochi, abloclo = two_prod(ablo, c)
    lo = ablochi + (abhiclo + abloclo)
    return hi, lo
end

"""
    two_inv(a)

Computes `hi = fl(inv(a))` and `lo = err(inv(a))`.
- Unchecked Precondition: !(isinf(a))
"""
@inline function two_inv(a::T) where {T}
     hi = inv(a)
     lo = fma(-hi, a, one(T))
     lo /= a
     return hi, lo
end

"""
    two_div(a, b)

Computes `hi = fl(a/b)` and `lo = err(a/b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_div(a::T, b::T) where {T}
     hi = a / b
     lo = fma(-hi, b, a)
     lo /= b
     return hi, lo
end

"""
    two_sqrt(a)

Computes `hi = fl(sqrt(a))` and `lo = err(sqrt(a))`.
- Unchecked Precondition: !(isinf(a))
"""
@inline function two_sqrt(a::T) where {T}
    hi = sqrt(a)
    lo = fma(-hi, hi, a)
    lo /= 2
    lo /= hi
    return hi, lo
end

Base.:(+)(xhi::T, xlo::T, y::T) where {T<:IEEEFloat} = two_sum(xhi, y, xlo)
Base.:(-)(xhi::T, xlo::T, y::T) where {T<:IEEEFloat} = two_sum(xhi, xlo, -y)

Base.:(+)(x::Tuple{T,T}, y::T) where {T<:IEEEFloat} = (+)(x[1], x[2], y)
Base.:(-)(x::Tuple{T,T}, y::T) where {T<:IEEEFloat} = (+)(x[1], x[2], y)

# Algorithm 12 from Tight and rigourous error bounds.  relative error <= 5u²
function Base.:(*)(xhi::T, xlo::T, y::T) where {T<:IEEEFloat}
    hi, lo = two_prod(xhi, y)
    t = fma(xlo, y, lo)
    hi, lo = two_hilo_sum(hi, t)
    return hi, lo
end

function Base.:(/)(xhi::T, xlo::T, y::T) where {T<:IEEEFloat}
    hi = xhi / y
    uh, ul = two_prod(hi, y)
    lo = ((((xhi - uh) - ul) + xlo))/y
    hi,lo = two_hilo_sum(hi, lo)
    return hi, lo
end

Base.:(*)(x::Tuple{T,T}, y::T) where {T<:IEEEFloat} = (*)(x[1], x[2], y)
Base.:(/)(x::Tuple{T,T}, y::T) where {T<:IEEEFloat} = (/)(x[1], x[2], y)

# exchange sort for magnitudes, 2 or 3 values

@inline magnitude_minmax(a, b) = abs(a) < abs(b) ? (a, b) : (b, a)

@inline function magnitude_mintomax(a, b, c)
    b, c = magnitude_minmax(b, c)
    a, c = magnitude_minmax(a, c)
    a, b = magnitude_minmax(a, b)
    return a, b, c
end