text:
using RecipesBase

import LinearAlgebra: norm, dot, cross
import StaticArrays: SVector

export VortexFilament, Segment, getfreevertices, getboundvertices, inducevelocity, issemiinf

const Segment = SVector{2,AbstractArray}
const eps = 1e-10

"""
$(TYPEDEF)

Defines a vortex filament of strength `Γ`, discretized by vertices in `vertices` and segments in `segments`, which connect the vertices. Each segment in `segments` should connect two subsequent entries of `vertices`. A vertex can be marked as a bound (e.g. to a wing) if its index appears in `boundidx`. Similarly, a vortex is marked as a free vertex if its index appears in `freeidx`.

# Fields

$(FIELDS)
"""
struct VortexFilament{TV<:AbstractVector{<:AbstractVector},TI<:AbstractVector{Int}}
    """Γ: Strength of the vortex filament."""
    Γ::Real
    """vertices: Array of vertices of the vortex filament."""
    vertices::TV
    """segments: Array of segments of the vortex filament."""
    segments::Vector{Segment}
    """freeidx: Array of indices of the free vertices of the vortex filament."""
    freeidx::TI
    """boundidx: Array of indices of the bound vertices of the vortex filament."""
    boundidx::TI
end

"""
$(SIGNATURES)

Constructs a `VortexFilament` with strength `Γ` and vertices `v`. `boundidx` can be used to provide the indices of vertices that are considered as 'bound'. The constructor then set `freeidx` of the `VortexFilament` as the complementary list of indices. If the keyword `isclosed` is `true`, the
`segments` of the `VortexFilament` will contain a segment that connects the first and last vertex in `v`, provided that these do not contain an infinite coordinate.
"""
function VortexFilament(Γ::Real, v::TV, boundidx::TI=Int64[]; isclosed::Bool=true) where {TV<:AbstractVector{<:AbstractVector}, TI<:AbstractVector{Int}}
    infbools = map(v->in(Inf,abs.(v)),v) # boolean array with the i-th element true if the i-th element of v lies at infinity

    segments = Segment.(v[1:end-1],v[2:end])

    if !any(infbools[[1,end]]) && isclosed # if there is no vertex at infinity, close the loop by adding a segment that connects the last and first vertex
        push!(segments,Segment(v[end],v[1]))
    end

    boundidx = unique(boundidx) # remove any duplicates from boundidx
    freeidx = setdiff(1:length(v),boundidx)

    return VortexFilament{TV,TI}(Γ,v,segments,freeidx,boundidx)
end

"""
$(TYPEDSIGNATURES)

Returns the free vertices of the vortex filament `vf`.
"""
function getfreevertices(vf::VortexFilament)
    v = vf.vertices[vf.freeidx]
    return v
end

"""
$(TYPEDSIGNATURES)

Returns the bound vertices of the vortex filament `vf`.
"""
function getboundvertices(vf::VortexFilament)
    v = vf.vertices[vf.boundidx]
    return v
end

"""
$(TYPEDSIGNATURES)

Checks if `vf` is an infinite vortex filament.
"""
@inline function Base.:isinf(vf::VortexFilament)
    any(isinf,vf.segments) && return true
    semiinfsegs = findall(issemiinf,vf.segments)
    return length(semiinfsegs) > 1
end

"""
$(TYPEDSIGNATURES)

Checks if `s` is an infinite segment.
"""
@inline function Base.:isinf(s::Segment)
    return any(isinf,s[1]) && any(isinf,s[2])
end

"""
$(TYPEDSIGNATURES)

Checks if `vf` is a semi-infinite vortex filament.
"""
@inline function issemiinf(vf::VortexFilament)
    semiinfsegs = findall(issemiinf,vf.segments)
    return length(semiinfsegs) == 1
end

"""
$(TYPEDSIGNATURES)

Checks if `s` is a semi-infinite segment.
"""
@inline function issemiinf(s::Segment)
    return xor(any(isinf,s[1]),any(isinf,s[2]))
end

"""
$(TYPEDSIGNATURES)

Returns the index of the coordinate for which the first and/or last vertex of the vortex filament `vf` has an infinite value if it has one.
"""
function infdir(vf::VortexFilament)
    if isinf(vf) || issemiinf(vf)
        segdirs = infdir.(vf.segments)
        uniquesegdirs = unique(filter!(dir->dir≠nothing,segdirs))
        length(uniquesegdirs) > 1 && throw(ErrorException("Inf values in different coordinates"))
        return uniquesegdirs[1]
    else
        return nothing
    end
end

"""
$(TYPEDSIGNATURES)

Returns the index of the coordinate for which the first and/or second vertex of the segment `s` has an infinite value if it has one.
"""
function infdir(s::Segment)
    !all(length.(findall.(isinf,s)) .< 2) && throw(ErrorException("Inf values in different coordinates"))
    if issemiinf(s)
        Inf in abs.(s[1]) ? infidx = 1 : infidx = 2
        return findnext(isinf,s[infidx],1)
    elseif isinf(s)
        dir1 = findnext(isinf,s[1],1)
        dir2 = findnext(isinf,s[2],1)
        dir1 == dir2 ? dir = dir1 : throw(ErrorException("Inf values in different coordinates"))
        sign(s[1][dir]) == sign(s[2][dir]) && throw(ErrorException("Infinite segment has non-opposite infinite vertices"))
        return dir
    else
        return nothing
    end
end

"""
$(TYPEDSIGNATURES)

Computes the induced velocity by one segment `s` of a unit-strength vortex filament on the evaluation point `xeval`. In case both vertices of the segment have an infinite coordinate (which have to be in the same direction), then the non-infinite coordinates are averaged over the two vertices to calculate the distance from xeval to the segment.
"""
function inducevelocity(s::Segment,xeval)
    if s[1] == s[2]
        return zeros(3)
    elseif issemiinf(s)
        #=
        In this case we calculate the induced velocity by splitting the semi-infinite segment into two segments using the point `v`, which is the closest point to xeval on the axis of the semi-infinite segment. At the end we sum the velocities induced by the two segments using two multipliers `K1` and `K2` to determine the sign of the velocities.

        The first segment `vfiniteseg` is the segment from the finite end of `s` to `v`. Because the sign of the velocity induced by this segment flips if `v` lies outside of `s` (because then this velocity has to be subtracted) and if the order of the two vertices of `s` is flipped, we multiply it by `K1` and `K2` when computing `vnet`.

        The second segment `vseminfseg` is the segment from the infinite end of `s` to `v`. Because the sign of the velocity induced by this segment flips if the order of the two vertices of `s` is flipped, we multiply it by `K1` when computing `vnet`.
        =#

        infidx = findnext(v->in(Inf,abs.(v)),s,1)
        infidx == 1 ? (finidx = 2; K1 = -1) : (finidx = 1; K1 = 1)

        dir = infdir(s)
        edir = zeros(3)
        edir[dir] = sign(s[infidx][dir])

        v,t = _closestpointonfinitesegmentaxis(Segment(s[finidx],s[finidx]+edir),xeval)
        t < 0 ? K2 = -1 : K2 = 1 # if t < 0, v does not lie on the segment
        vfiniteseg = inducevelocity(Segment(v,s[finidx]),xeval)

        w = xeval - v

        if norm(w) < eps
            vnet = zeros(3)
        else
            vseminfseg = -1/(4π)*cross(w,edir)/dot(w,w)
            vnet = K1*K2*vfiniteseg + K1*vseminfseg
        end
        return vnet

    elseif isinf(s)
        dir = infdir(s)
        edir = zeros(3)
        edir[dir] = sign(s[2][dir])

        centerpoint = (s[1]+s[2])/2
        centerpoint[dir] = 0.0

        v,t = _closestpointonfinitesegmentaxis(Segment(centerpoint,centerpoint+edir),xeval)

        w = xeval - v

        if norm(w) < eps
            vinfseg = zeros(3)
        else
            vinfseg = -1/(2π)*cross(w,edir)/dot(w,w)
        end
        return vinfseg
    else
        r0 = s[2] - s[1]
        r1 = xeval - s[1]
        r2 = xeval - s[2]
        r1norm = norm(r1)
        r2norm = norm(r2)

        r1r2 = cross(r1,r2)
        r1r2sq = dot(r1r2,r1r2)

        if r1norm < eps || r2norm < eps || r1r2sq < eps
            vseg = zeros(3)
        else
            r0r1 = dot(r0,r1)
            r0r2 = dot(r0,r2)
            K = 1.0/(4.0*π*r1r2sq)*(r0r1/r1norm - r0r2/r2norm)
            vseg = K*r1r2
        end

        return vseg
    end
end

"""
$(TYPEDSIGNATURES)

Computes the induced velocity by the vortex filament `vf` on the evaluation point `xeval`.
"""
function inducevelocity(vf::VortexFilament,xeval)
    vel = sum(inducevelocity.(vf.segments,Ref(xeval)))
    vel *= vf.Γ
end

@inline function _closestpointonfinitesegmentaxis(s::Segment,p)
    a = s[1]
    b = s[2]
    ab = b-a
    ap = p-a
    t = dot(ap,ab)/dot(ab,ab)
    return a + t*ab, t
end

lims(dir) = (:xlims, :ylims, :zlims)[dir]

@recipe function f(vf::VortexFilament)
    segments = vf.segments
    if isinf(vf) || issemiinf(vf)
        segments = deepcopy(vf.segments)
        for s in segments
            if isinf(s) || issemiinf(s)
                inftydir = infdir(s)
                noninftydir = filter!(x->x≠inftydir,[1,2,3])
                inflims = get(plotattributes,lims(inftydir),:none)
                inflims == :none && throw(ArgumentError("The filament has a vertex at infinity, please provide axis limits"))
                if isinf(s)
                    centerpoint = 0.5*(s[1] + s[2])
                    s[1][noninftydir] .= centerpoint[noninftydir]
                    s[2][noninftydir] .= centerpoint[noninftydir]
                end
                if isinf(s[1][inftydir])
                    sign(s[1][inftydir]) == -1 ? limitidx = 1 : limitidx = 2
                    s[1][inftydir] = inflims[limitidx]
                    s[1][noninftydir] .= s[2][noninftydir]
                end
                if isinf(s[2][inftydir])
                    sign(s[2][inftydir]) == -1 ? limitidx = 1 : limitidx = 2
                    s[2][inftydir] = inflims[limitidx]
                    s[2][noninftydir] .= s[1][noninftydir]
                end
            end
        end
    end
    vertices = vcat((s->s[1]).(segments),[segments[end][end]])
    return (v->v[1]).(vertices), (v->v[2]).(vertices), (v->v[3]).(vertices)
end

Base.:(==)(a::VortexFilament, b::VortexFilament) = isequal(a.Γ, b.Γ) && isequal(a.vertices, b.vertices) &&  isequal(a.segments, b.segments) && isequal(a.freeidx, b.freeidx) && isequal(a.boundidx, b.boundidx)

function Base.show(io::IO, vf::VortexFilament)
    isinf(vf) ? type = "Infinite" : issemiinf(vf) ? type = "Semi-infinite" : type = "Finite"
    nv = length(vf.vertices)
    ns = length(vf.segments)

    println(io, "$type vortex filament with $((nv == 1) ? "1 vertex" : "$nv vertices"), $((ns == 1) ? "1 segment" : "$ns segments"), and strength Γ = $(vf.Γ)")
end


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
using DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear,
       prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear

probArr = Vector{ODEProblem}(undef, 2)
bigprobArr = Vector{ODEProblem}(undef, 2)

probArr[1] = prob_ode_linear
probArr[2] = prob_ode_2Dlinear
bigprobArr[1] = prob_ode_bigfloatlinear
bigprobArr[2] = prob_ode_bigfloat2Dlinear
setprecision(400)
Random.seed!(100)
## Convergence Testing
println("Convergence Test on Linear")
dts = 1 .//2 .^(8:-1:4)
testTol = 0.3
superduperbool = Vector{Bool}(undef, 2)

for i = 1:2 # 1 = num, 2 = ExplicitRK
  global dts
  if i>1
    prob = probArr[2]
    bigprob = bigprobArr[2]
  else
    prob = probArr[1]
    bigprob = bigprobArr[1]
  end

  # Order 2

  tabalg = ExplicitRK(tableau=constructHeun())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  tabalg = ExplicitRK(tableau=constructRalston())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  tabalg = ExplicitRK(tableau=constructSSPRK22())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-2) < testTol

  # Order 3

  tabalg = ExplicitRK(tableau=constructBogakiShampine3())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-3) < testTol

  tabalg = ExplicitRK(tableau=constructSSPRK33())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-3) < testTol

  tabalg = ExplicitRK(tableau=constructSSPRK43())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-3) < testTol

  # Order 4

  tabalg = ExplicitRK(tableau=constructRKF4())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-4) < testTol

  tabalg = ExplicitRK(tableau=constructSSPRK104())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-4) < testTol

  # Order 5

  dts = 1 .//2 .^(7:-1:4)
  tabalg = ExplicitRK(tableau=constructRKF5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1 .//2 .^(7:-1:4)
  tabalg = ExplicitRK(tableau=constructDormandPrince())
  sim = test_convergence(dts,prob,tabalg)
  sim2 = test_convergence(dts,prob,DP5())
  @test (abs(sim.𝒪est[:l∞]-5) < testTol && (maximum(sim[end][end]-sim2[end][end]) < 1e-10))

  tabalg = ExplicitRK(tableau=constructCashKarp())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1 .//2 .^(7:-1:4)
  tabalg = ExplicitRK(tableau=constructRungeFirst5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructCassity5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructLawson5())
  sim = test_convergence(dts,prob,tabalg) #10
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructLutherKonen5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructLutherKonen52())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructLutherKonen53())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructPapakostasPapaGeorgiou5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructPapakostasPapaGeorgiou52())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1 .//2 .^(6:-1:4)
  tabalg = ExplicitRK(tableau=constructTsitouras5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  dts = 1 .//2 .^(6:-1:4)
  tabalg = ExplicitRK(tableau=constructBogakiShampine5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  tabalg = ExplicitRK(tableau=constructSharpSmart5())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5) < testTol

  # Order 6

  dts = 1 .//2 .^(6:-1:4)
  tabalg = ExplicitRK(tableau=constructButcher6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructButcher62())
  sim = test_convergence(dts,prob,tabalg) #20
  @test abs(sim.𝒪est[:l∞]-6) < testTol # Less stringent

  dts = 1 .//2 .^(6:-1:4)
  tabalg = ExplicitRK(tableau=constructButcher63())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructDormandPrince6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol+0.1 # Better on linear

  tabalg = ExplicitRK(tableau=constructSharpVerner6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructVerner916())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructVerner9162())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructVernerRobust6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructVernerEfficient6(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6.6) < testTol

  tabalg = ExplicitRK(tableau=constructPapakostas6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructLawson6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1 .//2 .^(3:-1:1)
  tabalg = ExplicitRK(tableau=constructTsitourasPapakostas6())
  sim = test_convergence(dts,prob,tabalg) #30
  @test abs(sim.𝒪est[:l∞]-6.7) < testTol # Better on linear

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructDormandLockyerMcCorriganPrince6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructTanakaKasugaYamashitaYazaki6D())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructTanakaKasugaYamashitaYazaki6C())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructTanakaKasugaYamashitaYazaki6B())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructTanakaKasugaYamashitaYazaki6A())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructMikkawyEisa())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6.53) < testTol # Odd behavior

  tabalg = ExplicitRK(tableau=constructChummund6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructChummund62())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructHuta6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-5.5) < testTol # Low convergence, error noted in Stone notes

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructHuta62())
  sim = test_convergence(dts,prob,tabalg)#40
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructVerner6())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6.7) < testTol # Better on linear

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructDverk())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  tabalg = ExplicitRK(tableau=constructClassicVerner6())
  sim = test_convergence(dts,prob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6) < testTol

  # Order 7

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructButcher7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1 .//2 .^(5:-1:2)
  tabalg = ExplicitRK(tableau=constructClassicVerner7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  tabalg = ExplicitRK(tableau=constructVernerRobust7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructEnrightVerner7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7.15) < testTol # Better on linear

  tabalg = ExplicitRK(tableau=constructTanakaYamashitaStable7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7.3) < testTol

  tabalg = ExplicitRK(tableau=constructTanakaYamashitaEfficient7(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1 .//2 .^(8:-1:3)
  tabalg = ExplicitRK(tableau=constructSharpSmart7(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg) #50
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  dts = 1 .//2 .^(3:-1:1)
  tabalg = ExplicitRK(tableau=constructSharpVerner7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-6.5) < testTol # Coefficients aren't accurate enough, drop off error

  tabalg = ExplicitRK(tableau=constructVernerEfficient7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  tabalg = ExplicitRK(tableau=constructVerner7())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-7) < testTol

  # Order 8
  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructClassicVerner8())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructVerner8(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructCooperVerner8())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol #Coefficients not accurate enough

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructCooperVerner82())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol #Coefficients not accurate enough

  tabalg = ExplicitRK(tableau=constructTsitourasPapakostas8(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructdverk78(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructEnrightVerner8(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructCurtis8())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  dts = 1 .//2 .^(4:-1:1)
  tabalg = ExplicitRK(tableau=constructRKF8(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8) < testTol

  tabalg = ExplicitRK(tableau=constructDormandPrince8(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8.4) < testTol

  dts = 1 .//2 .^(3:-1:1)
  tabalg = ExplicitRK(tableau=constructDormandPrince8_64bit(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-8.4) < testTol

  # Order 9

  tabalg = ExplicitRK(tableau=constructVernerRobust9(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-9) < testTol

  tabalg = ExplicitRK(tableau=constructVernerEfficient9(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-9) < testTol

  dts = 1 .//2 .^(3:-1:1)
  tabalg = ExplicitRK(tableau=constructSharp9(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-9) < testTol #Only works to Float64 precision

  dts = 1 .//2 .^(2:-1:1)
  tabalg = ExplicitRK(tableau=constructTsitouras9(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10.5) < testTol #Only works to Float64

  dts = 1 .//2 .^(3:-1:1)
  tabalg = ExplicitRK(tableau=constructTsitouras92(BigFloat))
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-9) < testTol  #Only works to Float64

  ## Order 10

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructCurtis10())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  tabalg = ExplicitRK(tableau=constructOno10())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  dts = 1 .//2 .^(5:-1:1)
  tabalg = ExplicitRK(tableau=constructFeagin10Tableau())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  tabalg = ExplicitRK(tableau=constructCurtis10())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10) < testTol

  dts = 1 .//2 .^(6:-1:1)
  tabalg = ExplicitRK(tableau=constructBaker10())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10.8) < testTol

  tabalg = ExplicitRK(tableau=constructHairer10())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-10.8) < testTol

  ## Order 12

  dts = 1 .//2 .^(6:-1:1)
  tabalg = ExplicitRK(tableau=constructFeagin12Tableau())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-12.6) < testTol

  tabalg = ExplicitRK(tableau=constructOno12())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-11.6) < testTol

  ## Order 14

  dts = 1 .//2 .^(6:-1:1)
  tabalg = ExplicitRK(tableau=constructFeagin14Tableau())
  sim = test_convergence(dts,bigprob,tabalg)
  @test abs(sim.𝒪est[:l∞]-15.5) < testTol
end

ext: jl


text:
#################################################################################################
# Functions needed for the manipulation of datacell labels
#
# Desired functionality:
#     - adding labels to existing DataCell
#     - removing labels from existing DataCell
#     - changing some/all labels from existing DataCell
#     - remap labels to something else : e.g. from "apple" to "bear" etc using regexpz
#     - extract column(s) from data and transform them into labels
#
# Observation: All functions return new objects (they do not modify the existing ones)
#################################################################################################

# Remove labels from DataCells
# - Can be done by creating a new datacell(old_datacell.x) or datacell(old_datacell.x, old_datacell.y[:,...])

# Change labels from existing datacell
# - Extract current labels, modify accordingly and create new datacell re-using the data and new labels.

"""
	addlabels(x, labels)

Creates a new datacell by adding `labels<:AbstractArray` to the existing labels
of data cell `x`.
"""
addlabels(x::T where T<:CellData, labels::S where S<:AbstractArray) = datacell(getx!(x), dcat(gety!(x), labels), name=x.tinfo)

"""
	unlabel(data)

Removes labels (if any) from the data cell and returns a new datacell with the
data contents only.
"""
unlabel(x::T where T<:CellData) = datacell(getx(x); name = x.tinfo)

"""
	labelize(x, f [,idx ,remove])
	labelize(x, idx [,remove])
	labelize(x, v [,remove])

Adds labels (also known as 'targets') to `x::DataCell`. If `x` is not provided, returns a fixed `Function cell` that when piped a DataCell into,
it will add labels/targets and return it.

# Arguments
  * `f` is a `targets`-related function (e.g. `targets(f,...)` ) that it is applied to the existing labels of `x` if `idx` is not specified or to the
  variables of `x` indicated by `idx`
  * `idx` specifies variables in `x`; can be anything that can be used as a variable index in a `DataCell`
  * `v` a vector that will become the new targets of `x`
  * `remove` defaults to `true` and specifies whether any existing labels/targets are to be kept or, if `idx` is present, whether to remove the
  variables from `x` from which the new labels/targets were obtained.
"""
labelize(f::T where T<:Function, remove::Bool=true) = FunctionCell(labelize, (f,remove), "Data labeler: f=$f, remove=$remove")
labelize(v::T where T<:AbstractArray, remove::Bool=true) = FunctionCell(labelize, (v,remove), "Data labeler: preloaded targets, remove=$remove")
labelize(idx::T where T, remove::Bool=true) = FunctionCell(labelize, (idx,remove), "Data labeler: idx=$idx remove=$remove")
labelize(f::T where T<:Function, idx::S where S, remove::Bool=true) = FunctionCell(labelize, (f,idx,remove), "Data labeler: f=$f idx=$idx remove=$remove")

labelize(x::T where T<:CellDataU, f::Function, remove::Bool=true) =
	error("[labelize] Targets or the indices of variables from which to create targets required.")

labelize(x::T where T<:CellData, f::Function, remove::Bool=true) = begin
	if remove
		datacell(getx(x), targets(f, gety(x)), name=x.tinfo) 			# replace labels
	else
		datacell(getx(x), dcat(gety(x), targets(f, gety(x))), name=x.tinfo)	# add to existing labels
	end
end

labelize(x::T where T<:CellData, v::S where S<:AbstractArray, remove::Bool=true) = begin
	if remove
		return datacell( getx(x), getobs(v), name=x.tinfo) 			# replace labels
	else
		return datacell( getx(x), dcat(gety(x), getobs(v)), name=x.tinfo ) 	# add to existing labels
	end
end

labelize(x::T where T<:CellData, idx::S where S, remove::Bool=true) = labelize(x, identity, idx, remove)

labelize(x::T where T<:CellData, f::Function, idx::S where S, remove::Bool=true) = begin
	labels = targets(f, getx!(varsubset(x,idx)))
	if remove
		return datacell( getobs(_variable_(x, setdiff(1:nvars(x),idx))), labels, name=x.tinfo)
	else
		return datacell( getx(x), labels, name=x.tinfo )
	end
end


text:
"""
	FunctionCell{T,S,U,V,W}(x::T, y::S, f::U, fargs::V<:Tuple, fkwargs::W<:Tuple, layer::Int, tinfo::String, oinfo::String)

Constructs the basic object for processing data, the `FunctionCell`. The constructor is parametric so new varieties of
`FunctionCells` can be aliased to cover new functionality.

# Main `FunctionCell` constructors
  * `FunctionCell(f::Function, fargs::Tuple, tinfo=string(f); fkwargs...)` creates a `fixed` FunctionCell or transform.
  * `FunctionCell(f::Function, fargs::Tuple, mprops::ModelProperties, tinfo=string(f); fkwargs...)` creates an `untrained` FunctionCell.
  * `FunctionCell(f::Function, m::Model, tinfo=string(f); fkwargs...)` creates an `trained` FunctionCell

where, `m` is the model, an object of type Model{T} (which is in itself a simple wrapper).

  In the above constructors, `f` is an overloaded function that, depending on circumstances, is applied on the input data piped
or called by the cell, `fargs` are `f`'s arguments, `tinfo` is the title and `fkwargs` are the keyword arguments of `f`. For
fixed function cells, `f` is the the transform applied on the input; for untrained function cells, `f` is the training algorithm which
generates a model and some arbitrary properties of the model, such as priors, expected input/output dimensions etc.
generated in training. The body of `f` for untrained cells should return a trained cell where `f` is the model execution function,
`m` is the model and `mprops` the model properties.
"""
#################
# FunctionCells #
#################
struct FunctionCell{T,S,U, V<:Tuple, W<:Tuple} <: AbstractCell{T,S,U}
	x::T                                                                                # Model
    	y::S                                                                                # Model properties
    	f::U                                                                                # Function for transforming/training/predicting
	fargs::V										# Function arguments
	fkwargs::W										# Function keyword arguments
	layer::Int										# Layer
    	tinfo::String										# Title or name
    	oinfo::String                                                                        	# Other information
end

"""
	ModelProperties([idim ,odim, labels, other])

Basic type used to describe the properties of a model.

# Arguments
  * `idim::Int` is the required input data dimension (0 stands for any dimension)
  * `odim::Int` is the expected dimension of the output of the model (0 stands for any dimension)
  * `labels::LabelEncoding` the encoding of the labels
  * `other::Dict` are other properties of the data i.e. priors, custom execution arguments
"""
struct ModelProperties{T<:LabelEncoding, S<:Dict}
	idim::Int
	odim::Int
	labels::T
	other::S
end

ModelProperties(idim::Int=0, odim::Int=0, enclabels::LabelEncoding=LabelEnc.NativeLabels([],Val{0})) =
	ModelProperties(idim, odim, enclabels, Dict())

labelencn(labels::T where T<:AbstractVector{S}) where S = begin
	ulabels::Vector{S} = sort(unique(labels))
 	return LabelEnc.NativeLabels(ulabels)
end

ModelProperties(idim::Int, odim::Int, labels::Void, other::Dict=Dict()) =
	ModelProperties(idim, odim, LabelEnc.NativeLabels([],Val{0}), other)

ModelProperties(idim::Int, odim::Int, labels::T where T<:AbstractVector, other::Dict=Dict()) =	# Known dimensions and labels
	ModelProperties(idim, odim, labelencn(labels), other)

ModelProperties(idim::Int, odim::Int, labels::T where T<:AbstractMatrix, other::Dict=Dict()) =	# Known dimensions and labels
	ModelProperties(idim, odim, labelenc(labels), other)

Base.show(io::IO, mprop::ModelProperties{T,S}) where {T,S} = begin
	println("Model properties")
	println("`- I/O dimensions: $(mprop.idim)/$(mprop.odim)")
	println("`- label encoding: $(typeof(mprop.labels))")
	print("`- other: $(mprop.other)")
end

"""
	Model(data [,properties])

Basic wrapper around model data and propertiesm, containing all the data needed for the execution of a model,
given an execution function of the form:
	`f(x, m)`
where `x` is the input data and `m::Model` the model.

# Arguments
  * `data` is the input data on which the model is applied
  * `properties::ModelProperties` are the properties of the model, containing parameters necessary for its correct operation
"""
struct Model{T, S<:ModelProperties}
	data::T
	properties::S
end

Model() = Model(nothing, ModelProperties()) 							# Empty model

Model(data) = Model(data, ModelProperties())							# Model without properties

Model(properties::S where S<:ModelProperties) = Model(nothing, properties)			# Model without data

Model(model::T where T<:Model, properties::S where S<:ModelProperties) =			# Grabdata from a model and replace properties
	Model(model.data, properties)

Base.show(io::IO, m::Model) = begin
	println("Model")
	println("`- data: $(typeof(m.data))")
	print("`- properties: Int/Int/$(typeof(m.properties.labels))/$(typeof(m.properties.other))")
end

# Constructor for fixed FunctionCells e.g. simple transforms
FunctionCell(f::U where U<:Function, fargs::V where V<: Tuple, tinfo::String = string(f); fkwargs...) =
	FunctionCell(nothing, nothing, f, fargs, ntuple(i->fkwargs[i], Val{length(fkwargs)}), 0, tinfo, oinfoglobal)

# Constructor for untrained FunctionCells e.g. untrained classifiers or trainable transforms
FunctionCell(f::U where U<:Function, fargs::V where V<:Tuple, mprops::ModelProperties, tinfo::String = string(f); fkwargs...) =
	FunctionCell(nothing, mprops, f, fargs, ntuple(i->fkwargs[i], Val{length(fkwargs)}), 0, tinfo, oinfoglobal)

# Constructor for trained FunctionCells (no arguments, keyword arguments, the model should fully describe the operation)
FunctionCell(f::U where U<:Function, m::T where T<:Model, tinfo::String = string(f)) =
	FunctionCell(m, nothing, f, (), (),  0, tinfo, oinfoglobal)

########################
# FunctionCell aliases #
########################
const CellFunF{U<:Function} = 			 FunctionCell{Void, Void, U}			# Fixed function cell: simple transforms
const CellFunU{S<:ModelProperties, U<:Function}= FunctionCell{Void, S, U}			# Untrained function cell: untrained classifier or transform
const CellFunT{T<:Model, S<:Void, U<:Function} = FunctionCell{T, S, U}				# Trained function cell: trained classifier, transform
const CellFun{T,S, U<:Function} = 		 FunctionCell{T, S, U}

##############################################################################################################################
# Operators [Function cells]																								 #
##############################################################################################################################

# Pipe operators for function cells (a different execution method for each FunctionCell variety: fixed, untrained and trained)
|>(x::T where T, c::S where S<:CellFunF) = getf!(c)(x, c.fargs...; c.fkwargs...)
|>(x::T where T, c::S where S<:CellFunU) = getf!(c)(x, c.fargs...; c.fkwargs...)
|>(x::T where T, c::S where S<:CellFunT) = begin
	# Check input dimension (if input dimension is 0, skip the check, input dimension can vary)
	getx!(c).properties.idim != 0 && getx!(c).properties.idim !=nvars(x) &&
		error("Expected $(getx!(c).properties.idim) input variables(s), got $(nvars(x)).")
	return getf!(c)(x, getx!(c))
end

"""
	functioncell([x,y,z])

Calls the `FunctionCell` constructor for one, two or three arguments. If no arguments are provided,
a `FunctionCell` that expects as input a collection with contents compatible with the arguments of
the `FunctionCell` constructor is returned.
"""
functioncell() = FunctionCell((args)->FunctionCell(args...), (), "FunctionCell creator")
functioncell(x::T) where T = FunctionCell(x)					# calls to convert methods
functioncell(x::T, y::S) where {T,S}= FunctionCell(x, y)			# Calls any 2 argument parametric constructor
functioncell(x::T, y::S, z::U) where {T,S,U}= FunctionCell(x, y, z)		# Calls any 3 argument parametric constructor


text:
function test1()
  filename = "test1.asc"
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)
  @static if Sys.iswindows()
    (circuitdir,_filename) = splitdir(circuitpath(sim))
    @test _filename == filename
    (logdir,logfilename) = splitdir(logpath(sim))
    @test logdir == circuitdir
    @test executablepath(sim) == ""
  end
  @test parameternames(sim) == ("Vin","load")
  @test measurementnames(sim) == ("Current",)
  @test stepnames(sim) == ()
  @test parametervalues(sim) == [5.00, 2.00]
  @test measurementvalues(sim) == [2.50]
  show(IOBuffer(),sim)
  @test stepvalues(sim) == ()
  @test sim["Vin"] == 5
  @test sim["load"] == 2
  @test sim["Current"] == 2.5
  @test keys(sim) == ["Vin","load","Current"]
  @test values(sim) == [5.0,2.0,2.5]
  for key in keys(sim)
    @test haskey(sim,key)
    @test sim[key] == get(sim,key,NaN)
  end
  @test_throws KeyError sim["not a valid key"]
  @test_throws KeyError sim["not a valid key"] = 1.0
  @test isnan(get(sim,"not a valid key",NaN))
  @test eltype(sim) == Float64
  @test length(sim) == 3
  sim["Vin"] = 1.0
  @test sim["Vin"] == 1.0
  sim["Vin"] = 5.0
end
test1()


text:
using Documenter, InfrastructureSystems
import DataStructures: OrderedDict
using Literate

if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

pages = OrderedDict(
        "Welcome Page" => "index.md",
        "Developer Guide" => Any[
            "dev_guide/components_and_container.md",
            "dev_guide/auto_generation.md",
            "dev_guide/time_series.md",
            "dev_guide/recorder.md",
            "dev_guide/tests.md",
            "dev_guide/logging.md",
        ],
        "Style Guide" => "style.md",
        "API" => "InfrastructureSystems.md"
)

# This code performs the automated addition of Literate - Generated Markdowns. The desired
# section name should be the name of the file for instance network_matrices.jl -> Network Matrices
# This code is generic to all SIIP documentation
julia_file_filter = x -> occursin(".jl", x)
folders = Dict(
    "Developer Guide" => filter(julia_file_filter, readdir("docs/src/dev_guide")),
)

for (section, folder) in folders
    for file in folder
        section_folder_name = lowercase(replace(section, " " => "_"))
        outputdir = joinpath(pwd(), "docs", "src", "$section_folder_name")
        inputfile = joinpath("$section_folder_name", "$file")
        outputfile = string("generated_", replace("$file", ".jl" => ""))
        Literate.markdown(joinpath(pwd(), "docs", "src", inputfile),
                          outputdir;
                          name = outputfile,
                          credit = false,
                          execute=true)
        subsection = titlecase(replace(split(file, ".")[1], "_" => " "))
        push!(pages[section], ("$subsection" =>  joinpath("$section_folder_name", "$(outputfile).md")))
    end
end

makedocs(
    modules = [InfrastructureSystems],
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS"),),
    sitename = "InfrastructureSystems.jl",
    pages = Any[p for p in pages]
)

deploydocs(
    repo = "github.com/NREL-SIIP/InfrastructureSystems.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
    push_preview = true,
)


text:
module TableIO

export read_table, write_table!, read_sql, list_tables

using TableIOInterface
using Tables, Requires, Suppressor
using DataFrames # required for multiple file types, therefore currently not optional

# specify if a reader accepts an io buffer as input or if creation of a temp file is required
supports_io_input(::TableIOInterface.AbstractFormat) = false

supports_io_input(::TableIOInterface.CSVFormat) = true
supports_io_input(::TableIOInterface.JSONFormat) = true
supports_io_input(::TableIOInterface.ArrowFormat) = true

# definition of the required packages for the specific formats
# if a format requires multiple packages, define them as a list
const PACKAGE_REQUIREMENTS = Dict{DataType, Union{Symbol, Vector{Symbol}}}(
    TableIOInterface.CSVFormat => :CSV,
    TableIOInterface.ZippedFormat => :ZipFile,
    TableIOInterface.JDFFormat => :JDF,
    TableIOInterface.ParquetFormat => :Parquet,
    TableIOInterface.ExcelFormat => :XLSX,
    TableIOInterface.SQLiteFormat => :SQLite,
    TableIOInterface.StataFormat => :StatFiles,
    TableIOInterface.SPSSFormat => :StatFiles,
    TableIOInterface.SASFormat => :StatFiles,
    TableIOInterface.JSONFormat => :JSONTables,
    TableIOInterface.ArrowFormat => :Arrow,
    TableIOInterface.PostgresFormat => [:LibPQ, :CSV],
    TableIOInterface.HDF5Format => :Pandas,
    TableIOInterface.JLD2Format => :JLD2,
)

## Dispatching on file extensions

"""
    read_table(filename:: AbstractString; kwargs...)

`filename`: path and filename of the input file
`kwargs...`: keyword arguments passed to the underlying file reading function (e.g. `CSV.File`)

Returns a Tables.jl interface compatible object.

Example:

    df = DataFrame(read_table("my_data.csv"); copycols=false)

"""
function read_table(filename:: AbstractString, args...; kwargs...)
    data_type = get_file_type(filename)
    try
        # to speed up the standard case (format specific package is imported), this is tried first
        return read_table(data_type, filename, args...; kwargs...)
    catch ex
        if ex isa MethodError
            # import format specific package and invoke latest version of the function to avoid world age issues
            _import_package(data_type)
            return Base.invokelatest(read_table, data_type, filename, args...; kwargs...)
        else
            rethrow()
        end
    end

end

"""
    read_table(file_picker:: Dict, args...; kwargs...)

Reading tabular data from a PlutoUI.jl FilePicker.

Usage (in a Pluto.jl notebook):

    using PlutoUI, TableIO, DataFrames
    using XLSX # import the packages required for the uploaded file formats
    @bind f PlutoUI.FilePicker()
    df = DataFrame(read_table(f); copycols=false)

"""
function read_table(file_picker:: Dict, args...; kwargs...)
    filename, data = _get_file_picker_data(file_picker)
    data_type = get_file_type(filename)
    data_buffer = IOBuffer(data)

    _import_package(data_type)

    if supports_io_input(data_type)
        data_object = data_buffer # if it is supported by the corresponding package, creation of a temporary file is avoided and the IOBuffer is used directly
    else
        tmp_file = joinpath(mktempdir(), filename)
        write(tmp_file, data_buffer)
        data_object = tmp_file
    end

    try
        # to speed up the standard case (format specific package is imported), this is tried first
        read_table(data_type, data_object, args...; kwargs...)
    catch ex
        if ex isa MethodError
            # import format specific package and invoke latest version of the function to avoid world age issues
            _import_package(data_type)
            return Base.invokelatest(read_table, data_type, data_object, args...; kwargs...)
        else
            rethrow()
        end
    end
end

"""
    list_tables(filename:: AbstractString)

Returns a list of all tables inside a file.
"""
function list_tables(filename:: AbstractString)
    data_type = get_file_type(filename)
    TableIOInterface.multiple_tables(data_type) || error("The data type $data_type does not support multiple tables per file.")
    try
        # to speed up the standard case (format specific package is imported), this is tried first
        return list_tables(data_type, filename)
    catch ex
        if ex isa MethodError
            # import format specific package and invoke latest version of the function to avoid world age issues
            _import_package(data_type)
            return Base.invokelatest(list_tables, data_type, filename)
        else
            rethrow()
        end
    end

end

"""
    write_table!(filename:: AbstractString, table; kwargs...):: AbstractString

`filename`: path and filename of the output file
`table`: a Tables.jl compatible object (e.g. a DataFrame) for storage
`kwargs...`: keyword arguments passed to the underlying file writing function (e.g. `CSV.write`)

Example:

    write_table!("my_output.csv", df)

"""
function write_table!(filename:: AbstractString, table, args...; kwargs...)
    data_type = get_file_type(filename)
    try
        # to speed up the standard case (format specific package is imported), this is tried first
        write_table!(data_type, filename, table, args...; kwargs...)
    catch ex
        if ex isa MethodError
            # import format specific package and invoke latest version of the function to avoid world age issues
            _import_package(data_type)
            return Base.invokelatest(write_table!, data_type, filename, table, args...; kwargs...)
        else
            rethrow()
        end
    end
    nothing
end

"""
    read_sql(db, sql:: AbstractString)

Returns the result of the SQL query as a Tables.jl compatible object.
"""
function read_sql end

include("julia.jl")

## conditional dependencies

function __init__()
    @require CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" begin
        include("csv.jl")
        @require LibPQ = "194296ae-ab2e-5f79-8cd4-7183a0a5a0d1" include("postgresql.jl")
    end
    @require ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" include("zip.jl")
    @require JDF = "babc3d20-cd49-4f60-a736-a8f9c08892d3" include("jdf.jl")
    @require Parquet = "626c502c-15b0-58ad-a749-f091afb673ae" include("parquet.jl")
    @require XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0" include("xlsx.jl")
    @require StatFiles = "1463e38c-9381-5320-bcd4-4134955f093a" include("stat_files.jl")
    @require SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" include("sqlite.jl")
    @require JSONTables = "b9914132-a727-11e9-1322-f18e41205b0b" include("json.jl")
    @require Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" include("arrow.jl")
    @require Pandas = "eadc2687-ae89-51f9-a5d9-86b5a6373a9c" include("pandas.jl")
    @require JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" include("jld2.jl")
end

## Utilities

get_package_requirements(::T) where {T <: TableIOInterface.AbstractFormat} = PACKAGE_REQUIREMENTS[T]
get_package_requirements(filename:: AbstractString) = get_package_requirements(get_file_type(filename))

function _import_package(t::T) where {T <: TableIOInterface.AbstractFormat}
    pkg_name = get_package_requirements(t)
    _import_package(pkg_name)
end

function _import_package(pkg_names:: Vector{Symbol})
    for pkg_name ∈ pkg_names
        _import_package(pkg_name)
    end
end

function _import_package(pkg_name:: Symbol)
    # A warning is raised if a package is imported which is not a dependency of TableIO. This warning is suppressed.
    # If the package is not installed, an error message is raised.
    try
        @suppress @eval import $pkg_name
    catch ex
        # If the package is not installed, the error message is swallowed by @suppress, but the warning message for a missing TableIO dependeny is raised.
        # To get back the more helpful error message for a not installed package, it is regenerated below.
        if ex isa ArgumentError
            throw(ArgumentError("""
                ERROR: ArgumentError: Package $pkg_name not found in current path:
                - Run `import Pkg; Pkg.add("$pkg_name")` to install the $pkg_name package.
                """))
        else
            rethrow()
        end
    end

    # note that it is required to use Base.invokelatest for calling any functionality depending on the imported package, unless one returns to global scope before.
end

_checktable(table) = Tables.istable(typeof(table)) || error("table has no Tables.jl compatible interface")

# poor man's approach to prevent SQL injections / garbage inputs
_checktablename(tablename) = match(r"^[a-zA-Z0-9_]*$", tablename) === nothing && error("tablename must only contain alphanumeric characters and underscores")

function _get_file_picker_data(file_picker:: Dict)
    data = file_picker["data"]:: Vector{UInt8} # brings back type stability
    length(data) == 0 && error("no file selected yet")
    filename = file_picker["name"]:: String
    return filename, data
end

end


text:
using Base.Test, EPInference, ExpFamily

# Proposing samples from a wide spherical gaussian.

srand(123)

C  = [1.0 0.5; 0.5 1.0]
P  = 2/3 * [2 -1; -1 2]
m  = [1.0;-0.5]

gNP     = GaussianNatParam(mean=m, cov=C)
gNPprop = ones(GaussianNatParam, 2)/20 # division INFLATES the variance

(samples,w) = impsampling(x->loglik(gNP,x), gNPprop, 10000)

@test isapprox(sum(samples[:,i]*w[i] for i in 1:length(w)), m, rtol=0.05)


text:
#
#  Cardinality assertions.
#

function ExpectQuery(base::Query, optional::Bool, plural::Bool)
    dom = domain(output(base))
    q = RecordQuery(istub(base), base)
    return (
        q >>
        Query(
            ExpectSig(optional, plural),
            Input(domain(output(q))),
            Output(dom, OutputMode(optional, plural))))
end

ExpectOneQuery(base::Query) = ExpectQuery(base, false, false)

ExpectAtMostOneQuery(base::Query) = ExpectQuery(base, true, false)

ExpectAtLeastOneQuery(base::Query) = ExpectQuery(base, false, true)

immutable ExpectSig <: AbstractPrimitive
    optional::Bool
    plural::Bool
end

function ev(sig::ExpectSig, dv::DataVector)
    icol = column(dv, 1)
    col = column(dv, 2)
    if !sig.optional
        expect_not_optional_impl!(icol, col)
    end
    if !sig.plural
        expect_not_plural_impl!(icol, col)
    end
    return Column{sig.optional,sig.plural}(col.offs, col.vals)
end

function expect_not_optional_impl!(icol::Column, col::Column)
    if !isoptional(col)
        return
    end
    icr = cursor(icol)
    cr = cursor(col)
    while !done(col, cr)
        next!(icol, icr)
        next!(col, cr)
        if length(cr) < 1
            error("expected at least one value: $(icr[1]) -> $cr")
        end
    end
end

function expect_not_plural_impl!(icol::Column, col::Column)
    if !isplural(col)
        return
    end
    icr = cursor(icol)
    cr = cursor(col)
    while !done(col, cr)
        next!(icol, icr)
        next!(col, cr)
        if length(cr) > 1
            error("expected at most one value: $(icr[1]) -> $cr")
        end
    end
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using Match
using Printf

import Base: read

function mapperstep!(console::Console{T}) where {T <: Mapper}
end

# Mapper 1

mutable struct Mapper1 <: Mapper
  shiftRegister::UInt8
  control::UInt8
  prgMode::UInt8
  chrMode::UInt8
  prgBank::UInt8
  chrBank0::UInt8
  chrBank1::UInt8
  prgOffsets::Vector{Int32}
  chrOffsets::Vector{Int32}

  function Mapper1(cartridge::Cartridge)
    m = new(0x10, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, zeros(Int32, 2), zeros(Int32, 2))
    @inbounds m.prgOffsets[2] = prgbankoffset(m, cartridge, Int32(-1))
    m
  end
end

function read(m::Mapper1, cartridge::Cartridge, address::UInt16)::UInt8
  if address < 0x2000
    bank = address ÷ 0x1000
    offset = address & 0x0FFF
    @inbounds cartridge.CHR[m.chrOffsets[bank + 1] + offset + 1]
  elseif address >= 0x8000
    address = address - 0x8000
    bank = address ÷ 0x4000
    offset = address & 0x3FFF
    @inbounds cartridge.PRG[m.prgOffsets[bank + 1] + offset + 1]
  elseif (address >= 0x6000)
    @inbounds cartridge.SRAM[address - 0x6000 + 1]
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper1 read at address: 0x%04X", address)
    0x00
  end
end

function write!(m::Mapper1, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address < 0x2000
    bank = address ÷ 0x1000
    offset = address & 0x0FFF
    @inbounds cartridge.CHR[m.chrOffsets[bank + 1] + offset + 1] = val
  elseif address >= 0x8000
    loadregister!(m, cartridge, address, val)
  elseif address >= 0x6000
    @inbounds cartridge.SRAM[address - 0x6000 + 1] = val
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper1 write at address: 0x%04X", address)
  end
end

function loadregister!(m::Mapper1, cartridge::Cartridge, address::UInt16, val::UInt8)
  if (val & 0x80) == 0x80
    m.shiftRegister = 0x10
    writecontrol!(m, cartridge, m.control | 0x0C)
  else
    complete = (m.shiftRegister & 0x01) == 0x01
    m.shiftRegister = m.shiftRegister >> 0x01
    m.shiftRegister = m.shiftRegister | ((val & 0x01) << 0x04)
    if complete
      writeregister!(m, cartridge, address, m.shiftRegister)
      m.shiftRegister = 0x10
    end
  end
end

function writeregister!(m::Mapper1, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address <= 0x9FFF
    writecontrol!(m, cartridge, val)
  elseif address <= 0xBFFF
    writechrbank0!(m, cartridge, val)
  elseif address <= 0xDFFF
    writechrbank1!(m, cartridge, val)
  else
    writeprgbank!(m, cartridge, val)
  end
end

# Control (internal, $8000-$9FFF)
function writecontrol!(m::Mapper1, cartridge::Cartridge, val::UInt8)
  m.control = val
  m.chrMode = (val >> 0x04) & 0x01
  m.prgMode = (val >> 0x02) & 0x03
  mirror = val & 0x03
  @match mirror begin
    0 => begin cartridge.mirror = UInt8(Single0::MirrorModes) end
    1 => begin cartridge.mirror = UInt8(Single1::MirrorModes) end
    2 => begin cartridge.mirror = UInt8(Vertical::MirrorModes) end
    3 => begin cartridge.mirror = UInt8(Horizontal::MirrorModes) end
    _ => throw(ErrorException("Unreachable code."))
  end
  updateoffsets!(m, cartridge)
end

# CHR bank 0 (internal, $A000-$BFFF)
function writechrbank0!(m::Mapper1, cartridge::Cartridge, val::UInt8)
  m.chrBank0 = val
  updateoffsets!(m, cartridge)
end

# CHR bank 1 (internal, $C000-$DFFF)
function writechrbank1!(m::Mapper1, cartridge::Cartridge, val::UInt8)
  m.chrBank1 = val
  updateoffsets!(m, cartridge)
end

# PRG bank (internal, $E000-$FFFF)
function writeprgbank!(m::Mapper1, cartridge::Cartridge, val::UInt8)
  m.prgBank = val & 0x0F
  updateoffsets!(m, cartridge)
end

function prgbankoffset(m::Mapper1, cartridge::Cartridge, index::Int32)::Int32
  if index >= 0x80
    index -= 0x100
  end
  index %= length(cartridge.PRG) >> 14
  offset = index * 0x4000
  if offset < 0
    offset += length(cartridge.PRG)
  end
  offset
end

function chrbankoffset(m::Mapper1, cartridge::Cartridge, index::Int32)::Int32
  if index >= 0x80
    index -= 0x100
  end
  index %= length(cartridge.CHR) >> 12
  offset = index * 0x1000
  if offset < 0
    offset += length(cartridge.CHR)
  end
  offset
end

# PRG ROM bank mode (0, 1: switch 32 KB at $8000, ignoring low bit of bank number;
#                    2: fix first bank at $8000 and switch 16 KB bank at $C000;
#                    3: fix last bank at $C000 and switch 16 KB bank at $8000)
# CHR ROM bank mode (0: switch 8 KB at a time; 1: switch two separate 4 KB banks)
function updateoffsets!(m::Mapper1, cartridge::Cartridge)
  @match m.prgMode begin
    0 => begin
      @inbounds m.prgOffsets[1] = prgbankoffset(m, cartridge::Cartridge, Int32(m.prgBank & 0xFE))
      @inbounds m.prgOffsets[2] = prgbankoffset(m, cartridge::Cartridge, Int32(m.prgBank | 0x01))
    end
    1 => begin
      @inbounds m.prgOffsets[1] = prgbankoffset(m, cartridge::Cartridge, Int32(m.prgBank & 0xFE))
      @inbounds m.prgOffsets[2] = prgbankoffset(m, cartridge::Cartridge, Int32(m.prgBank | 0x01))
    end
    2 => begin
      @inbounds m.prgOffsets[1] = 0
      @inbounds m.prgOffsets[2] = prgbankoffset(m, cartridge::Cartridge, Int32(m.prgBank))
    end
    3 => begin
      @inbounds m.prgOffsets[1] = prgbankoffset(m, cartridge::Cartridge, Int32(m.prgBank))
      @inbounds m.prgOffsets[2] = prgbankoffset(m, cartridge::Cartridge, Int32(-1))
    end
    _ => throw(ErrorException("Unreachable code."))
  end
  @match m.chrMode begin
    0 => begin
      @inbounds m.chrOffsets[1] = chrbankoffset(m, cartridge::Cartridge, Int32(m.chrBank0 & 0xFE))
      @inbounds m.chrOffsets[2] = chrbankoffset(m, cartridge::Cartridge, Int32(m.chrBank0 | 0x01))
    end
    1 => begin
      @inbounds m.chrOffsets[1] = chrbankoffset(m, cartridge::Cartridge, Int32(m.chrBank0))
      @inbounds m.chrOffsets[2] = chrbankoffset(m, cartridge::Cartridge, Int32(m.chrBank1))
    end
    _ => throw(ErrorException("Unreachable code."))
  end
end

# Mapper 2

mutable struct Mapper2 <: Mapper
  prgBanks::Int32
  prgBank1::Int32
  prgBank2::Int32

  function Mapper2(cartridge::Cartridge)
    m = new(0, 0, 0)
    m.prgBanks = length(cartridge.PRG) ÷ 0x4000
    m.prgBank2 = m.prgBanks - Int32(1)
    m
  end
end

function read(m::Mapper2, cartridge::Cartridge, address::UInt16)::UInt8
  if address < 0x2000
    cartridge.CHR[address + 1]
  elseif address >= 0xC000
    index = m.prgBank2 * Int32(0x4000) + Int32(address - 0xC000)
    cartridge.PRG[index + 1]
  elseif address >= 0x8000
    index = m.prgBank1 * Int32(0x4000) + Int32(address - 0x8000)
    cartridge.PRG[index + 1]
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1]
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper2 read at address: 0x%04X", address)
    0x00
  end
end

function write!(m::Mapper2, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address < 0x2000
    cartridge.CHR[address + 1] = val
  elseif address >= 0x8000
    m.prgBank1 = Int32(val) % m.prgBanks
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1] = val
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper2 write at address: 0x%04X", address)
  end
end

# Mapper 225

mutable struct Mapper225 <: Mapper
  chrBank::Int32
  prgBank1::Int32
  prgBank2::Int32

  Mapper225() = new(0, 0, 0)
end

function read(m::Mapper225, cartridge::Cartridge, address::UInt16)::UInt8
  if address < 0x2000
    index = m.chrBank * Int32(0x2000) + Int32(address)
    cartridge.CHR[index + 1]
  elseif address >= 0xC000
    index = m.prgBank2 * Int32(0x4000) + Int32(address-0xC000)
    cartridge.PRG[index + 1]
  elseif address >= 0x8000
    index = m.prgBank1 * Int32(0x4000) + Int32(address - 0x8000)
    cartridge.PRG[index + 1]
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1]
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper225 read at address: 0x%04X", address)
    0
  end
end

function write!(m::Mapper225, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address < 0x8000
    return
  end

  A = Int32(address)
  bank = (A >> 14) & 1
  m.chrBank = (A & 0x3f) | (bank << 6)
  prg = ((A >> 6) & 0x3f) | (bank << 6)
  mode = (A >> 12) & 1
  if mode == 1
    m.prgBank1 = prg
    m.prgBank2 = prg
  else
    m.prgBank1 = prg
    m.prgBank2 = prg + 1
  end
  mirr = (A >> 13) & 1
  if mirr == 1
    cartridge.mirror = UInt8(Horizontal::MirrorModes)
  else
    cartridge.mirror = UInt8(Vertical::MirrorModes)
  end
end

# Mapper 3

mutable struct Mapper3 <: Mapper
  chrBank::Int32
  prgBank1::Int32
  prgBank2::Int32

  function Mapper3(cartridge::Cartridge)
    m = new(0, 0, 0)
    m.prgBank2 = length(cartridge.PRG) ÷ 0x4000 - 1;
    m
  end
end

function read(m::Mapper3, cartridge::Cartridge, address::UInt16)::UInt8
  if address < 0x2000
    index = m.chrBank * Int32(0x2000) + address
    cartridge.CHR[index + 1]
  elseif address >= 0xC000
    index = m.prgBank2 * Int32(0x4000) + Int32(address - 0xC000)
    cartridge.PRG[index + 1]
  elseif address >= 0x8000
    index = m.prgBank1 * Int32(0x4000) + Int32(address - 0x8000)
    cartridge.PRG[index + 1]
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1]
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper3 read at address: 0x%04X", address)
    0
  end
end

function write!(m::Mapper3, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address < 0x2000
    index = m.chrBank * Int32(0x2000) + address
    cartridge.CHR[index + 1] = val
  elseif address >= 0x8000
    m.chrBank = val & 3
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1] = val
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper3 write at address: 0x%04X", address)
  end
end

# Mapper4

mutable struct Mapper4 <: Mapper
  reg::UInt8
  registers::Vector{UInt8}
  prgMode::UInt8
  chrMode::UInt8
  prgOffsets::Vector{Int32}
  chrOffsets::Vector{Int32}
  reload::UInt8
  counter::UInt8
  irqEnable::Bool

  function Mapper4(cartridge::Cartridge)
    m = new(0x00, zeros(UInt8, 8), 0x00, 0x00, zeros(Int32, 4), zeros(Int32, 8), 0x00, 0x00, false)
    m.prgOffsets[1] = prgbankoffset(m, cartridge, Int32(0))
    m.prgOffsets[2] = prgbankoffset(m, cartridge, Int32(1))
    m.prgOffsets[3] = prgbankoffset(m, cartridge, Int32(-2))
    m.prgOffsets[4] = prgbankoffset(m, cartridge, Int32(-1))
    m
  end
end

function mapperstep!(console::Console{Mapper4})
  m = console.mapper
  ppu = console.ppu
  if ppu.cycle != 280 # TODO: this *should* be 260
    return
  end
  if (ppu.scanline > 239) && (ppu.scanline < 261)
    return
  end
  if (ppu.flagShowBackground == 0) && (ppu.flagShowSprites == 0)
    return
  end
  handlescanline(m, console)
end

function read(m::Mapper4, cartridge::Cartridge, address::UInt16)::UInt8
  if address < 0x2000
    bank = address >> 10
    offset = address & 0x03FF
    cartridge.CHR[m.chrOffsets[bank + 1] + offset + 1]
  elseif address >= 0x8000
    address = address - 0x8000
    bank = address >> 13
    offset = address & 0x1FFF
    cartridge.PRG[m.prgOffsets[bank + 1] + offset + 1]
  elseif address >= 0x6000
    cartridge.SRAM[address - 0x6000 + 1]
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper4 read at address: 0x%04X", address)
    0
  end
end

function write!(m::Mapper4, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address < 0x2000
    bank = address >> 10
    offset = address & 0x03FF
    cartridge.CHR[m.chrOffsets[bank + 1] + offset + 1] = val
  elseif address >= 0x8000
    writeregister!(m, cartridge, address, val)
  elseif address >= 0x6000
    cartridge.SRAM[address - 0x6000 + 1] = val
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper4 write at address: 0x%04X", address)
  end
end

function handlescanline!(m::Mapper4, console::Console)
  if m.counter == 0
    m.counter = m.reload
  else
    m.counter -= 1
    if (m.counter == 0) && m.irqEnable
      triggerirq(console.cpu)
    end
  end
end

function writeregister!(m::Mapper4, cartridge::Cartridge, address::UInt16, val::UInt8)
  if (address <= 0x9FFF) && (address & 1 == 0)
    writebankselect!(m, cartridge, val)
  elseif (address <= 0x9FFF) && (address & 1 == 1)
    writebankdata!(m, cartridge, val)
  elseif (address <= 0xBFFF) && (address & 1 == 0)
    writemirror!(m, cartridge, val)
  elseif (address <= 0xBFFF) && (address & 1 == 1)
    writeprotect!(m, val)
  elseif (address <= 0xDFFF) && (address & 1 == 0)
    writeirqlatch!(m, val)
  elseif (address <= 0xDFFF) && (address & 1 == 1)
    writeirqreload!(m, val)
  elseif address & 1 == 0
    writeirqdisable!(m, val)
  elseif address & 1 == 1
    writeirqenable!(m, val)
  end
end

function writebankselect!(m::Mapper4, cartridge::Cartridge, val::UInt8)
  m.prgMode = (val >> 6) & 1
  m.chrMode = (val >> 7) & 1
  m.reg = val & 7
  updateoffsets!(m, cartridge)
end

function writebankdata!(m::Mapper4, cartridge::Cartridge, val::UInt8)
  m.registers[reg] = val
  updateoffsets!(m, cartridge)
end

function writemirror!(m::Mapper4, cartridge::Cartridge, val::UInt8)
  c = val & 1
  @match c begin
    0 => begin
      cartridge.mirror = UInt8(Vertical::MirrorModes)
    end
    1 => begin
      cartridge.mirror = UInt8(Horizontal::MirrorModes)
    end
    _ => throw(ErrorException("Unreachable code."))
  end
end

function writeprotect!(m::Mapper4, val::UInt8)
  # Not sure why this function exists...
  return
end

function writeirqlatch!(m::Mapper4, val::UInt8)
  m.reload = val
end

function writeirqreload!(m::Mapper4, val::UInt8)
  m.counter = 0
end

function writeirqdisable!(m::Mapper4, val::UInt8)
  m.irqEnable = false
end

function writeirqenable!(m::Mapper4, val::UInt8)
  m.irqEnable = true
end

function prgbankoffset(m::Mapper4, cartridge::Cartridge, index::Int32)::Int32
  if index >= 0x80
    index -= 0x100
  end
  index %= length(cartridge.PRG) >> 13
  offset = index * Int32(0x2000)
  if offset < 0
    offset += length(cartridge.PRG)
  end
  offset
end

function chrbankoffset(m::Mapper4, cartridge::Cartridge, index::Int32)::Int32
  if index >= 0x80
    index -= 0x100
  end
  index %= length(cartridge.CHR) >> 10
  offset = index * Int32(0x0400)
  if offset < 0
    offset += length(cartridge.CHR)
  end
  offset
end

function updateoffsets!(m::Mapper4, cartridge::Cartridge)
  @match m.prgMode begin
    0 => begin
      m.prgOffsets[1] = prgbankoffset(m, cartridge::Cartridge, Int32(m.registers[7]))
      m.prgOffsets[2] = prgbankoffset(m, cartridge::Cartridge, Int32(m.registers[8]))
      m.prgOffsets[3] = prgbankoffset(m, cartridge::Cartridge, Int32(-2))
      m.prgOffsets[4] = prgbankoffset(m, cartridge::Cartridge, Int32(-1))
    end
    1 => begin
      m.prgOffsets[1] = prgbankoffset(m, cartridge::Cartridge, Int32(-2))
      m.prgOffsets[2] = prgbankoffset(m, cartridge::Cartridge, Int32(m.registers[8]))
      m.prgOffsets[3] = prgbankoffset(m, cartridge::Cartridge, Int32(m.registers[7]))
      m.prgOffsets[4] = prgbankoffset(m, cartridge::Cartridge, Int32(-1))
    end
    _ => throw(ErrorException("Invalid prgMode."))
  end
  @match m.chrMode begin
    0 => begin
      m.chrOffsets[1] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[1] & 0xFE))
      m.chrOffsets[2] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[1] | 0x01))
      m.chrOffsets[3] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[2] & 0xFE))
      m.chrOffsets[4] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[2] | 0x01))
      m.chrOffsets[5] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[3]))
      m.chrOffsets[6] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[4]))
      m.chrOffsets[7] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[5]))
      m.chrOffsets[8] = chrbankoffset(m, cartridge::Cartridge, Int32(m.registers[6]))
    end
    1 => begin
      m.chrOffsets[1] = chrBankOffset(m, cartridge::Cartridge, Int32(m.registers[3]))
      m.chrOffsets[2] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[4]))
      m.chrOffsets[3] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[5]))
      m.chrOffsets[4] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[6]))
      m.chrOffsets[5] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[1] & 0xFE))
      m.chrOffsets[6] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[1] | 0x01))
      m.chrOffsets[7] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[2] & 0xFE))
      m.chrOffsets[8] = chrBankOffset(m, cartridge::Cartridge,  Int32(m.registers[2] | 0x01))
    end
    _ => throw(ErrorException("Invalid prgMode."))
  end
end

# Mapper7

mutable struct Mapper7 <: Mapper
  prgBank::Int32

  Mapper7() = new(0)
end

function read(m::Mapper7, cartridge::Cartridge, address::UInt16)::UInt8
  if address < 0x2000
    cartridge.CHR[address + 1]
  elseif address >= 0x8000
    index = m.prgBank * Int32(0x8000) + Int32(address - 0x8000)
    cartridge.PRG[index + 1]
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1]
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper7 read at address: 0x%04X", address)
    0
  end
end

function write!(m::Mapper7, cartridge::Cartridge, address::UInt16, val::UInt8)
  if address < 0x2000
    cartridge.CHR[address + 1] = val
  elseif address >= 0x8000
    m.prgBank = val & 7
    c = val & 0x10
    @match c begin
      0x00 => begin
        cartridge.mirror = UInt8(Single0::MirrorModes)
      end
      0x10 => begin
        cartridge.mirror = UInt8(Single1::MirrorModes)
      end
      _ => throw(ErrorException("Unreachable code"))
    end
  elseif address >= 0x6000
    index = address - 0x6000
    cartridge.SRAM[index + 1] = val
  else
    # TODO: add logging
    # log.Fatalf("unhandled mapper7 write at address: 0x%04X", address)
  end
end

function createmapper(cartridge::Cartridge)::Mapper
  @match cartridge.mapper begin
    0 => Mapper2(cartridge)
    1 => Mapper1(cartridge)
    2 => Mapper2(cartridge)
    3 => Mapper3(cartridge)
    4 => Mapper4(cartridge)
    7 => Mapper7()
    225 => Mapper225()
    _ => throw(ErrorException(@sprintf("Unimplemented mapper: %d", cartridge.mapper)))
  end
end


text:
using ...Network
import ...Network: act, state_dict, accept
using ...Signals

struct DelayerAgent{T <: AbstractFloat, U <: Unsigned, V <: Signal}
    delay::T # delay time
    w::AbstractFloat
    stack::Vector{Tuple{T, V}} # (time left, signal)
    donor::Address{U}
    acceptor::Address{U}
end

function DelayerAgent(
    delay::T,
    w::T,
    donor::Address{U},
    acceptor::Address{U},
    signal_name::String,
    network::Dict{String, Population{U}}
) where {T <: AbstractFloat, U <: Unsigned, V <: ForwardSignal}
    donor_check = can_add_acceptor(get_agent(network, donor), signal_name)
    acceptor_check = can_add_donor(get_agent(network, acceptor), signal_name)
    if donor_check && acceptor_check
        DelayerAgent{T, U, V}(delay, w, [], donor, acceptor)
    else
        error(string("DelayerAgent cannot be created for $(signal_name) on: ",
                     donor_check ? "" : "donor $(donor)",
                     ! donor_check && ! acceptor_check ? "and" : "",
                     acceptor_check ? "" : "acceptor $(acceptor)",
                     "."))
    end
end

function update(agent::DelayerAgent{T, U, V},
                updates::Vector{Tuple{T, V}}) where {T <: AbstractFloat, U <: Signal, V <: ForwardSignal}
    agent.stack = updates
end

## state_dict is for recording states, not for network-serialization. temporarily commented.
# function state_dict(agent::DelayerAgent{T}) where {T <: AbstractFloat, U <: ForwardSignal}
#     Dict(["stack" => map(state_dict, agent.stack)])
# end

function act(address::Address{U},
             agent::DelayerAgent{T, V},
             dt::T,
             trigger,
             push_signal) where {T <: AbstractFloat, U <: Unsigned, V <: ForwardSignal}

    new_stack, due_stack = reduce(
        agent.stack;
        Init=(Vector{(T, V)}(undef, 0), Vector{V}(undef, 0))
    ) do (acc, pair)
        new_t = pair[1] - dt
        if new_t >= 0
            ([acc[1]..., (new_t, pair[2])], acc[2])
        else
            (acc[1], [acc[2]..., amplify(pair[2], agent.w)])
        end
    end

    foreach(s -> push_signal(agent.acceptor, s), due_stack)
    agent.stack - new_stack

    if length(new_stack) > 0
        trigger(address)
    end

end

function accept(agent::DelayerAgent{T, U, V}, signal::U) where {T <: AbstractFloat, U <: Unsigned, V <: ForwardSignal}
    push!(agent.stack, (agent.delay, signal))
end

function can_add_donor(agent::DelayerAgent{T, U, V},
                       signal_name::String) where {T <: AbstractFloat, U <: Unsigned, V <: ForwardSignal}
    signal_name == name_dc_instruction
end

function accept(agent::DelayerAgent{T, U, V},
                signal::BackSpike) where {T <: AbstractFloat, U <: Unsigned, V <: ForwardSignal}
    # STDP
end


text:
module Geom

using Color
using Compose
using DataFrames
using Gadfly

import Compose.combine # Prevent DataFrame.combine from taking over.
import Gadfly.render, Gadfly.element_aesthetics, Gadfly.inherit, Gadfly.escape_id
import Iterators.cycle, Iterators.product, Iterators.distinct

# Geometry that renders nothing.
type Nil <: Gadfly.GeometryElement
end

const nil = Nil()

function render(geom::Nil, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
end

# Catchall
function default_statistic(::Gadfly.GeometryElement)
    Gadfly.Stat.identity
end

# Geometry which displays points at given (x, y) positions.
type PointGeometry <: Gadfly.GeometryElement
end

const point = PointGeometry()

function element_aesthetics(::PointGeometry)
    [:x, :y, :size, :color]
end

# Generate a form for a point geometry.
#
# Args:
#   geom: point geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::PointGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    default_aes.size = Measure[theme.default_point_size]
    aes = inherit(aes, default_aes)

    lw0 = convert(Compose.SimpleMeasure{Compose.MillimeterUnit}, theme.line_width)
    lw1 = 10 * lw0
    compose(circle(aes.x, aes.y, aes.size),
            fill(aes.color),
            stroke([theme.highlight_color(c) for c in aes.color]),
            linewidth(theme.line_width),
            d3embed(@sprintf(".on(\"mouseover\", geom_point_mouseover(%0.2f), false)",
                             lw1.value)),
            d3embed(@sprintf(".on(\"mouseout\", geom_point_mouseover(%0.2f), false)",
                             lw0.value)),
            svgclass([@sprintf("geometry color_%s", escape_id(aes.color_label(c)))
                      for c in aes.color]))
end

# Line geometry connects (x, y) coordinates with lines.
type LineGeometry <: Gadfly.GeometryElement
end

const line = LineGeometry()

function element_aesthetics(::LineGeometry)
    [:x, :y, :color]
end

# Render line geometry.
#
# Args:
#   geom: line geometry.
#   theme: the plot's theme.
#   aes: aesthetics.
#
# Returns:
#   A compose Form.
#
function render(geom::LineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.point", aes, :x, :y)
    Gadfly.assert_aesthetics_equal_length("Geom.point", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    if length(aes.color) == 1
        points = {(x, y) for (x, y) in zip(aes.x, aes.y)}
        sort!(points)
        form = lines(points...) <<
               stroke(aes.color[1]) <<
               svgclass("geometry")
    else
        # group points by color
        points = Dict{ColorValue, Array{(Float64, Float64),1}}()
        for (x, y, c) in zip(aes.x, aes.y, cycle(aes.color))
            if !haskey(points, c)
                points[c] = Array((Float64, Float64),0)
            end
            push!(points[c], (x, y))
        end

        forms = Array(Any, length(points))
        for (i, (c, c_points)) in enumerate(points)
            sort!(c_points)
            forms[i] = lines({(x, y) for (x, y) in c_points}...) <<
                            stroke(c) <<
                            svgclass(@sprintf("geometry color_%s", escape_id(aes.color_label(c))))
        end
        form = combine(forms...)
    end

    form << fill(nothing) << linewidth(theme.line_width)
end

# Bar geometry summarizes data as vertical bars.
type BarGeometry <: Gadfly.GeometryElement
    # How bars should be positioned if they are grouped by color.
    # Valid options are:
    #   :stack -> place bars on top of each other (default)
    #   :dodge -> place bar next to each other
    position::Symbol
end

const bar = BarGeometry(:stack)
const histogram = BarGeometry(:stack)

function element_aesthetics(::BarGeometry)
    [:x, :y, :color]
end

function default_statistic(::BarGeometry)
    Gadfly.Stat.histogram
end

# Render a single color bar chart
function render_colorless_bar(geom::BarGeometry,
                              theme::Gadfly.Theme,
                              aes::Gadfly.Aesthetics)
    bar_form = empty_form
    for (x_min, x_max, y) in zip(aes.x_min, aes.x_max, aes.y)
        geometry_id = Gadfly.unique_svg_id()
        bar_form |= compose(rectangle(x_min*cx + theme.bar_spacing/2, 0.0,
                                      (x_max - x_min)*cx - theme.bar_spacing, y),
                            svgid(geometry_id), svgclass("geometry"))
    end
    bar_form
end

# Render a bar chart grouped by discrete colors and stacked.
function render_colorful_stacked_bar(geom::BarGeometry,
                                     theme::Gadfly.Theme,
                                     aes::Gadfly.Aesthetics)
    bar_form = empty_form

    stack_height = Dict()
    for x_min in aes.x_min
        stack_height[x_min] = 0.0
    end

    for i in sortperm(aes.color.refs, rev=true)
        x_min, x_max, y, c = aes.x_min[i], aes.x_max[i], aes.y[i], aes.color[i]
        geometry_id = Gadfly.unique_svg_id()
        bar_form |= compose(rectangle(x_min*cx + theme.bar_spacing/2,
                                      stack_height[x_min],
                                      (x_max - x_min)*cx - theme.bar_spacing,
                                      y),
                            fill(c),
                            svgid(geometry_id),
                            svgclass("geometry"))
        stack_height[x_min] += y
    end
    bar_form
end

# Render a bar chart grouped by discrete colors and stuck next to each other.
function render_colorful_dodged_bar(geom::BarGeometry,
                                    theme::Gadfly.Theme,
                                    aes::Gadfly.Aesthetics)
    bar_form = empty_form
    dodge_width = Dict()
    dodge_count = Dict()
    for x_min in aes.x_min
        dodge_width[x_min] = 0.0cx
        dodge_count[x_min] = get(dodge_count, x_min, 0) + 1
    end

    for i in sortperm(aes.color.refs, rev=true)
        x_min, x_max, y, c = aes.x_min[i], aes.x_max[i], aes.y[i], aes.color[i]
        geometry_id = Gadfly.unique_svg_id()
        barwidth = ((x_max - x_min)*cx - theme.bar_spacing) / dodge_count[x_min]
        bar_form |= compose(rectangle(x_min*cx + dodge_width[x_min] + theme.bar_spacing/2,
                                      0.0,
                                      barwidth,
                                      y),
                            fill(c),
                            svgid(geometry_id),
                            svgclass("geometry"))
        dodge_width[x_min] += barwidth
    end
    bar_form
end

# Render bar geometry.
#
# Args:
#   geom: bar geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::BarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.bar", aes, :y)
    if (is(aes.x_min, nothing) || is(aes.x_max, nothing)) && is(aes.x, nothing)
        error("Geom.bar required \"x\" to be bound or both \"x_min\" and \"x_max\".")
    end

    if aes.x_min === nothing
        aes2 = Gadfly.Aesthetics()
        aes2.x_min = Array(Float64, length(aes.x))
        aes2.x_max = Array(Float64, length(aes.x))
        for (i, x) in enumerate(aes.x)
            aes2.x_min[i] = x - 0.5
            aes2.x_max[i] = x + 0.5
        end
        aes = inherit(aes, aes2)
    end

    if aes.color === nothing
        form = render_colorless_bar(geom, theme, aes)
    elseif geom.position == :stack
        form = render_colorful_stacked_bar(geom, theme, aes)
    elseif geom.position == :dodge
        form = render_colorful_dodged_bar(geom, theme, aes)
    else
        error("$(geom.position) is not a valid position for the bar geometry.")
    end

    compose(canvas(units_inherited=true),
            form,
            svgattribute("shape-rendering", "crispEdges"),
            fill(theme.default_color),
            stroke(nothing))
end

type RectangularBinGeometry <: Gadfly.GeometryElement
end

const rectbin = RectangularBinGeometry()

function element_aesthetics(::RectangularBinGeometry)
    [:x, :y, :x_min, :x_max, :y_min, :y_max, :color]
end

# Render a rectbin geometry with continuous x_min/x_max y_min/y_max coordinates.
function render_continuous_rectbin(geom::RectangularBinGeometry,
                                   theme::Gadfly.Theme,
                                   aes::Gadfly.Aesthetics)
    n = length(aes.x_min)
    forms = Array(Compose.Form, 0)

    for (i, c) in zip(1:n, cycle(aes.color))
        if !isna(c)
            push!(forms, rectangle(aes.x_min[i], aes.y_min[i],
                                  (aes.x_max[i] - aes.x_min[i])*cx - theme.bar_spacing,
                                  (aes.y_max[i] - aes.y_min[i])*cy + theme.bar_spacing) <<
                             fill(c) << svgclass("geometry"))
        end
    end

    compose(combine(forms...),
            stroke(nothing),
            svgattribute("shape-rendering", "crispEdges"))
end

# Rendere a rectbin geometry with discrete x/y coordinaes.
function render_discrete_rectbin(geom::RectangularBinGeometry,
                                   theme::Gadfly.Theme,
                                   aes::Gadfly.Aesthetics)
    n = length(aes.x)
    forms = Array(Compose.Form, 0)
    for (i, c) in zip(1:n, cycle(aes.color))
        if !isna(c)
            x, y = aes.x[i], aes.y[i]
            push!(forms, compose(rectangle(x - 0.5, y - 0.5, 1.0, 1.0),
                                 fill(c),
                                 svgclass("geometry")))
        end
    end

    compose(combine(forms...),
            stroke(nothing),
            svgattribute("shape-rendering", "crispEdges"))
end

# Render rectangular bin (e.g., heatmap) geometry.
#
# Args:
#   geom: rectbin geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::RectangularBinGeometry,
                theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    if aes.x_min === nothing
        Gadfly.assert_aesthetics_defined("Geom.bar", aes, :x, :y)
        Gadfly.assert_aesthetics_equal_length("Geom.bar", aes, :x, :y)
        render_discrete_rectbin(geom, theme, aes)
    else
        Gadfly.assert_aesthetics_defined("Geom.bar",
                                         aes, :x_min, :x_max, :y_min, :y_max)
        Gadfly.assert_aesthetics_equal_length("Geom.bar",
                                              aes, :x_min, :x_max,
                                              :y_min, :y_max)
        render_continuous_rectbin(geom, theme, aes)
    end
end

function default_statistic(::RectangularBinGeometry)
    Gadfly.Stat.rectbin
end

type BoxplotGeometry <: Gadfly.GeometryElement
end

const boxplot = BoxplotGeometry()

element_aesthetics(::BoxplotGeometry) = [:x, :y, :color]

default_statistic(::BoxplotGeometry) = Gadfly.Stat.boxplot

function render(geom::BoxplotGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.bar", aes,
                                     :lower_fence, :lower_hinge, :middle,
                                     :upper_hinge, :upper_fence, :outliers)
    Gadfly.assert_aesthetics_equal_length("Geom.bar", aes,
                                     :lower_fence, :lower_hinge, :middle,
                                     :upper_hinge, :upper_fence, :outliers)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    default_aes.x = Float64[1]
    aes = inherit(aes, default_aes)

    aes_iter = zip(aes.lower_fence,
                   aes.lower_hinge,
                   aes.middle,
                   aes.upper_hinge,
                   aes.upper_fence,
                   aes.outliers,
                   cycle(aes.x),
                   cycle(aes.color.refs))

    forms = Compose.Form[]
    middle_forms = Compose.Form[]

    r = theme.default_point_size
    bw = 1.0cx - theme.boxplot_spacing

    # TODO: handle color non-nothing color

    for (lf, lh, mid, uh, uf, outliers, x, cref) in aes_iter
        c = aes.color.pool[cref]
        sc = theme.highlight_color(c) # stroke color
        mc = theme.middle_color(c) # middle color

        # Middle
        push!(middle_forms, compose(lines((x - 1/6, mid), (x + 1/6, mid)),
                                    linewidth(theme.line_width), stroke(mc)))

        # Box
        push!(forms, compose(rectangle(x*cx - bw/2, lh, bw, uh - lh),
                            fill(c), stroke(sc),
                            linewidth(theme.highlight_width)))

        # Whiskers
        push!(forms, compose(lines((x, lh), (x, lf)),
                            linewidth(theme.line_width), stroke(sc)))

        push!(forms, compose(lines((x, uh), (x, uf)),
                            linewidth(theme.line_width), stroke(sc)))

        # Fences
        push!(forms, compose(lines((x - 1/6, lf), (x + 1/6, lf)),
                            linewidth(theme.line_width), stroke(sc)))

        push!(forms, compose(lines((x - 1/6, uf), (x + 1/6, uf)),
                            linewidth(theme.line_width), stroke(sc)))

        # Outliers
        if !isempty(outliers)
            push!(forms, compose(combine([circle(x, y, r) for y in outliers]...),
                                fill(c), stroke(sc)))
        end
    end

    compose(canvas(units_inherited=true),
            (canvas(units_inherited=true), combine(forms...)),
            (canvas(units_inherited=true, order=1), combine(middle_forms...)),
            svgclass("geometry"))
end

type LabelGeometry <: Gadfly.GeometryElement
end

element_aesthetics(::LabelGeometry) = [:x, :y, :label]

default_statistic(::LabelGeometry) = Gadfly.Stat.identity

const label = LabelGeometry()

# A deferred canvas function for labeling points in a plot. Optimizing label
# placement depends on knowing the absolute size of the containing canvas.
function deferred_label_canvas(aes, theme, box, unit_box)

    # Label layout is non-trivial problem. Quite a few papers and at least one
    # Phd thesis has been written on the topic. The approach here is pretty
    # simple. A label may be placed anywhere surrounding a point or hidden.
    # Simulated annealing is used to try to minimize a penalty, which is equal
    # to the number of overlapping or out fo bounds labels, plus terms
    # penalizing hidden labels.
    #
    # TODO:
    # Penalize to prefer certain label positions over others.

    canvas_width, canvas_height = box.width, box.height

    # This should maybe go in theme? Or should we be using Aesthetics.size?
    padding = 3mm

    point_positions = Array(Tuple, 0)
    for (x, y) in zip(aes.x, aes.y)
        x = absolute_measure(x*cx, unit_box, box)
        y = absolute_measure(y*cy, unit_box, box)
        push!(point_positions, (x, y))
    end

    extents = [text_extents(theme.point_label_font,
                            theme.point_label_font_size,
                            label)
               for label in aes.label]

    extents = [(width + padding, height + padding)
               for (width, height) in extents]

    positions = Gadfly.Maybe(BoundingBox)[]
    for (i, (text_width, text_height)) in enumerate(extents)
        x, y = point_positions[i]
        push!(positions, BoundingBox(x, y, text_width, text_height))
    end

    # TODO: use Aesthetics.size and/or theme.default_point_size
    for (x, y) in point_positions
        push!(positions, BoundingBox(x - 0.5mm, y - 0.5mm, 1mm, 1mm))
        push!(extents, (1mm, 1mm))
    end

    n = length(aes.label)

    # Return a box containing every point that the label could possibly overlap.
    function max_extents(i)
        BoundingBox(positions[i].x0 - extents[i][1],
                    positions[i].y0 - extents[i][2],
                    2*extents[i][1],
                    2*extents[i][2])
    end

    # True if two boxes overlap
    function overlaps(a, b)
        if a === nothing || b === nothing
            return false
        end

        a.x0 + a.width  >= b.x0 && a.x0 <= b.x0 + b.width &&
        a.y0 + a.height >= b.y0 && a.y0 <= b.y0 + b.height
    end

    # True if a is fully contained in box.
    function box_contains(a)
        if a === nothing
            return true
        end

        0mm < a.x0 && a.x0 + a.width < box.width &&
        0mm < a.y0 - a.height && a.y0 < box.height
    end

    # Checking for label overlaps is O(n^2). To mitigate these costs, we build a
    # sparse overlap matrix. This also costs O(n^2), but we only have to do it
    # once, rather than every iteration of annealing.
    possible_overlaps = [Array(Int, 0) for _ in 1:length(positions)]

    for j in 1:n
        for i in (j+1):n
            if overlaps(max_extents(i), max_extents(j))
                push!(possible_overlaps[i], j)
                push!(possible_overlaps[j], i)
            end
        end

        for i in (n+1):length(positions)
            # skip the point box corresponding to label
            if i == j + n
                continue
            end

            if overlaps(positions[i], max_extents(j))
                push!(possible_overlaps[i], j)
                push!(possible_overlaps[j], i)
            end
        end
    end

    # This variable holds the value of the objective function we wish to
    # minimize. A label overlap is a penalty of 1. Other penaties (out of bounds
    # labels, hidden labels) or calibrated to that.
    total_penalty = 0

    for i in 1:n
        if !box_contains(positions[i])
            total_penalty += theme.label_out_of_bounds_penalty
        end
    end

    for j in 1:n
        for i in possible_overlaps[j]
            if i > j && overlaps(positions[i], positions[j])
                total_penalty += 1
            end
        end
    end

    num_iterations = n * theme.label_placement_iterations
    for k in 1:num_iterations
        if total_penalty == 0
            break
        end
        j = rand(1:n)

        new_total_penalty = total_penalty

        # Propose flipping the visibility of the label.
        if !is(positions[j], nothing) && rand() < theme.label_visibility_flip_pr
            pos = nothing
            new_total_penalty += theme.label_hidden_penalty

        # Propose a change to label placement.
        else
            if positions[j] === nothing
                new_total_penalty -= theme.label_hidden_penalty
            end

            r = rand()
            point_x, point_y = point_positions[j]
            xspan = extents[j][1]
            yspan = extents[j][2]

            if rand() < 0.5
                xpos = Gadfly.lerp(rand(),
                                   (point_x - 7xspan/8).value,
                                   (point_x - 6xspan/8).value) * mm
            else
                xpos = Gadfly.lerp(rand(),
                                   (point_x - 2xspan/8).value,
                                   (point_x - 1xspan/8).value) * mm
            end

            ypos = Gadfly.lerp(rand(),
                               (point_y - 3yspan/4).value,
                               (point_y - 1yspan/4).value) * mm

            # choose a side
            if r < 0.25 # top
                pos = BoundingBox(xpos, point_y - extents[j][2],
                                  extents[j][1], extents[j][2])
            elseif 0.25 <= r < 0.5 # right
                pos = BoundingBox(point_x, ypos,
                                  extents[j][1], extents[j][2])
            elseif 0.5 <= r < 0.75 # bottom
                pos = BoundingBox(xpos, point_y,
                                  extents[j][1], extents[j][2])
            else # left
                pos = BoundingBox(point_x - extents[j][1], ypos,
                                  extents[j][1], extents[j][2])
            end
        end

        if !box_contains(positions[j])
            new_total_penalty -= theme.label_out_of_bounds_penalty
        end

        if !box_contains(pos)
            new_total_penalty += theme.label_out_of_bounds_penalty
        end

        for i in possible_overlaps[j]
            if overlaps(positions[i], positions[j])
                new_total_penalty -= 1
            end

            if overlaps(positions[i], pos)
                new_total_penalty += 1
            end
        end

        improvement = total_penalty - new_total_penalty
        T = 0.5 * (1.0 - (k / (1 + num_iterations)))
        if improvement >= 0 || rand() < exp(improvement / T)
            positions[j] = pos
            total_penalty = new_total_penalty
        end
    end

    forms = Array(Any, 0)

    # Quite useful for visually debugging this stuff:
    #for position in positions
        #if position === nothing
            #continue
        #end

        #push!(forms,
             #rectangle(position.x0, position.y0, position.width, position.height)
             #<< stroke("red") << fill(nothing) << linewidth(0.1mm))
    #end

    for i in 1:n
        if !is(positions[i], nothing)
            # Padding? The direction depends on what side we are on.
            point_x, point_y = point_positions[i]
            x, y = positions[i].x0, positions[i].y0

            x += extents[i][1] / 2
            y += extents[i][2] / 2

            push!(forms, text(x, y, aes.label[i], hcenter, vcenter) <<
                    svgclass("geometry"))
        end
    end

    #println("total_penalty = ", total_penalty)

    compose(canvas(unit_box=convert(Units, unit_box)),
            combine(forms...),
            font(theme.point_label_font),
            fontsize(theme.point_label_font_size),
            fill(theme.point_label_color),
            stroke(nothing))
end

function render(geom::LabelGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.Label", aes, :label, :x, :y)
    deferredcanvas((box, unit_box) -> deferred_label_canvas(aes, theme, box, unit_box))
end

end # module Geom


text:
using SynapseUploadUtils
uploadfolder(ARGS)


text:
include("lpvalue.jl")

function solutionofbase(X,C,A,b,VB)
    CB = hcat(map(x->C[findfirst(y->y==x,X)],VB)...)
    B = hcat(map(x->A[:,findfirst(y->y==x,X)[2]],VB)...)
    Binv = inv(B)
    Ap = Binv * A
    bp = Binv * b
    Cp = CB * Binv * A - C
    Z = CB * Binv * b
    return Cp, Z, Ap, bp
end

function findnewbase(X, VB, Cp, Ap, bp)
    possibleenterings = findall(x->x==min(Cp...)&&x<0,Cp)
    if length(possibleenterings) == 0 return :solved end
    #enteringid = possibleenterings[rand(1:length(possibleenterings))][2]
    enteringid = possibleenterings[1][2]
    reasons = ((x,y)->y <= 0 ? -1 : x/y).(bp, Ap[:,enteringid])
    acceptablereasons = filter(x->x>0,reasons)
    if length(acceptablereasons) == 0 return :unbounded end
    possibleexitings = findall(x->x==min(acceptablereasons...),reasons)
    #exitingid = possibleexitings[rand(1:length(possibleexitings))][1]
    exitingid = possibleexitings[1][1]
    return VB[exitingid] => X[enteringid]
end

function finddualbase(X, VB, Cp, Ap, bp)
    acceptableexitings = filter(x->x<0,bp)
    possibleexitings = findall(x->x<0&&abs(x)==max(abs.(acceptableexitings)...),bp)
    if length(possibleexitings) == 0 return :solved end
    exitingid = possibleexitings[1][1]
    reasons = ((x,y)->y < 0 ? x/y : 1).(Cp, transpose(Ap[exitingid,:]))
    acceptablereasons = filter(x->x<0, reasons)
    if length(acceptablereasons) == 0 return :unbounded end
    possibleenterings = findall(x->abs(x)==min(abs.(acceptablereasons)...),
        reasons)
    enteringid = possibleenterings[1][2]
    return VB[exitingid] => X[enteringid]
end

function solvesimplex(X,C,A,b,VB,logger,basemethod)
    Cp, Z, Ap, bp = solutionofbase(X,C,A,b,VB)
    status = :starting
    logger(X,C,A,b,VB,Cp,Z,Ap,bp,status)
    while true
        status = basemethod(X,VB,Cp,Ap,bp)
        if status == :unbounded || status == :solved break end
        VB = replace(VB,status)
        Cp, Z, Ap, bp = solutionofbase(X,C,A,b,VB)
        logger(X,C,A,b,VB,Cp,Z,Ap,bp,status)
    end
    logger(X,C,A,b,VB,Cp,Z,Ap,bp,status)
    return VB,Cp,Z,Ap,bp,status
end

function standardizerestriction(X,C,A,R,i)
    nX = copy(X)
    nC = copy(C)
    nA = copy(A)
    nR = copy(R)
    if R[i,1] == "<="
        nX = hcat(X, [ "S"*string(i) ])
        nC = hcat(C, [ 0 ])
        nA = hcat(A, [ k == i ? 1 : 0 for k=1:size(A)[1], l=[1] ])
        nR[i,1] = "="
    elseif R[i,1] == "="
        nX = hcat(X, [ "A"*string(i) ])
        nC = hcat(C, [ LPValue(1,0) ])
        nA = hcat(A, [ k == i ? 1 : 0 for k=1:size(A)[1], l=[1] ])
    elseif R[i,1] == ">="
        nX = hcat(X, [ "E"*string(i) "A"*string(i) ])
        nC = hcat(C, [ 0 LPValue(-1,0) ])
        nA = hcat(A, [ k == i ? -1 : 0 for k=1:size(A)[1], l=[1] ],
                     [ k == i ? 1 : 0 for k=1:size(A)[1], l=[1] ])
    end
    return nX, nC, nA, nR
end

function solveproblem(D,C,A,R,b,CR,twophases,dual)
    X = [ "X"*string(i) for j=[ 1 ], i=1:size(C)[2] ]
    basemethod = findnewbase
    inverted = false
    if !dual
        if D == "min"
            C *= -1
            inverted = true
        end
    else
        C *= -1
        A *= -1
        b *= -1
        R = map(x->x==">=" ? "<=" : (x=="=" ? "=" : ">="),R)
        basemethod = finddualbase
        inverted = true
    end
    for i in 1:size(A)[1]
        X, C, A, R = standardizerestriction(X, C, A, R, i)
    end
    VB = hcat([X[:,i] for j=[1],
        i=(findall(x->sum(x)==1,collect(eachcol(A))))]...)

    function logger(x,c,a,b,vb,cp,z,ap,bp,status)
        println(status)
        names = hcat([ "Z" ], x, "ld")
        funcobj = hcat([ inverted ? -1 : 1 ], cp, z)
        li = transpose(vb)
        lower = hcat(li, ap, bp)
        final = vcat(names, funcobj, lower)
        final = map(x->x isa Number && x < 10^-10 && x > -10^-10 ? 0 : x, final)
        display(final)
    end
    if twophases
        firstphasec = map(x->x[1]=='A' ? -1 : 0,X)
        firstsolution = solvesimplex(X, firstphasec, A, b, VB, logger,
            basemethod)[1]
        artificialindices = map(y->y[2],findall(x->x[1]=='A',X))
        newindices = transpose(filter(x->!(x in artificialindices),
            collect(1:size(A)[2])))
        X = hcat([X[i] for j = [ 1 ], i=newindices]...)
        A = hcat([A[:,i] for j = [ 1 ], i=newindices]...)
        C = hcat([C[i] for j = [ 1 ], i=newindices]...)
        solvesimplex(X, C, A, b, firstsolution, logger, basemethod)
    else
        solvesimplex(X, C, A, b, VB, logger, basemethod)
    end
end

max_stars_repo_name: mkualquiera/egypt

max_issues_repo_name: mkualquiera/egypt

max_forks_repo_name: mkualquiera/egypt