text:
using LightGraphs
using IterTools

function load_net(dir::AbstractString)
    nodeFfiles = readdir(dir)
    filter!(x->ismatch(r"\w+.txt$",x),nodeFfiles)
    # node_ids_byname = Dict{AbstractString,Integer}()
    node_ids_bynum = Dict{Integer,AbstractString}()
    node_tables = Dict{AbstractString,Array{AbstractString,1}}()

    for name in	nodeFfiles
        f = open(joinpath(dir,name))
        try
            node_id = split(name,['.','_'])[1:2]
            node_tables[node_id[2]] = readlines(f)
            node_number = parse(Int,node_id[1])
            # node_ids_byname[node_id[2]] = node_number
            node_ids_bynum[node_number] = String(node_id[2])
        finally
            close(f)
        end
    end
    allunique(values(node_ids_bynum))||error("Node identifiers must be unique")
    ids = [node_id for (node_num,node_id) in sort(collect(node_ids_bynum))]
    nodes = Vector{BNode2{Int8}}(length(node_ids_bynum))
    graph = DiGraph(length(node_ids_bynum))
    tau_matrix::SparseMatrixCSC{Int,Int} = spzeros(nv(graph),nv(graph))
    for (nsigma,id) in enumerate(ids)
        test_pattern = matchall(r"\d+", strip(node_tables[id][1]))
        !isempty(test_pattern)||error("Invalid value: First line of node file $id must contain a numeric value to set default state (optional extra field, separated by space, order number for semi-synchronous updating)")
        if length(test_pattern)==1
            def_state,sync_order = parse(Int,test_pattern[1]),1
        elseif length(test_pattern)==2
            def_state,sync_order = map(x->parse(Int,x),test_pattern)
        else
            error("Incorrect number of arguments: Less or equal than 2 arguments expected")
        end
        test_pattern = match(r"\d+", strip(node_tables[id][2]))
        max_states = test_pattern != nothing ? parse(Int,test_pattern.match) : error("Invalid value: Second line of node file $id must contain the maximum number of possible states")
        state_range = collect(range(0, max_states)) # All possible states is given by a sequence of consecutive integer numbers starting from 0
        sigma = nodes[nsigma] = BNode2{Int8}(id,def_state,state_range,sync_order)
        push!(sigma.thresholds,(max_states-1)/max_states*sigma.state_range[2:end]...)
        reg_ids = map(x->split(x,['_','.']),matchall(r"(\d+_)?[[:alnum:]]+(\.\d+)?",strip(node_tables[id][3])))
        reg_indices = Int[]
        for reg_id in reg_ids
            if length(reg_id) == 1 # If only the id is typed,
                reg_id_name = reg_id[1]
            elseif length(reg_id) == 2
                if ismatch(r"^\d+$",reg_id[1])
                    reg_id_name = reg_id[2]
                    tau = 0
                else
                    reg_id_name = reg_id[1]
                    tau = parse(reg_id[2])
                end
            elseif length(reg_id) == 3
                reg_id_name = reg_id[2]
                tau = parse(Int,reg_id[3])
            else
                error("Format error for list of regulators for node ",sigma.id)
            end
            index_reg = findfirst(ids,String(reg_id_name))
            index_reg!=0 || error("Unregistered node in list of regulators: There is no file for node $reg_id_name")
            push!(reg_indices,index_reg)
            tau_matrix[index_reg,nsigma] = tau
            add_edge!(graph,index_reg,nsigma)
        end
        sorted_reg_indices = sortperm(reg_indices)
        for reg_phrase in node_tables[sigma.id][4:end]
            r = matchall(r"\d+",reg_phrase)
            if !isempty(r)
                sigma.rule[[parse(Int,ri) for ri in r[1:(end-1)]][sorted_reg_indices]] = parse(Int,r[end])
            end
        end
    end
    if isempty(nonzeros(tau_matrix))
        return Net2{Int8}(nodes,graph)
    else
        return Net2{Int8}(nodes,graph,tau_matrix)
    end
end

function to_boolean!(net::Net2{Int8})
    netstateB = Vector{BNode2{Bool}}()
    ids = [sigma.id for sigma in net.nodes]
    #Add auxiliary nodes
    for (ni,sigma) in enumerate(net.nodes)
        max_state = length(sigma.state_range)
        if max_state>2
            def_state_index = findfirst(sort(sigma.state_range),sigma.def_state)
            max_digits = length(digits(max_state,2)) #transform max_state to base two to determine maximum nodes needed to reencode multistate nodes
            for s in max_digits:-1:1
                push!(netstateB,BNode2{Bool}(join([sigma.id,s-1,s]),digits(def_state_index-1,2,max_digits)[s],[false,true],sigma.sync_order))
                push!(sigma.thresholds,0.5) #Watch out, this might result in a big difference with respect to the original multistate node
            end
        else
            push!(netstate,BNode2{Bool}(sigma.id,sigma.def_state==1,[false,true],sigma.sync_order))
        end
    end
    graphB = DiGraph(length(netstateB))
    #Rewire network
    idsB = [sigma.id for sigma in netstateB]
    indexB = 0
    for (ni,sigma) in enumerate(net.nodes)
        max_state = length(sigma.state_range)
        max_digits = length(digits(max_state,2))
        in_neighbors(net.graph,ni)

        for reg_id in reg_ids
            if length(reg_id) == 1 # If only the id is typed,
                reg_id_name = reg_id[1]
            elseif length(reg_id) == 2
                if ismatch(r"^\d+$",reg_id[1])
                    reg_id_name = reg_id[2]
                    tau = 0
                else
                    reg_id_name = reg_id[1]
                    tau = parse(reg_id[2])
                end
            elseif length(reg_id) == 3
                reg_id_name = reg_id[2]
                tau = parse(Int,reg_id[3])
            else
                error("Format error for list of regulators for node ",sigma.id)
            end
            index_reg = findfirst(ids,reg_id_name)
            index_reg!=0 || error("Unregistered node in list of regulators: There is no file for node $reg_id_name")
            push!(reg_indices,index_reg)
            tau_matrix[index_reg,nsigma] = tau
            add_edge!(graph,index_reg,nsigma)
        end
        sorted_reg_indices = sortperm(reg_indices)
        for reg_phrase in node_tables[sigma.id][4:end]
            r = matchall(r"\d+",reg_phrase)
            if !isempty(r)
                sigma.rule[[parse(Int,ri) for ri in r[1:(end-1)]][sorted_reg_indices]] = parse(Int,r[end])
            end
        end
        indexB = indexB + max_digits
    end

    netB = Net2{Bool}(netstate,graph)

    return netB
end

function generate_RBN2(N::Integer,K::Integer;p_scheme::Symbol=:fix,tau_scheme::Symbol=:fix,params...)
    findfirst([:fix,:prob,:semifix],p_scheme)!=0 || error("Invalid scheme name: Choose from :fix,:prob,:semifix")
    p = 0.5
    tau = 0
    taus = [0]
    for (id,value) in params
        if id==:p
            p = value
        elseif id==:tau
            tau = value
        elseif id==:taus
            taus = value
        end
    end
    net_graph = random_regular_digraph(N, K; dir=:in, seed=-1)
    net_nodes = Vector{BNode2{Bool}}(N)
    for n in vertices(net_graph)
        net_nodes[n] = BNode2{Bool}(join(["n",string(n)]),rand([false,true]),[false,true])
    end
    println("Succesfully created ",length(net_nodes)," nodes")
    fill_tables!(net_nodes,net_graph,p_scheme,p)
    tau_matrix = fill_taus(net_graph,tau_scheme,tau=tau,taus=taus)
    return Net2{Bool}(net_nodes,net_graph,tau_matrix)
end

function fill_tables!(net_nodes::Vector{BNode2{Bool}},net_graph::DiGraph,p_scheme::Symbol,p::AbstractFloat)
    A = adjacency_matrix(net_graph)
    rows = rowvals(A)
    if p_scheme==:semifix
        for (ns,sigma) in enumerate(netnet_nodes)
            itp = product((net_nodes[sigmak].state_range for sigmak in rows[nzrange(A,ns)])...)
            table_size = length(itp)
            n_ones = round(Int,p*table_size)
            table_output = trues(n_ones)
            push!(table_output,bitrand(table_size-n_ones)...)
            shuffle!(table_output)
            sigma.rule = Dict(BitArray(collect(reg_phrase)) => table_output[i] for (i,reg_phrase) in enumerate(itp))::reg_function
        end
    elseif p_scheme==:prob
        for (ns,sigma) in enumerate(net_nodes)
            itp = product((net_nodes[sigmak].state_range for sigmak in rows[nzrange(A,ns)])...)
            table_output = rand(length(itp)) .< p
            sigma.rule = Dict(BitArray(collect(reg_phrase)) => table_output[i] for (i,reg_phrase) in enumerate(itp))::reg_function
        end
    elseif p_scheme==:fix
        for (ns,sigma) in enumerate(net_nodes)
            itp = product((net_nodes[sigmak].state_range for sigmak in rows[nzrange(A,ns)])...)
            table_size = length(itp)
            table_output = falses(table_size)
            n_ones = round(Int,p*table_size)
            table_output[1:n_ones] = true
            shuffle!(table_output)
            sigma.rule = Dict(BitArray(collect(reg_phrase)) => table_output[i] for (i,reg_phrase) in enumerate(itp))::reg_function
        end
    else
        error("Invalid scheme name: Choose from :fix,:prob,:semifix")
    end
    return nothing
end

function fill_taus(net_graph::DiGraph,tau_scheme::Symbol;tau_params...)
    tau = 0
    taus = [0]
    found = false
    tau_matrix::SparseMatrixCSC{Int,Int} = spzeros(nv(net_graph),nv(net_graph))
    if tau_scheme==:fix
        for (id,value) in tau_params
            if id==:tau
                tau = value
                found = true
                break
            end
        end
        found || error("Expected tau parameter")
        if tau != 0
            tau_matrix = adjacency_matrix(net_graph)
            tau_vals = nonzeros(tau_matrix)
            tau_vals[:] = tau
        end
        # for (nsigma,sigma) in enumerate(netstate)
        #     push!(sigma.tau_in,fill(tau,length(sigma.sigma_in))...)
        #     push!(sigma.sigma_in_state,zeros(Int8,length(sigma.sigma_in))...)
        #     sigma.sigma_in_state_mu = Vector{Float16}(length(sigma.sigma_in))
        # end
    elseif tau_scheme==:random
        for (id,value) in tau_params
            if id==:taus
                taus = value
                found = true
                break
            end
        end
        found || error("Expected taus parameter")
        tau_matrix = adjacency_matrix(net_graph)
        tau_vals = nonzeros(tau_matrix)
        tau_vals[:] = rand(taus,length(tau_vals))
        # for (nsigma,sigma) in enumerate(netstate)
        #     push!(sigma.tau_in,rand(taus,length(sigma.sigma_in))...)
        #     push!(sigma.sigma_in_state,zeros(Int8,length(sigma.sigma_in))...)
        #     sigma.sigma_in_state_mu = Vector{Float16}(length(sigma.sigma_in))
        # end
    else
        error("Invalid scheme name: Choose from :fix,:random")
    end
    # mem_size = maximum(vcat([sigma.tau_in for sigma in netstate]...)) + 1
    # for sigma in netstate
    #     sigma.state = zeros(Int8,mem_size)
    # end
    return tau_matrix
end

function knock_out!(net::Net2,nodeKO_id::Union{AbstractString,Int})
    index_KO = typeof(nodeKO_id)<:Int ? nodeKO_id : findfirst(sigma->sigma.id==nodeKO_id,net.nodes)
	0 < index_KO <= length(net.nodes) || error("Unexisting node: $nodeKO_id")
    nodeKO = net.nodes[index_KO]
    net.sim_buffer[index_KO,:] = nodeKO.def_state
    # nodeKO.thresholds = [nodeKO.thresholds[findfirst(sort([nodeKO.def_state,nodeKO.thresholds...]),nodeKO.def_state)]]
    nodeKO.state_range = [nodeKO.def_state]
    for k in keys(nodeKO.rule)
        nodeKO.rule[k] = nodeKO.def_state
    end
end

function knock_down!(net::Net2,nodeKD_id::Union{AbstractString,Int})
    index_KD = typeof(nodeKD_id)<:Int ? nodeKD_id : findfirst(sigma->sigma.id==nodeKD_id,net.nodes)
	0 < index_KD <= length(net.nodes) || error("Unexisting node: $nodeKD_id")
    nodeKD = net.nodes[index_KD]
    net.sim_buffer[index_KD,:] = nodeKD.def_state
    for k in keys(nodeKD.rule)
        nodeKD.rule[k] = nodeKD.def_state
    end
end

function fix_state!(net::Net2,const_node_id::Union{AbstractString,Int},const_state::Integer)
    index_const = typeof(const_node_id)<:Int ? const_node_id : findfirst(sigma->sigma.id==const_node_id,net.nodes)
	0 < index_const <= length(net.nodes) || error("Unexisting node: $const_node_id")
    node_const = net.nodes[index_const]
    const_state in node_const.state_range || error("Trying to set node $(node_const.id) to a forbidden state")
    net.sim_buffer[index_const,:] = const_state
    for k in keys(node_const.rule)
        node_const.rule[k] = const_state
    end
end


text:
using GenomicIntersections
using Test

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

ext: jl


text:
module BISC195Bioinformatics1

export normalizeDNA,
        composition,
        gc_content,
        complement,
        reverse_complement,
        parse_fasta

# # uncomment the following line if you intend to use BioSequences types
# using BioSequences

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("Invalid base $base")
    end
    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

# Your code here.
# Don't forget to export your functions!

"""
    composition(sequence)

Counts the number of each type of base (including 'N')
in a DNA sequence and returns a Dict

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> composition("ACCGGGTTTTN")
    Dict{Char,Int64} with 5 entries:
        'A' => 1
        'G' => 3
        'T' => 4
        'N' => 1
        'C' => 2

    julia> composition("AAX")
    ERROR: Invalid base, X
"""
function composition(sequence)
    sequence = normalizeDNA(sequence)
    comp = Dict('A' => 0, 'C' => 0, 'G' => 0, 'T' => 0, 'N' => 0)
    for base in sequence
        if base == 'A'
            comp['A'] = comp['A'] + 1
        elseif base == 'C'
            comp['C'] = comp['C'] + 1
        elseif base == 'G'
            comp['G'] = comp['G'] + 1
        elseif base == 'T'
            comp['T'] = comp['T'] + 1
        elseif base == 'N'
            comp['N'] = comp['N'] + 1
        end
    end
    return comp
end

"""
    gc_content(sequence)

Calculates the GC ratio of a DNA sequence.
The GC ratio is the total number of G and C bases divided by the total length of the sequence.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> gc_content("AATG")
    0.25

    julia> gc_content("cccggg") * 100
    100.0

    julia> gc_content("ATta")
    0.0

    julia> gc_content("ATty")
    Error: Invalid base Y encountered

    julia> gc_content("ATCGN")
    0.4
"""
function gc_content(sequence)
    basecomp = composition(sequence)
    return (basecomp['G'] + basecomp['C']) / length(sequence)
end

"""
    complement(base::Char)

Get the DNA complement of the provided base:

    A <-> T
    G <-> C

Accepts uppercase or lowercase `String` or `Char`,
but always returns an uppercase `Char`.
If a valid base is not provided, the function throws an error.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> complement('C')
    'G': ASCII/Unicode U+0047 (category Lu: Letter, uppercase)
"""
function complement(base::Char)
    complements = Dict("A" => 'T',
                       "T" => 'A',
                       "G" => 'C',
                       "C" => 'G',
                       "N" => 'N')
    base = uppercase(string(base))
    !(base in keys(complements)) && error("Invalid base $base")
    return complements[base]
end

"""
    complement(sequence::AbstractString)

Takes a DNA sequence and returns the complement
of that sequence.

Takes lowercase or uppercase sequences,
but always returns uppercase.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> complement("ATTN")
    "TAAN"

    julia> complement("ATTAGC")
    "TAATCG"
"""
function complement(sequence::AbstractString)
    seq = normalizeDNA(sequence)
    comp = map(complement, seq)
    return comp
end

"""
    reverse_complement(sequence)

Takes a DNA sequence and returns the reverse complement
of that sequence.

Takes lowercase or uppercase sequences,
but always returns uppercase.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> reverse_complement("AAATTT")
    "AAATTT"

    julia> reverse_complement("GCAT")
    "ATGC"

    julia> rc = reverse_complement("TTGGG");

    julia> println(rc)
    CCCAA

    julia> reverse_complement("ATTAGC")
    "GCTAAT"

    julia> reverse_complement("ATN")
    "NAT"
"""
function reverse_complement(sequence)
    seq = normalizeDNA(sequence)
    comp = map(complement, seq)
    reversecomp = reverse(comp)
    return reversecomp
end

"""
    function parse_fasta(path)

Reads a fasta-formated file and returns 2 vectors,
one containing headers,
the other containing the entire sequence as a `String`.

Examples
≡≡≡≡≡≡≡≡≡

    julia> ex1 = parse_fasta("data/ex1.fasta");

    julia> ex1[1]
    2-element Array{String,1}:
    "ex1.1 | easy"
    "ex1.2 | multiline"

    julia> ex1[2]
    2-element Array{String,1}:
    "AATTATAGC"
    "CGCCCCCCAGTCGGATT"

    julia> ex2 = parse_fasta("data/ex2.fasta");
    ERROR: invalid base H
"""
function parse_fasta(path)
    headers = String[]
    sequences = String[]
    sequence = []
    for line in eachline(path)
        if occursin('>', line)
            push!(headers, chop(line, head = 1, tail = 0))
            push!(sequences, join(sequence))
            sequence = []
        else
            push!(sequence, line)
        end
    end
    push!(sequences, join(sequence))
    popfirst!(sequences)
    normalizeDNA.(sequences)
    return (headers, sequences)
end

end # module BISC195Bioinformatics1


text:
type Katex
  val::String
  displayMode::Bool
end

Katex(ex::String) = Katex(ex, false)

function showmsg(obj::Katex, opts::Dict=Dict())
  args = Dict(:newnid   => getnid(),
              :compname => "katex",
              :params   => Dict(:expr => obj.val,
                                :options => Dict(:displayMode => obj.displayMode)),
              :deco     => opts)

  sendcurrent("append", args)
  nothing
end

function loadmsg(::Type{Katex})
  comppath = normpath(joinpath(@__DIR__, "../../client/katex/katex.js"))
  send(currentSession, 0,
       "load",
       Dict{Symbol,Any}(:assetname => "katex",
   		                  :assetpath => comppath))
  sleep(3) # give time to load
end

export Katex


text:
# `Reference`s will automatically be unpacked.
mutable struct Reference{T}
    data::T
end

@inline dereference(r::Reference) = getfield(r, :data)
@inline dereference(r) = r
@inline function ThreadingUtilities.load(p::Ptr{UInt}, ::Type{Reference{T}}, i) where {T}
    i, ref = ThreadingUtilities.load(p, Ptr{Reference{T}}, i)
    i, getfield(Base.unsafe_pointer_to_objref(ref)::Reference{T}, :data)::T
    # i, getfield(unsafe_load(ref)::Reference{T}, :data)::T
end
@inline function ThreadingUtilities.store!(p::Ptr{UInt}, r::Reference, i)
    ThreadingUtilities.store!(p + i, reinterpret(UInt, pointer_from_objref(r)))
    i + sizeof(UInt)
end
@inline function ThreadingUtilities.load(p::Ptr{UInt}, ::Type{PtrArray{S,D,T,N,C,B,R,X,O}}, i) where {S,D,T,N,C,B,R,X,O}
    (i, sp) = ThreadingUtilities.load(p, StridedPointer{T,N,C,B,R,X,O}, i)
    (i, sz) = ThreadingUtilities.load(p, S, i)
    i, PtrArray(sp, sz, Val{D}())
end
@inline function ThreadingUtilities.store!(p::Ptr{UInt}, A::PtrArray, i)
    i = ThreadingUtilities.store!(p, stridedpointer(A), i)
    ThreadingUtilities.store!(p, size(A), i)
end

@generated _isbits(::T) where {T} = Expr(:call, isbitstype(T) ? :True : :False)

@inline object_and_preserve(x) = _object_and_preserve(x, _isbits(x))
@inline _object_and_preserve(x, ::True) = (x, nothing)
@inline function _object_and_preserve(x, ::False)
    r = Reference(x)
    r, r
end

@inline object_and_preserve(A::AbstractArray{T}) where {T<:VectorizationBase.NativeTypes} = array_object_and_preserve(ArrayInterface.device(A), A)
@inline object_and_preserve(A::AbstractArray) = _object_and_preserve(A, False())
@inline array_object_and_preserve(::ArrayInterface.CPUPointer, A::AbstractArray) = (PtrArray(A), preserve_buffer(A))
@inline array_object_and_preserve(_, A::AbstractArray) = _object_and_preserve(A, False())

@inline object_and_preserve(r::Base.RefValue{T}) where {T} = (Base.unsafe_convert(Ptr{T}, r), r)


text:
abstract type Data end

struct IncludeExcludePairs{V}
    mode::Symbol
    point_ids::Vector{V}
    include_pairs::Matrix{V}
end
function IncludeExcludePairs(V)
    IncludeExcludePairs(:undef, V[], Matrix{V}(undef,0,0))
end

struct NetworkData{T,V} <: Data
    coords::Tuple{Vector{V},Vector{V},Vector{T}}
    fp::Vector{V}
    source_map::Matrix{T}
    ground_map::Matrix{T}
end

struct RasterMeta
    ncols::Int
    nrows::Int
    xllcorner::Float64
    yllcorner::Float64
    cellsize::Float64
    nodata::Float64
    file_type::Int
    crs::String
    affine_map::AffineMap
end
function RasterMeta()
    RasterMeta(0,0,0,0,0,0,0,"",AffineMap([0 0; 0 0], [0, 0]))
end

struct RasData{T,V} <: Data
    cellmap::Matrix{T}
    polymap::Matrix{V}
    source_map::Matrix{T}
    ground_map::Matrix{T}
    points_rc::Tuple{Vector{V},Vector{V},Vector{V}}
    strengths::Matrix{T}
    included_pairs::IncludeExcludePairs{V}
    hbmeta::RasterMeta
end

function load_graph(V, gpath::String, ::Type{T}) where {T}
    g = readdlm(gpath, T)
    i = zeros(V, size(g, 1))
    j = zeros(V, size(g, 1))
    v = zeros(T, size(g, 1))
    for iter = 1:size(g, 1)
        i[iter] = g[iter,1] + 1
        j[iter] = g[iter,2] + 1
        v[iter] = g[iter,3]
    end
    i,j,v
end

read_focal_points(V, path::String) = V.(vec(readdlm(path)) .+ 1)

function read_point_strengths(T, path::String, inc = true)
    a = readdlm(path, T)
    if inc
        @. a[:,1] = a[:,1] + 1
    end
    a
end

function read_cellmap(habitat_file::String, is_res::Bool, ::Type{T}) where {T}

    if is_geotiff(habitat_file)
        cell_map, rastermeta = read_geotiff(T, habitat_file)
    else
        cell_map, rastermeta = read_ascii_grid(T, habitat_file)
    end

    gmap = similar(cell_map)
    ind = findall(x -> x == -9999, cell_map)
    if is_res
        if count(x -> x == 0, cell_map) > 0
            throw("Error: zero resistance values are not currently supported for habitat maps. Use a short-circuit region file instead.")
        else
            for i in eachindex(cell_map)
                gmap[i] = 1 ./ cell_map[i]
            end
            gmap[ind] .= 0
        end
    else
        copyto!(gmap, cell_map)
        gmap[ind] .= 0
    end
    gmap, rastermeta
end

function _ascii_grid_reader(T, file)
    f = endswith(file, ".gz") ? GZip.open(file, "r") : open(file, "r")
    rastermeta = _ascii_grid_read_header(file, f)
    c = Matrix{T}(undef,0,0)
    ss = 6
    if rastermeta.nodata == -Inf
        ss = 5
    end
    try
        c = readdlm(f, T; skipstart = ss)
    catch
        seek(f, 0)
        try
            d = readdlm(f; skipstart = ss)
            d = d[:, 1:end-1]
            c = map(T, d)
        catch
            error("Failed to read habitat map. There may be errors in your file.")
        end
    end
    map!(x -> x == rastermeta.nodata ? -9999. : x , c, c)
    close(f)
    c, rastermeta
end

function _ascii_grid_read_header(habitat_file, f)
    file_type = _guess_file_type(habitat_file, f)
    ncols = parse(Int, split(readline(f))[2])
    nrows = parse(Int, split(readline(f))[2])
    xllcorner = parse(Float64, split(readline(f))[2])
    yllcorner = parse(Float64, split(readline(f))[2])
    cellsize = parse(Float64, split(readline(f))[2])
    nodata = -Inf
    s = split(readline(f))
    if occursin("NODATA", s[1]) || occursin("nodata", s[1])
        nodata = parse(Float64, s[2])
    end
    seek(f, 0)
    RasterMeta(ncols, nrows, xllcorner, yllcorner, cellsize, nodata, file_type ,"", AffineMap([0 0; 0 0], [0, 0]))
end

#=function _guess_file_type(filename, f)
    s = readline(f)
    seek(f, 0)

    if startswith(s, "min")
        return PAIRS_AAGRID
    elseif startswith(s, "mode")
        return PAIRS_LIST
    elseif occursin(".asc", filename)
        return AAGRID
    elseif endswith(filename, ".txt")
        return TXTLIST
    else
        throw("Check file format")
    end

end=#

function is_geotiff(filename)
    endswith(filename,"tif") | endswith(filename,"tiff")
end

function is_aagrid(filename)
    endswith(filename,"asc")
end

function read_geotiff(T, f)
    gt=GeoArrays.read(f)
    ncols = size(gt,1)
    nrows = size(gt,2)
    xll = gt.f.translation[1]
    yll= gt.f.translation[2]-gt.f.linear[1,1]-(gt.f.linear[1,1]*(size(gt,2)-1))
    cellsize = gt.f.linear[1,1]
    crs = gt.crs
    affine_map = gt.f
    file_type = FILE_TYPE_GEOTIFF
    nodata = -Inf
    Array{T}(permutedims(gt.A[:,:,1])), RasterMeta(ncols, nrows, xll, yll, cellsize, nodata, file_type, crs, affine_map)
end

function read_ascii_grid(T, f)
    gt=GeoArrays.read(f)
    ncols = size(gt,1)
    nrows = size(gt,2)
    xll = gt.f.translation[1]
    yll= gt.f.translation[2]-gt.f.linear[1,1]-(gt.f.linear[1,1]*(size(gt,2)-1))
    cellsize = gt.f.linear[1,1]
    crs = gt.crs
    affine_map = gt.f
    file_type = FILE_TYPE_AAGRID
    nodata = -Inf
    Array{T}(permutedims(gt.A[:,:,1])), RasterMeta(ncols, nrows, xll, yll, cellsize, nodata, file_type, crs, affine_map)
end

function _guess_file_type(filename, f)

    if is_geotiff(filename)
        filetype = FILE_TYPE_GEOTIFF
    elseif is_aagrid(filename)
        filetype = FILE_TYPE_AAGRID
    else
        hdr = readline(f)
        seek(f, 0)

        if startswith(hdr, FILE_HDR_NPY)
            filetype = FILE_TYPE_NPY
        elseif startswith(lowercase(hdr), FILE_HDR_AAGRID)
            filetype = FILE_TYPE_AAGRID
        elseif startswith(hdr, FILE_HDR_INCL_PAIRS_AAGRID)
            filetype = FILE_TYPE_INCL_PAIRS_AAGRID
        elseif startswith(hdr, FILE_HDR_INCL_PAIRS)
            filetype = FILE_TYPE_INCL_PAIRS
        else
            filetype = FILE_TYPE_TXTLIST
        end
    end

    return filetype
end

function read_polymap(T, file::String, habitatmeta;
                            nodata_as = 0, resample = true)

    if is_geotiff(file)
        polymap, rastermeta = read_geotiff(T, file)
    else
        polymap, rastermeta = read_ascii_grid(T, file)
    end

    ind = findall(x -> x == rastermeta.nodata, polymap)
    if nodata_as != -1
        polymap[ind] .= nodata_as
    end

    if(is_geotiff(file) || is_aagrid(file))
        if rastermeta.yllcorner != habitatmeta.yllcorner
            cswarn("yllcorner is not the same")
        elseif rastermeta.xllcorner != habitatmeta.xllcorner
            cswarn("xllcorner is not the same")
        elseif rastermeta.ncols != habitatmeta.ncols
            cswarn("ncols is not the same")
        elseif rastermeta.nrows != habitatmeta.nrows
            cswarn("nrows is not the same")
        elseif rastermeta.cellsize != habitatmeta.cellsize
            cswarn("cellsize is not the same")
        elseif rastermeta.crs != habitatmeta.crs
            cswarn("reference system is not the same")
        end
    else
        if rastermeta.cellsize != habitatmeta.cellsize
            cswarn("cellsize is not the same")
        elseif rastermeta.ncols != habitatmeta.ncols
            cswarn("ncols is not the same")
        elseif rastermeta.nrows != habitatmeta.nrows
            cswarn("nrows is not the same")
        elseif rastermeta.yllcorner != habitatmeta.yllcorner
            cswarn("yllcorner is not the same")
        elseif rastermeta.xllcorner != habitatmeta.xllcorner
            cswarn("xllcorner is not the same")
        end
    end

    polymap
end

function read_point_map(V, file, habitatmeta)

    # Advanced mode
    if file == "none"
        return (V[], V[], V[])
    end

    f = endswith(file, ".gz") ? GZip.open(file, "r") : open(file, "r")
    filetype = _guess_file_type(file, f)
    _points_rc = filetype == FILE_TYPE_TXTLIST ? readdlm(file) :
                        read_polymap(V, file, habitatmeta)

    i = V[]
    j = V[]
    v = V[]
    if filetype == FILE_TYPE_TXTLIST
        I = _points_rc[:,2]
        J = _points_rc[:,3]
        v = _points_rc[:,1]
        i  = ceil.(V, habitatmeta.nrows .- (J .- habitatmeta.yllcorner) ./ habitatmeta.cellsize)
        j = ceil.(V, (I .- habitatmeta.xllcorner) ./ habitatmeta.cellsize)
    else
        _I = findall(!iszero, _points_rc)
        (i,j,v) =  (getindex.(_I, 1), getindex.(_I, 2), _points_rc[_I])
    end

    ind = findall(x -> x < 0, v)

    # Get rid of negative resistances
    for index in ind
        deleteat!(i, index)
        deleteat!(j, index)
        deleteat!(v, index)
    end

    # Sort them
    idx = sortperm(v)
    i = i[idx]
    j = j[idx]
    v = v[idx]

    if (minimum(i) < 0) || (minimum(j) < 0) ||
            (maximum(i) > (habitatmeta.nrows)) ||
            (maximum(j) > (habitatmeta.ncols))
         throw("At least one focal node location falls outside of habitat map")
    end

    if size(unique(v),1) < 2
        throw("Less than two valid focal nodes found. Please check focal node location file.")
    end

    close(f)
    V.(i), V.(j), V.(v)
end

function read_source_and_ground_maps(T, V, source_file, ground_file, habitatmeta,
                                        is_res)

    ground_map = Matrix{T}(undef,0,0)
    source_map = Matrix{T}(undef,0,0)

    f = endswith(ground_file, "gz") ? Gzip.open(ground_file, "r") : open(ground_file, "r")
    filetype = _guess_file_type(ground_file, f)

    if filetype == FILE_TYPE_AAGRID || filetype == FILE_TYPE_GEOTIFF
        ground_map = read_polymap(T, ground_file, habitatmeta; nodata_as = -1)
        ground_map = map(T, ground_map)
    else
        rc = readdlm(ground_file, V)
        ground_map = -9999 * ones(T, habitatmeta.nrows, habitatmeta.ncols)
        ground_map[rc[:,2], rc[:,3]] = rc[:,1]
    end
    close(f)

    f = endswith(source_file, "gz") ? Gzip.open(source_file, "r") : open(source_file, "r")
    filetype = _guess_file_type(source_file, f)

    if filetype == FILE_TYPE_AAGRID || filetype == FILE_TYPE_GEOTIFF
        source_map = read_polymap(T, source_file, habitatmeta)
        source_map = map(T, source_map)
    else
        rc = readdlm(source_file, V)
        source_map = -9999 * ones(T, habitatmeta.nrows, habitatmeta.ncols)
        source_map[rc[:,2], rc[:,3]] = rc[:,1]
    end
    close(f)

    if is_res
        ind = findall(x -> x == -9999, ground_map)
        ground_map = 1 ./ ground_map
        ground_map[ind] .= 0
    else
        ind = findall(x -> x == -9999, ground_map)
        ground_map[ind] .= 0
    end

    source_map, ground_map
end

function read_included_pairs(V, filename)

    f = endswith(filename, "gz") ? Gzip.open(filename, "r") : open(filename, "r")
    filetype = _guess_file_type(filename, f)
    minval = 0
    maxval = 0
    mode = :undef

    if filetype == FILE_TYPE_INCL_PAIRS_AAGRID
        open(filename, "r") do fV
            minval = parse(Float64, split(readline(f))[2])
            maxval = parse(Float64, split(readline(f))[2])
        end
        included_pairs = readdlm(filename, skipstart=2)
        point_ids = V.(included_pairs[:,1])
        deleteat!(point_ids, 1)
        included_pairs = included_pairs[2:end, 2:end]
        map!(x -> x > maxval ? 0 : x, included_pairs, included_pairs)
        idx = findall(x -> x >= minval, included_pairs)
        mode = :include
        bin = map(x -> x >= minval ? V(1) : V(0), included_pairs)
        return IncludeExcludePairs(mode, point_ids, bin)
    elseif filetype == FILE_TYPE_INCL_PAIRS
        open(filename, "r") do f
            mode = Symbol(split(readline(f))[2])
        end
        included_pairs = readdlm(filename, skipstart = 1)
        if size(included_pairs, 1) == 1
            pl = zeros(V, 1,2)
            pl[1,:] = included_pairs
            included_pairs = pl
        end
        point_ids = V.(sort!(unique(included_pairs)))
        idx = findall(x -> x == 0 , point_ids)
        if length(idx) > 0
            deleteat!(point_ids, idx)
            cswarn("Code to include pairs is activated, some entries did not match with focal node file. Some focal nodes may have been dropped")
        end

        mat = zeros(V, size(point_ids, 1), size(point_ids, 1))

        for i = 1:size(included_pairs, 1)
            idx1 = findfirst(x -> x == included_pairs[i, 1], point_ids)
            idx2 = findfirst(x -> x == included_pairs[i, 2], point_ids)
            if idx1 != nothing && idx2 != nothing
                mat[idx1,idx2] = 1
                mat[idx2,idx1] = 1
            end
        end

        return IncludeExcludePairs(mode, point_ids, mat)
    end

end

function get_network_data(T, V, cfg)::NetworkData{T,V}

    hab_is_res = cfg["habitat_map_is_resistances"] in TRUELIST
    hab_file = cfg["habitat_file"]
    fp_file = cfg["point_file"]
    source_file = cfg["source_file"]
    ground_file = cfg["ground_file"]

    is_pairwise = cfg["scenario"] in PAIRWISE

    i,j,v = load_graph(V, hab_file, T)
    if hab_is_res
        v = 1 ./ v
    end

    if is_pairwise
        fp = read_focal_points(V, fp_file)
    else
        fp = V[]
    end

    if !is_pairwise
        source_map = read_point_strengths(T, source_file)
        ground_map = read_point_strengths(T, ground_file)
    else
        source_map = Matrix{T}(undef,0,0)
        ground_map = Matrix{T}(undef,0,0)
    end

    NetworkData((i,j,v), fp, source_map, ground_map)
end

function load_raster_data(T, V, cfg)::RasData{T,V}

    # Habitat file
    hab_file = cfg["habitat_file"]
    hab_is_res = cfg["habitat_map_is_resistances"] in TRUELIST

    # Polygons
    use_polygons = cfg["use_polygons"] in TRUELIST
    polygon_file = cfg["polygon_file"]

    # Mask file
    use_mask = cfg["use_mask"] in TRUELIST
    mask_file = cfg["mask_file"]

    # Point file
    point_file = cfg["point_file"]

    # Variable source strengths
    use_var_source = cfg["use_variable_source_strengths"] in TRUELIST
    var_source_file = cfg["variable_source_file"]

    # Included Pairs
    use_inc_pairs = cfg["use_included_pairs"] in TRUELIST
    inc_pairs_file = cfg["included_pairs_file"]

    # Advanced mode
    is_pairwise = cfg["scenario"] in PAIRWISE
    is_advanced = cfg["scenario"] in ADVANCED
    source_file = cfg["source_file"]
    ground_file = cfg["ground_file"]
    ground_is_res = cfg["ground_file_is_resistances"] in TRUELIST

    csinfo("Reading maps")

    # Read cell map
    cellmap, hbmeta = read_cellmap(hab_file, hab_is_res, T)
    c = count(x -> x > 0, cellmap)
    csinfo("Resistance/Conductance map has $c nodes")

    # Read polymap
    if use_polygons
        polymap = read_polymap(V, polygon_file, hbmeta)
    else
        polymap = Matrix{V}(undef,0,0)
    end

    # Read and update cellmap with mask file
    if use_mask
        update!(cellmap, mask_file, hbmeta)
        sum(cellmap) == 0 && throw("Mask file deleted everything!")
    end

    # Read point file
    if !is_advanced
        points_rc = read_point_map(V, point_file, hbmeta)
    else
        points_rc = (V[], V[], V[])
    end

    # Advanced mode reading
    if is_advanced
        source_map, ground_map =
        read_source_and_ground_maps(T, V, source_file, ground_file,
                                    hbmeta, ground_is_res)
    else
        source_map, ground_map = Matrix{T}(undef,0,0), Matrix{T}(undef,0,0)
    end

    # Included Pairs
    if use_inc_pairs
        included_pairs = read_included_pairs(V, inc_pairs_file)
    else
        included_pairs = IncludeExcludePairs(V)
    end

    # Variable source strengths
    if use_var_source
        strengths = read_point_strengths(T, var_source_file)
    else
        strengths = Matrix{T}(undef, 0,0)
    end

    RasData(cellmap, polymap, source_map, ground_map, points_rc, strengths,
                    included_pairs, hbmeta)
end

function update!(cellmap::Matrix{T}, m::String, hbmeta) where {T}
    mask = read_polymap(T, m, hbmeta)
    map!(x -> x > 0 ? 1 : 0, mask, mask)
    cellmap .= cellmap .* mask
end


text:
"""
   HeistRecordSpec is an object that holds information about how to access a
   data product from gallery.

# Arguments

   You construct it by passing in the following *String* arguments
   with keywords (! are required):

         namespace is the namespace of the product, e.g. "gm2truth"

       ! productType is the C++ type of the product. It can be a typedef to a vector like
            CalorimeterArtRecordCollection

       ! header is the path to the header for the product. It can be a relative path if you have issued
         `addHeaderDir` (that's been done for `gm2dataproducts`).

       ! moduleLabel is the module label name associated with the product

         instanceLabel is the instance label name associated with the product

# Examples

```julia
trackRS = HeistRecordSpec(
             namespace     = "gm2truth",
             productType   = "TrackingActionArtRecordCollection",
             header        = "gm2dataproducts/mc/actions/track/TrackingActionArtRecord.hh",
             moduleLabel   = "artg4"
           )
```
"""

struct HeistRecordSpec
   productName::String
   artInputTag::Cxx.CppValue

   function HeistRecordSpec(; namespace="", productType="", header="", moduleLabel="", instanceLabel="")
     @assert productType != ""
     @assert moduleLabel != ""
     @assert header != ""

     cxxparse("#include \"" * header * "\"")
     productName = namespace == "" ? productType : namespace * "::" * productType
     inputTagString = instanceLabel == "" ? moduleLabel : moduleLabel * ":" * instanceLabel
     new(productName, icxx"""art::InputTag it($inputTagString); return it;""")
   end

end


text:
import Base: linspace

# Multidimensional Linspace
function linspace{D, T}(start ::NTuple{D, T}, stop ::NTuple{D, T}, n=50)
  zip([linspace(start[d], stop[d], n+1)[1:n] for d in 1:D]...)
end

function flatten{T}(a :: Array{T})
  return reshape(a, prod(size(a)))
end

import Iterators

function linpath{D, T <: Number}(points ::Vector{NTuple{D, T}}, n=50)
  segments = []
  for ip in 2:length(points)
    push!(segments, linspace(points[ip-1], points[ip], n))
  end
  Iterators.chain(segments...)
end


text:
using HttpServer

http = HttpHandler() do req::Request, res::Response
    Response("Hello World")
end
http.events["listen"] = (saddr) -> println("Running on https://$saddr (Press CTRL+C to quit)")

server = Server(http)
run(server, host=IPv4(127,0,0,1), port=8000)


text:
using SparseArrays
using LinearAlgebra
using UnicodePlots
using WriteVTK
using DiagonalSBP

# flatten tuples to arrays
flatten_tuples(x) =
    reshape(collect(Iterators.flatten(x)), length(x[1]), length(x))

⊗(A, B) = kron(A, B)

const BC_DIRICHLET = 1
const BC_NEUMANN = 2
const BC_LOCKED_INTERFACE = 0
const BC_JUMP_INTERFACE = 7

# Transfinite Blend
function transfinite_blend(α1, α2, α3, α4, α1s, α2s, α3r, α4r, r, s)
    # +---4---+
    # |       |
    # 1       2
    # |       |
    # +---3---+
    @assert [α1(-1) α2(-1) α1(1) α2(1)] ≈ [α3(-1) α3(1) α4(-1) α4(1)]

    x =
        (1 .+ r) .* α2(s) / 2 +
        (1 .- r) .* α1(s) / 2 +
        (1 .+ s) .* α4(r) / 2 +
        (1 .- s) .* α3(r) / 2 -
        (
            (1 .+ r) .* (1 .+ s) .* α2(1) +
            (1 .- r) .* (1 .+ s) .* α1(1) +
            (1 .+ r) .* (1 .- s) .* α2(-1) +
            (1 .- r) .* (1 .- s) .* α1(-1)
        ) / 4

    xr =
        α2(s) / 2 - α1(s) / 2 +
        (1 .+ s) .* α4r(r) / 2 +
        (1 .- s) .* α3r(r) / 2 -
        (
            +(1 .+ s) .* α2(1) +
            -(1 .+ s) .* α1(1) +
            +(1 .- s) .* α2(-1) +
            -(1 .- s) .* α1(-1)
        ) / 4

    xs =
        (1 .+ r) .* α2s(s) / 2 + (1 .- r) .* α1s(s) / 2 + α4(r) / 2 -
        α3(r) / 2 -
        (
            +(1 .+ r) .* α2(1) +
            +(1 .- r) .* α1(1) +
            -(1 .+ r) .* α2(-1) +
            -(1 .- r) .* α1(-1)
        ) / 4

    return (x, xr, xs)
end

function transfinite_blend(α1, α2, α3, α4, r, s, p)
    (Nrp, Nsp) = size(r)
    (Dr, _, _, _) = DiagonalSBP.D1(p, Nrp - 1; xc = (-1, 1))
    (Ds, _, _, _) = DiagonalSBP.D1(p, Nsp - 1; xc = (-1, 1))

    α2s(s) = α2(s) * Ds'
    α1s(s) = α1(s) * Ds'
    α4r(s) = Dr * α4(r)
    α3r(s) = Dr * α3(r)

    transfinite_blend(α1, α2, α3, α4, α1s, α2s, α3r, α4r, r, s)
end

function transfinite_blend(
    v1::T1,
    v2::T2,
    v3::T3,
    v4::T4,
    r,
    s,
) where {T1 <: Number, T2 <: Number, T3 <: Number, T4 <: Number}
    e1(α) = v1 * (1 .- α) / 2 + v3 * (1 .+ α) / 2
    e2(α) = v2 * (1 .- α) / 2 + v4 * (1 .+ α) / 2
    e3(α) = v1 * (1 .- α) / 2 + v2 * (1 .+ α) / 2
    e4(α) = v3 * (1 .- α) / 2 + v4 * (1 .+ α) / 2
    e1α(α) = -v1 / 2 + v3 / 2
    e2α(α) = -v2 / 2 + v4 / 2
    e3α(α) = -v1 / 2 + v2 / 2
    e4α(α) = -v3 / 2 + v4 / 2
    transfinite_blend(e1, e2, e3, e4, e1α, e2α, e3α, e4α, r, s)
end

# connectivityarrays
function connectivityarrays(EToV, EToF)
    # number of elements
    nelems = size(EToV, 2)
    nfaces = maximum(maximum(EToF))

    # Determine secondary arrays
    # FToE : Unique Global Face to Element Number
    # FToLF: Unique Global Face to Element local face number
    # EToO : Element to Unique Global Faces Orientation
    # EToS : Element to Unique Global Face Side

    FToE = zeros(Int64, 2, nfaces)
    FToLF = zeros(Int64, 2, nfaces)
    EToO = Array{Bool, 2}(undef, 4, nelems)
    EToS = zeros(Int64, 4, nelems)

    # Local Face to Local Vertex map
    LFToLV = flatten_tuples(((1, 3), (2, 4), (1, 2), (3, 4)))
    for e in 1:nelems
        for lf in 1:4
            gf = EToF[lf, e]
            if FToE[1, gf] == 0
                @assert FToLF[1, gf] == 0
                FToE[1, gf] = e
                FToLF[1, gf] = lf
                EToO[lf, e] = true
                EToS[lf, e] = 1
            else
                @assert FToE[2, gf] == 0
                @assert FToLF[2, gf] == 0
                FToE[2, gf] = e
                FToLF[2, gf] = lf
                EToS[lf, e] = 2

                ne = FToE[1, gf]
                nf = FToLF[1, gf]

                nv = EToV[LFToLV[:, nf], ne]
                lv = EToV[LFToLV[:, lf], e]
                if nv == lv
                    EToO[lf, e] = true
                elseif nv[end:-1:1] == lv
                    EToO[lf, e] = false
                else
                    error("problem with connectivity")
                end
            end
        end
    end
    (FToE, FToLF, EToO, EToS)
end

function create_metrics(
    pm,
    Nr,
    Ns,
    xf = (r, s) -> (r, ones(size(r)), zeros(size(r))),
    yf = (r, s) -> (s, zeros(size(s)), ones(size(s))),
)
    Nrp = Nr + 1
    Nsp = Ns + 1
    Np = Nrp * Nsp

    # Derivative operators for the metric terms
    @assert pm <= 8
    pp = pm == 6 ? 8 : pm

    r = range(-1, stop = 1, length = Nrp)
    s = range(-1, stop = 1, length = Nsp)

    # Create the mesh
    r = ones(1, Nsp) ⊗ r
    s = s' ⊗ ones(Nrp)
    (x, xr, xs) = xf(r, s)
    (y, yr, ys) = yf(r, s)

    J = xr .* ys - xs .* yr
    @assert minimum(J) > 0

    JI = 1 ./ J

    rx = ys ./ J
    sx = -yr ./ J
    ry = -xs ./ J
    sy = xr ./ J

    # multiply by 2 since reference element is (-1, 1)
    hr = 2mapreduce(hypot, min, xr, yr) / Nr
    hs = 2mapreduce(hypot, min, xs, ys) / Ns
    hmin = min(hr, hs)

    # variable coefficient matrix components (BAE: these would only be so for
    # isotropic problem?)
    crr = J .* (rx .* rx + ry .* ry)
    crs = J .* (sx .* rx + sy .* ry)
    css = J .* (sx .* sx + sy .* sy)

    #
    # Block surface matrices
    #
    (xf1, yf1) = (view(x, 1, :), view(y, 1, :))
    nx1 = -ys[1, :]
    ny1 = xs[1, :]
    sJ1 = hypot.(nx1, ny1)
    nx1 = nx1 ./ sJ1
    ny1 = ny1 ./ sJ1

    (xf2, yf2) = (view(x, Nrp, :), view(y, Nrp, :))
    nx2 = ys[end, :]
    ny2 = -xs[end, :]
    sJ2 = hypot.(nx2, ny2)
    nx2 = nx2 ./ sJ2
    ny2 = ny2 ./ sJ2

    (xf3, yf3) = (view(x, :, 1), view(y, :, 1))
    nx3 = yr[:, 1]
    ny3 = -xr[:, 1]
    sJ3 = hypot.(nx3, ny3)
    nx3 = nx3 ./ sJ3
    ny3 = ny3 ./ sJ3

    (xf4, yf4) = (view(x, :, Nsp), view(y, :, Nsp))
    nx4 = -yr[:, end]
    ny4 = xr[:, end]
    sJ4 = hypot.(nx4, ny4)
    nx4 = nx4 ./ sJ4
    ny4 = ny4 ./ sJ4

    (
        coord = (x, y),
        facecoord = ((xf1, xf2, xf3, xf4), (yf1, yf2, yf3, yf4)),
        crr = crr,
        css = css,
        crs = crs,
        J = J,
        JI = JI,
        sJ = (sJ1, sJ2, sJ3, sJ4),
        nx = (nx1, nx2, nx3, nx4),
        ny = (ny1, ny2, ny3, ny4),
        rx = rx,
        ry = ry,
        sx = sx,
        sy = sy,
        hmin = hmin,
    )
end

function locoperator(
    p,
    Nr,
    Ns,
    metrics = create_metrics(p, Nr, Ns);
    crr = metrics.crr,
    css = metrics.css,
    crs = metrics.crs,
)
    csr = crs

    J = metrics.J
    JI = metrics.JI

    Nrp = Nr + 1
    Nsp = Ns + 1
    Np = Nrp * Nsp

    # Derivative operators for the rest of the computation
    (Dr, HrI, Hr, r) = DiagonalSBP.D1(p, Nr; xc = (-1, 1))
    Qr = Hr * Dr
    QrT = sparse(transpose(Qr))
    Dr = sparse(Dr)
    DrT = sparse(transpose(Dr))
    Hr = sparse(Hr)

    (Ds, HsI, Hs, s) = DiagonalSBP.D1(p, Ns; xc = (-1, 1))
    Qs = Hs * Ds
    QsT = sparse(transpose(Qs))
    Ds = sparse(Ds)
    DsT = sparse(transpose(Ds))
    Hs = sparse(Hs)

    # Identity matrices for the computation
    Ir = sparse(I, Nrp, Nrp)
    Is = sparse(I, Nsp, Nsp)

    # Set up the rr derivative matrix
    ISr0 = Array{Int64, 1}(undef, 0)
    JSr0 = Array{Int64, 1}(undef, 0)
    VSr0 = Array{Float64, 1}(undef, 0)
    ISrN = Array{Int64, 1}(undef, 0)
    JSrN = Array{Int64, 1}(undef, 0)
    VSrN = Array{Float64, 1}(undef, 0)

    Ae = DiagonalSBP.variable_D2(p, Nr, rand(Nrp)).A
    S0e = DiagonalSBP.variable_D2(p, Nr, ones(Nrp)).S0
    SNe = DiagonalSBP.variable_D2(p, Nr, ones(Nrp)).SN

    S0eT = sparse(transpose(S0e))
    SNeT = sparse(transpose(SNe))

    IArr = Array{Int64, 1}(undef, Nsp * length(Ae.nzval))
    JArr = Array{Int64, 1}(undef, Nsp * length(Ae.nzval))
    VArr = Array{Float64, 1}(undef, Nsp * length(Ae.nzval))
    stArr = 0

    for j in 1:Nsp
        rng = (j - 1) * Nrp .+ (1:Nrp)

        Ae = DiagonalSBP.variable_D2(p, Nr, crr[rng]).A
        (Ie, Je, Ve) = findnz(Ae)

        IArr[stArr .+ (1:length(Ve))] = Ie .+ (j - 1) * Nrp
        JArr[stArr .+ (1:length(Ve))] = Je .+ (j - 1) * Nrp
        VArr[stArr .+ (1:length(Ve))] = Hs[j, j] * Ve
        stArr += length(Ve)
    end

    Ãrr = sparse(IArr[1:stArr], JArr[1:stArr], VArr[1:stArr], Np, Np)

    # Set up the ss derivative matrix
    Ae = DiagonalSBP.variable_D2(p, Ns, rand(Nsp)).A

    IAss = Array{Int64, 1}(undef, Nrp * length(Ae.nzval))
    JAss = Array{Int64, 1}(undef, Nrp * length(Ae.nzval))
    VAss = Array{Float64, 1}(undef, Nrp * length(Ae.nzval))
    stAss = 0

    for i in 1:Nrp
        rng = i .+ Nrp * (0:Ns)
        Ae = DiagonalSBP.variable_D2(p, Ns, css[rng]).A
        R = Ae - Dr' * Hr * Diagonal(css[rng]) * Dr

        (Ie, Je, Ve) = findnz(Ae)
        IAss[stAss .+ (1:length(Ve))] = i .+ Nrp * (Ie .- 1)
        JAss[stAss .+ (1:length(Ve))] = i .+ Nrp * (Je .- 1)
        VAss[stAss .+ (1:length(Ve))] = Hr[i, i] * Ve
        stAss += length(Ve)
    end
    Ãss = sparse(IAss[1:stAss], JAss[1:stAss], VAss[1:stAss], Np, Np)

    # Set up the sr and rs derivative matrices
    Ãsr =
        (QsT ⊗ Ir) *
        sparse(1:length(crs), 1:length(crs), view(crs, :)) *
        (Is ⊗ Qr)
    Ãrs =
        (Is ⊗ QrT) *
        sparse(1:length(csr), 1:length(csr), view(csr, :)) *
        (Qs ⊗ Ir)

    Ã = Ãrr + Ãss + Ãrs + Ãsr

    #
    # Boundary point matrices
    #
    Er0 = sparse([1], [1], [1], Nrp, Nrp)
    ErN = sparse([Nrp], [Nrp], [1], Nrp, Nrp)
    Es0 = sparse([1], [1], [1], Nsp, Nsp)
    EsN = sparse([Nsp], [Nsp], [1], Nsp, Nsp)

    er0 = sparse([1], [1], [1], Nrp, 1)
    erN = sparse([Nrp], [1], [1], Nrp, 1)
    es0 = sparse([1], [1], [1], Nsp, 1)
    esN = sparse([Nsp], [1], [1], Nsp, 1)

    er0T = sparse([1], [1], [1], 1, Nrp)
    erNT = sparse([1], [Nrp], [1], 1, Nrp)
    es0T = sparse([1], [1], [1], 1, Nsp)
    esNT = sparse([1], [Nsp], [1], 1, Nsp)

    #
    # Store coefficient matrices on faces as matrices
    #TODO: switching 2 and 3 below??
    crr1 = sparse(Diagonal(crr[1 .+ Nrp * (0:Ns)]))
    crr3 = sparse(Diagonal(crr[1:Nrp]))
    crr2 = sparse(Diagonal(crr[Nrp .+ Nrp * (0:Ns)]))
    crr4 = sparse(Diagonal(crr[Nrp * Ns .+ (1:Nrp)]))

    crs1 = sparse(Diagonal(crs[1 .+ Nrp * (0:Ns)]))
    crs3 = sparse(Diagonal(crs[1:Nrp]))
    crs2 = sparse(Diagonal(crs[Nrp .+ Nrp * (0:Ns)]))
    crs4 = sparse(Diagonal(crs[Nrp * Ns .+ (1:Nrp)]))

    csr1 = sparse(Diagonal(csr[1 .+ Nrp * (0:Ns)]))
    csr3 = sparse(Diagonal(csr[1:Nrp]))
    csr2 = sparse(Diagonal(csr[Nrp .+ Nrp * (0:Ns)]))
    csr4 = sparse(Diagonal(csr[Nrp * Ns .+ (1:Nrp)]))

    css1 = sparse(Diagonal(css[1 .+ Nrp * (0:Ns)]))
    css3 = sparse(Diagonal(css[1:Nrp]))
    css2 = sparse(Diagonal(css[Nrp .+ Nrp * (0:Ns)]))
    css4 = sparse(Diagonal(css[Nrp * Ns .+ (1:Nrp)]))

    # Restriction operators: restrict the volume solution to the face
    L1 = Is ⊗ er0T
    L2 = Is ⊗ erNT
    L3 = es0T ⊗ Ir
    L4 = esNT ⊗ Ir

    #
    # Surface mass matrices
    #
    H1 = Hs
    H1I = HsI

    H2 = Hs
    H2I = HsI

    H3 = Hr
    H3I = HrI

    H4 = Hr
    H4I = HrI

    #
    # Penalty terms: α = θ_H
    #
    (l, θ_H, _, θ_R) = DiagonalSBP.D2_remainder_parameters(p)

    # check below
    #(θ_1, θ_2, θ_3, θ_4) = (H1[1,1], H2[Nsp, Nsp], H3[1,1], H4[Nrp, Nrp])
    #(ζ_1, ζ_2, ζ_3, ζ_4) = (β*H1[1,1], β*H2[Nsp, Nsp], β*H3[1,1], β*H4[Nrp, Nrp])

    #ψmin = reshape((crr + css - sqrt.((crr - css).^2 + 4crs.^2)) / 2, Nrp, Nsp)
    ψmin_r = reshape(crr, Nrp, Nsp)
    ψmin_s = reshape(css, Nrp, Nsp)
    @assert minimum(ψmin_r) > 0
    @assert minimum(ψmin_s) > 0

    hr = 2 / Nr
    hs = 2 / Ns

    ψ1 = ψmin_r[1, :]
    ψ2 = ψmin_r[Nrp, :]
    ψ3 = ψmin_s[:, 1]
    ψ4 = ψmin_s[:, Nsp]
    for k in 2:l
        ψ1 = min.(ψ1, ψmin_r[k, :])
        ψ2 = min.(ψ2, ψmin_r[Nrp + 1 - k, :])
        ψ3 = min.(ψ3, ψmin_s[:, k])
        ψ4 = min.(ψ4, ψmin_s[:, Nsp + 1 - k])
    end

    τR1 = (1 / (θ_R * hr)) * Is
    τR2 = (1 / (θ_R * hr)) * Is
    τR3 = (1 / (θ_R * hs)) * Ir
    τR4 = (1 / (θ_R * hs)) * Ir

    p1 = ((crr[1, :]) ./ ψ1)
    p2 = ((crr[Nrp, :]) ./ ψ2)
    p3 = ((css[:, 1]) ./ ψ3)
    p4 = ((css[:, Nsp]) ./ ψ4)

    P1 = sparse(1:Nsp, 1:Nsp, p1)
    P2 = sparse(1:Nsp, 1:Nsp, p2)
    P3 = sparse(1:Nrp, 1:Nrp, p3)
    P4 = sparse(1:Nrp, 1:Nrp, p4)
    # Penalty matrices:
    Γ1 = 1 * ((2 / (θ_H * hr)) * Is + τR1 * P1)
    Γ2 = 1 * ((2 / (θ_H * hr)) * Is + τR2 * P2)
    Γ3 = 1 * ((2 / (θ_H * hs)) * Ir + τR3 * P3)
    Γ4 = 1 * ((2 / (θ_H * hs)) * Ir + τR4 * P4)

    # Boundary matrices on faces
    Brr_1 = crr1 * L1 * (Is ⊗ S0e)
    Brr_2 = crr2 * L2 * (Is ⊗ SNe)
    Brr_3 = crr3 * L3 * (Is ⊗ Dr)
    Brr_4 = crr4 * L4 * (Is ⊗ Dr)
    Brr_1T = (Is ⊗ S0eT) * L1' * crr1
    Brr_2T = (Is ⊗ SNeT) * L2' * crr2
    Brr_3T = (Is ⊗ DrT) * L3' * crr3
    Brr_4T = (Is ⊗ DrT) * L4' * crr4

    Brs_1 = crs1 * L1 * (Ds ⊗ Ir) #TODO: START FIXING  HERE
    Brs_2 = crs2 * L2 * (Ds ⊗ Ir)
    Brs_3 = crs3 * L3 * (S0e ⊗ Ir)
    Brs_4 = crs4 * L4 * (SNe ⊗ Ir)
    Brs_1T = (DsT ⊗ Ir) * L1' * crs1
    Brs_2T = (DsT ⊗ Ir) * L2' * crs2
    Brs_3T = (S0eT ⊗ Ir) * L3' * crs3
    Brs_4T = (SNeT ⊗ Ir) * L4' * crs4

    Bsr_1 = csr1 * L1 * (Is ⊗ S0e)
    Bsr_2 = csr2 * L2 * (Is ⊗ SNe)
    Bsr_3 = csr3 * L3 * (Is ⊗ Dr)
    Bsr_4 = csr4 * L4 * (Is ⊗ Dr)
    Bsr_1T = (Is ⊗ S0eT) * L1' * csr1
    Bsr_2T = (Is ⊗ SNeT) * L2' * csr2
    Bsr_3T = (Is ⊗ DrT) * L3' * csr3
    Bsr_4T = (Is ⊗ DrT) * L4' * csr4

    Bss_1 = css1 * L1 * (Ds ⊗ Ir)
    Bss_2 = css2 * L2 * (Ds ⊗ Ir)
    Bss_3 = css3 * L3 * (S0e ⊗ Ir)
    Bss_4 = css4 * L4 * (SNe ⊗ Ir)
    Bss_1T = (DsT ⊗ Ir) * L1' * css1
    Bss_2T = (DsT ⊗ Ir) * L2' * css2
    Bss_3T = (S0eT ⊗ Ir) * L3' * css3
    Bss_4T = (SNeT ⊗ Ir) * L4' * css4

    HfI = (H1I, H2I, H3I, H4I)

    JH = sparse(1:Np, 1:Np, view(J, :)) * (Hs ⊗ Hr)
    JHI = sparse(1:Np, 1:Np, view(J, :)) * (HsI ⊗ HrI)
    J = sparse(1:Np, 1:Np, view(J, :))
    JI = sparse(1:Np, 1:Np, view(JI, :))
    HI = (HsI ⊗ HrI)
    H = (Hs ⊗ Hr)

    (
        Ã = Ã,
        coord = metrics.coord,
        facecoord = metrics.facecoord,
        JH = JH,
        JHI = JHI,
        sJ = metrics.sJ,
        nx = metrics.nx,
        ny = metrics.ny,
        Hf = (H1, H2, H3, H4),
        L = (L1, L2, L3, L4),
        Crrf = (crr1, crr2, crr3, crr4),
        Crsf = (crs1, crs2, crs3, crs4),
        Csrf = (csr1, csr2, csr3, csr4),
        Cssf = (css1, css2, css3, css4),
        Brr = (Brr_1, Brr_2, Brr_3, Brr_4, Brr_1T, Brr_2T, Brr_3T, Brr_4T),
        Brs = (Brs_1, Brs_2, Brs_3, Brs_4, Brs_1T, Brs_2T, Brs_3T, Brs_4T),
        Bsr = (Bsr_1, Bsr_2, Bsr_3, Bsr_4, Bsr_1T, Bsr_2T, Bsr_3T, Bsr_4T),
        Bss = (Bss_1, Bss_2, Bss_3, Bss_4, Bss_1T, Bss_2T, Bss_3T, Bss_4T),
        Γ = (Γ1, Γ2, Γ3, Γ4),
        HI = HI,
        H = H,
        J = J,
        JI = JI,
        HfI = (H1I, H2I, H3I, H4I),
    )
end

function rhsoperators(
    p,
    Nqr,
    Nqs,
    metrics,
    gDfun,
    gDdotfun,
    gNfun,
    F = nothing;
    rho = 1,
    crr = metrics.crr,
    css = metrics.css,
    crs = metrics.crs,
)

    # Build local operators

    lop = locoperator(
        p,
        Nqr - 1,
        Nqs - 1,
        metrics;
        crr = crr,
        css = css,
        crs = crs,
    )

    Np = Nqr * Nqs
    Ã = lop.Ã

    # Surface mass matrices
    H = lop.Hf

    (xf1, xf2, xf3, xf4) = metrics.facecoord[1]
    (yf1, yf2, yf3, yf4) = metrics.facecoord[2]

    #Surface normal vectors
    nx = metrics.nx
    ny = metrics.ny #DON'T assume normal vectors don't vary on face

    #hat normal vectors on 4 faces:
    nr = (-1, 1, 0, 0)
    ns = (0, 0, -1, 1)

    sJ = metrics.sJ

    # coefficients on the 4 faces
    crr = lop.Crrf
    crs = lop.Crsf
    csr = lop.Csrf
    css = lop.Cssf

    Z1 = sqrt.(rho .* (nx[1] .* nx[1] .+ ny[1] .* ny[1]))
    Z2 = sqrt.(rho .* (nx[2] .* nx[2] .+ ny[2] .* ny[2]))
    Z3 = sqrt.(rho .* (nx[3] .* nx[3] .+ ny[3] .* ny[3]))
    Z4 = sqrt.(rho .* (nx[4] .* nx[4] .+ ny[4] .* ny[4]))
    Z = (Z1, Z2, Z3, Z4)

    Ẑ1 = sJ[1] .* Z1
    Ẑ2 = sJ[2] .* Z2
    Ẑ3 = sJ[3] .* Z3
    Ẑ4 = sJ[4] .* Z4

    #Ẑ = (Ẑ1, Ẑ2, Ẑ3, Ẑ4)
    Ẑ = (
        spdiagm(0 => Ẑ1),
        spdiagm(0 => Ẑ2),
        spdiagm(0 => Ẑ3),
        spdiagm(0 => Ẑ4),
    )

    Ẑinv = (
        spdiagm(0 => 1 ./ Ẑ1),
        spdiagm(0 => 1 ./ Ẑ2),
        spdiagm(0 => 1 ./ Ẑ3),
        spdiagm(0 => 1 ./ Ẑ4),
    )

    # interpolation operators
    L = lop.L
    Γ = lop.Γ  # Need to fix these to be correct

    HI = lop.HI
    JH = lop.JH
    H̃ = lop.H
    JHI = lop.JHI
    H = lop.Hf
    J = lop.J
    JI = lop.JI
    JIHI = JI * HI

    # boundary deriv ops - include coefficients
    Brr = lop.Brr[1:4]
    BrrT = lop.Brr[5:8]
    Brs = lop.Brs[1:4]
    BrsT = lop.Brs[5:8]
    Bsr = lop.Bsr[1:4]
    BsrT = lop.Bsr[5:8]
    Bss = lop.Bss[1:4]
    BssT = lop.Bss[5:8]

    nCB = (
        nr[1] * (Brr[1] + Brs[1]) + ns[1] * (Bsr[1] + Bss[1]),
        nr[2] * (Brr[2] + Brs[2]) + ns[2] * (Bsr[2] + Bss[2]),
        nr[3] * (Brr[3] + Brs[3]) + ns[3] * (Bsr[3] + Bss[3]),
        nr[4] * (Brr[4] + Brs[4]) + ns[4] * (Bsr[4] + Bss[4]),
    )

    nCnΓ = (
        (
            nr[1] * crr[1] * nr[1] +
            nr[1] * crs[1] * ns[1] +
            ns[1] * csr[1] * nr[1] +
            ns[1] * css[1] * ns[1]
        ) * Γ[1],
        (
            nr[2] * crr[2] * nr[2] +
            nr[2] * crs[2] * ns[2] +
            ns[2] * csr[2] * nr[2] +
            ns[2] * css[2] * ns[2]
        ) * Γ[2],
        (
            nr[3] * crr[3] * nr[3] +
            nr[3] * crs[3] * ns[3] +
            ns[3] * csr[3] * nr[3] +
            ns[3] * css[3] * ns[3]
        ) * Γ[3],
        (
            nr[4] * crr[4] * nr[4] +
            nr[4] * crs[4] * ns[4] +
            ns[4] * csr[4] * nr[4] +
            ns[4] * css[4] * ns[4]
        ) * Γ[4],
    )

    nCnΓL = (
        (
            nr[1] * crr[1] * nr[1] +
            nr[1] * crs[1] * ns[1] +
            ns[1] * csr[1] * nr[1] +
            ns[1] * css[1] * ns[1]
        ) *
        Γ[1] *
        L[1],
        (
            nr[2] * crr[2] * nr[2] +
            nr[2] * crs[2] * ns[2] +
            ns[2] * csr[2] * nr[2] +
            ns[2] * css[2] * ns[2]
        ) *
        Γ[2] *
        L[2],
        (
            nr[3] * crr[3] * nr[3] +
            nr[3] * crs[3] * ns[3] +
            ns[3] * csr[3] * nr[3] +
            ns[3] * css[3] * ns[3]
        ) *
        Γ[3] *
        L[3],
        (
            nr[4] * crr[4] * nr[4] +
            nr[4] * crs[4] * ns[4] +
            ns[4] * csr[4] * nr[4] +
            ns[4] * css[4] * ns[4]
        ) *
        Γ[4] *
        L[4],
    )

    BtnCH = (
        (nr[1] * (BrrT[1] + BrsT[1]) + ns[1] * (BsrT[1] + BssT[1])) * H[1],
        (nr[2] * (BrrT[2] + BrsT[2]) + ns[2] * (BsrT[2] + BssT[2])) * H[2],
        (nr[3] * (BrrT[3] + BrsT[3]) + ns[3] * (BsrT[3] + BssT[3])) * H[3],
        (nr[4] * (BrrT[4] + BrsT[4]) + ns[4] * (BsrT[4] + BssT[4])) * H[4],
    )

    X1 =
        spdiagm(0 => 1 ./ diag(Γ[1])) *
        ((1 / nr[1]) * spdiagm(0 => 1 ./ diag(crr[1])) * (1 / nr[1]))
    X2 =
        spdiagm(0 => 1 ./ diag(Γ[2])) *
        ((1 / nr[2]) * spdiagm(0 => 1 ./ diag(crr[2])) * (1 / nr[2]))
    X3 =
        spdiagm(0 => 1 ./ diag(Γ[3])) *
        ((1 / ns[3]) * spdiagm(0 => 1 ./ diag(css[3])) * (1 / ns[3]))
    X4 =
        spdiagm(0 => 1 ./ diag(Γ[4])) *
        ((1 / ns[4]) * spdiagm(0 => 1 ./ diag(css[4])) * (1 / ns[4]))

    X = (X1, X2, X3, X4)

    ue_n =
        (t, e) -> (
            gNfun.(nx[1], ny[1], xf1[:], yf1[:], t, e),
            gNfun.(nx[2], ny[2], xf2[:], yf2[:], t, e),
            gNfun.(nx[3], ny[3], xf3[:], yf3[:], t, e),
            gNfun.(nx[4], ny[4], xf4[:], yf4[:], t, e),
        )

    gD =
        (t, e) -> (
            gDfun.(xf1[:], yf1[:], t, e),
            gDfun.(xf2[:], yf2[:], t, e),
            gDfun.(xf3[:], yf3[:], t, e),
            gDfun.(xf4[:], yf4[:], t, e),
        )
    gDdot =
        (t, e) -> (
            gDdotfun.(xf1[:], yf1[:], t, e),
            gDdotfun.(xf2[:], yf2[:], t, e),
            gDdotfun.(xf3[:], yf3[:], t, e),
            gDdotfun.(xf4[:], yf4[:], t, e),
        )

    sJgDdot =
        (t, e) -> (
            spdiagm(0 => sJ[1]) * gDdot(t, e)[1],
            spdiagm(0 => sJ[2]) * gDdot(t, e)[2],
            spdiagm(0 => sJ[3]) * gDdot(t, e)[3],
            spdiagm(0 => sJ[4]) * gDdot(t, e)[4],
        )

    gN =
        (t, e) -> (
            spdiagm(0 => sJ[1]) * ue_n(t, e)[1],
            spdiagm(0 => sJ[2]) * ue_n(t, e)[2],
            spdiagm(0 => sJ[3]) * ue_n(t, e)[3],
            spdiagm(0 => sJ[4]) * ue_n(t, e)[4],
        )

    if isnothing(F)
        myF = nothing
    else
        myF = (t, e) -> F.(metrics.coord[1][:], metrics.coord[2][:], t, e)
    end
    (
        Ã = Ã,
        nCB = nCB,
        nCnΓ = nCnΓ,
        nCnΓL = nCnΓL,
        gD = gD,
        gDdot = gDdot,
        Z = Ẑ,
        L = L,
        gN = gN,
        sJgDdot = sJgDdot,
        Zinv = Ẑinv,
        H = H,
        BtnCH = BtnCH,
        JIHI = JIHI,
        F = myF,
        JH = JH,
        sJ = sJ,
        X = X,
        H̃ = H̃,
        J = J,
    )
end

# Constructor for inp files
function read_inp_2d(T, S, filename::String; bc_map = 1:10000)
    # Read in the file
    f = try
        open(filename)
    catch
        error("InpRead cannot open \"$filename\" ")
    end
    lines = readlines(f)
    close(f)

    # Read in nodes
    str = "NSET=ALLNODES"
    linenum = SeekToSubstring(lines, str)
    linenum > 0 || error("did not find: $str")
    num_nodes = 0
    for l in (linenum + 1):length(lines)
        occursin(r"^\s*[0-9]*\s*,.*", lines[l]) ? num_nodes += 1 : break
    end
    Vx = fill(S(NaN), num_nodes)
    Vy = fill(S(NaN), num_nodes)
    Vz = fill(S(NaN), num_nodes)
    for l in linenum .+ (1:num_nodes)
        node_data = split(lines[l], r"\s|,", keepempty = false)
        (node_num, node_x, node_y, node_z) = try
            (
                parse(T, node_data[1]),
                parse(S, node_data[2]),
                parse(S, node_data[3]),
                parse(S, node_data[4]),
            )
        catch
            error("cannot parse line $l: \"$(lines[l])\" ")
        end

        Vx[node_num] = node_x
        Vy[node_num] = node_y
        Vz[node_num] = node_z
    end

    # Read in Elements
    str = "ELEMENT"
    linenum = SeekToSubstring(lines, str)
    num_elm = 0
    while linenum > 0
        for l in linenum .+ (1:length(lines))
            occursin(r"^\s*[0-9]*\s*,.*", lines[l]) ? num_elm += 1 : break
        end
        linenum = SeekToSubstring(lines, str; first = linenum + 1)
    end
    num_elm > 0 || error("did not find any element")

    EToV = fill(T(0), 4, num_elm)
    EToBlock = fill(T(0), num_elm)
    linenum = SeekToSubstring(lines, str)
    while linenum > 0
        foo = split(lines[linenum], r"[^0-9]", keepempty = false)
        B = parse(T, foo[end])
        for l in linenum .+ (1:num_elm)
            elm_data = split(lines[l], r"\s|,", keepempty = false)
            # read into z-order
            (elm_num, elm_v1, elm_v2, elm_v4, elm_v3) = try
                (
                    parse(T, elm_data[1]),
                    parse(T, elm_data[2]),
                    parse(T, elm_data[3]),
                    parse(T, elm_data[4]),
                    parse(T, elm_data[5]),
                )
            catch
                break
            end
            EToV[:, elm_num] = [elm_v1, elm_v2, elm_v3, elm_v4]
            EToBlock[elm_num] = B
        end
        linenum = SeekToSubstring(lines, str; first = linenum + 1)
    end

    # Determine connectivity
    EToF = fill(T(0), 4, num_elm)

    VsToF = Dict{Tuple{Int64, Int64}, Int64}()
    numfaces = 0
    for e in 1:num_elm
        for lf in 1:4
            if lf == 1
                Vs = (EToV[1, e], EToV[3, e])
            elseif lf == 2
                Vs = (EToV[2, e], EToV[4, e])
            elseif lf == 3
                Vs = (EToV[1, e], EToV[2, e])
            elseif lf == 4
                Vs = (EToV[3, e], EToV[4, e])
            end
            if Vs[1] > Vs[2]
                Vs = (Vs[2], Vs[1])
            end
            if haskey(VsToF, Vs)
                EToF[lf, e] = VsToF[Vs]
            else
                numfaces = numfaces + 1
                EToF[lf, e] = VsToF[Vs] = numfaces
            end
        end
    end

    # Read in side set info
    FToB = Array{T, 1}(undef, numfaces)
    fill!(FToB, BC_LOCKED_INTERFACE)
    linenum = SeekToSubstring(lines, "\\*ELSET")
    inp_to_zorder = [3, 2, 4, 1]
    while linenum > 0
        foo = split(lines[linenum], r"[^0-9]", keepempty = false)
        (bc, face) = try
            (parse(T, foo[1]), parse(T, foo[2]))
        catch
            error("cannot parse line $linenum: \"$(lines[linenum])\" ")
        end
        bc = bc_map[bc]
        face = inp_to_zorder[face]
        for l in (linenum + 1):length(lines)
            if !occursin(r"^\s*[0-9]+", lines[l])
                break
            end
            elms = split(lines[l], r"\s|,", keepempty = false)
            for elm in elms
                elm = try
                    parse(T, elm)
                catch
                    error("cannot parse line $linenum: \"$(lines[l])\" ")
                end
                if bc == 3
                    bc = BC_LOCKED_INTERFACE
                end
                FToB[EToF[face, elm]] = bc
                @assert (
                    bc == BC_DIRICHLET ||
                    bc == BC_NEUMANN ||
                    bc == BC_LOCKED_INTERFACE ||
                    abs(bc) >= BC_JUMP_INTERFACE
                )
            end
        end
        linenum = SeekToSubstring(lines, "\\*ELSET"; first = linenum + 1)
    end

    ([Vx Vy]', EToV, EToF, FToB, EToBlock)
end
read_inp_2d(filename; kw...) = read_inp_2d(Int64, Float64, filename; kw...)

function SeekToSubstring(lines, substring; first = 1)
    for l in first:length(lines)
        if occursin(Regex(".*$(substring).*"), lines[l])
            return l
        end
    end
    return -1
end

function plot_connectivity(verts, EToV)
    Lx = extrema(verts[1, :])
    Lx = (floor(Int, Lx[1]), ceil(Int, Lx[2]))
    Ly = extrema(verts[2, :])
    Ly = (floor(Int, Ly[1]), ceil(Int, Ly[2]))
    width = Lx[2] - Lx[1]
    height = Ly[2] - Ly[1]
    plt = Plot(
        BrailleCanvas(
            80,
            ceil(Int, 40 * height / width),
            origin_x = Lx[1],
            origin_y = Ly[1],
            width = width,
            height = height,
        ),
    )

    annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
    annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
    annotate!(plt, :bl, string(Lx[1]), color = :light_black)
    annotate!(plt, :br, string(Lx[2]), color = :light_black)
    for e in 1:size(EToV, 2)
        (v1, v2, v3, v4) = EToV[1:4, e]
        x = verts[1, [v1 v2 v4 v3 v1]][:]
        y = verts[2, [v1 v2 v4 v3 v1]][:]
        lineplot!(plt, x, y)
    end
    title!(plt, "connectivity")
    display(plt)
end

function plot_blocks(lop, FToB, EToF)
    Lx = (floatmax(), -floatmax())
    Ly = (floatmax(), -floatmax())
    for e in 1:length(lop)
        (x, y) = lop[e].coord
        Lxe = extrema(x)
        Lye = extrema(y)
        Lx = (min(Lx[1], Lxe[1]), max(Lx[2], Lxe[2]))
        Ly = (min(Ly[1], Lye[1]), max(Ly[2], Lye[2]))
    end

    Lx = (floor(Int, Lx[1]), ceil(Int, Lx[2]))
    Ly = (floor(Int, Ly[1]), ceil(Int, Ly[2]))

    width = Lx[2] - Lx[1]
    height = Ly[2] - Ly[1]
    plt = Plot(
        BrailleCanvas(
            80,
            ceil(Int, 40 * height / width),
            origin_x = Lx[1],
            origin_y = Ly[1],
            width = width,
            height = height,
        ),
    )

    annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
    annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
    annotate!(plt, :bl, string(Lx[1]), color = :light_black)
    annotate!(plt, :br, string(Lx[2]), color = :light_black)

    for e in 1:length(lop)
        (xf, yf) = lop[e].facecoord
        for lf in 1:length(xf)
            bc_type_face = FToB[EToF[lf, e]]
            if bc_type_face == BC_LOCKED_INTERFACE
                lineplot!(plt, xf[lf], yf[lf], color = :blue)
            elseif bc_type_face == BC_DIRICHLET
                lineplot!(plt, xf[lf], yf[lf], color = :green)
            elseif bc_type_face == BC_DIRICHLET
                lineplot!(plt, xf[lf], yf[lf], color = :yellow)
            else
                lineplot!(plt, xf[lf], yf[lf], color = :red)
            end
        end
    end
    title!(plt, "mesh")
    display(plt)
end

function write_vtk(
    output_file,
    metrics,
    q;
    cxx = nothing,
    cyy = nothing,
    cxy = nothing,
)
    vtk_multiblock(output_file) do vtm
        # How many blocks do we have
        nblocks = length(metrics)

        # loop over the blocks and output for this block

        # pointer for end of the last block
        n = 0
        for block in 1:nblocks
            # get the grid
            x = metrics[block].coord[1]
            y = metrics[block].coord[2]

            # get the number of points associated with this block
            Nr, Ns = size(x)
            Np = Nr * Ns

            # get the length of the q vector for this block
            q_len = 2(Np + Nr + Ns)

            # get the displacements

            vtk = vtk_grid(vtm, x, y)
            vtk["u"] = @view q[n .+ (1:Np)]
            vtk["v"] = @view q[(n + Np) .+ (1:Np)]

            if !isnothing(cxx)
                vtk["cxx"] = cxx.(x, y)
                vtk["cxy"] = cxy.(x, y)
                vtk["cyy"] = cyy.(x, y)
            end

            # update the last point of the block
            n += q_len
        end

        # We should have hit all the data by now
        @assert n == length(q)
    end
end

# This uses the fourth-order, low-storage, Runge--Kutta scheme of
# Carpenter and Kennedy (1994) (in their notation (5,4) 2N-Storage RK
#                               scheme).
#
# @TECHREPORT{CarpenterKennedy1994,
#             author = {M.~H. Carpenter and C.~A. Kennedy},
#             title = {Fourth-order {2N-storage} {Runge-Kutta}
#                      schemes},
#             institution = {National Aeronautics and Space
#                            Administration},
#             year = {1994},
#             number = {NASA TM-109112},
#             address = {Langley Research
#                        Center, Hampton, VA},
#            }
function timestep!(q, f!, p, dt, (t0, t1), Δq = similar(q), Δq2 = similar(q))
    T = eltype(q)

    RKA = (
        T(0),
        T(-567301805773 // 1357537059087),
        T(-2404267990393 // 2016746695238),
        T(-3550918686646 // 2091501179385),
        T(-1275806237668 // 842570457699),
    )

    RKB = (
        T(1432997174477 // 9575080441755),
        T(5161836677717 // 13612068292357),
        T(1720146321549 // 2090206949498),
        T(3134564353537 // 4481467310338),
        T(2277821191437 // 14882151754819),
    )

    RKC = (
        T(0),
        T(1432997174477 // 9575080441755),
        T(2526269341429 // 6820363962896),
        T(2006345519317 // 3224310063776),
        T(2802321613138 // 2924317926251),
    )

    nstep = ceil(Int, (t1 - t0) / dt)
    dt = (t1 - t0) / nstep

    fill!(Δq, 0)
    fill!(Δq2, 0)
    for step in 1:nstep
        t = t0 + (step - 1) * dt
        for s in 1:length(RKA)
            f!(Δq2, q, p, t + RKC[s] * dt)
            Δq .+= Δq2
            q .+= RKB[s] * dt * Δq
            Δq .*= RKA[s % length(RKA) + 1]
        end
    end

    nothing
end

# Based on the zero algorithm of  Richard P.  Brent, Algorithms for Minimization
# without Derivatives, Prentice-Hall, Englewood Cliffs, New Jersey, 1973, 195
# pp. (available online at
# https://maths-people.anu.edu.au/~brent/pub/pub011.html)
function brentdekker(
    f,
    a::T,
    b::T,
    xatol::T = eps(T),
    max_iterations = 100,
    xrtol::T = xatol,
    ftol::T = T(0),
) where {T <: AbstractFloat}

    # Evaluate the function
    fa, fb = f(a), f(b)

    # Not a bracket
    if fa * fb > 0
        return (b, false, -1)
    end

    # Initialize the other side of the bracket
    c, fc = a, fa

    # below d is the current step and e is previous step
    d = e = b - a

    macheps = eps(T)

    # solution is in bracket [c, b] and a is the previous b
    for i in 1:max_iterations

        # `b` should be best guess (even if that means we don't update value)
        # This also resets `a` (the previous `b`) to be the same as `c
        if abs(fc) < abs(fb)
            a, b, c = b, c, b
            fa, fb, fc = fb, fc, fb
        end

        # tol is used to make sure that we are properly scaled with respect to
        # the input
        tol = 2xrtol * abs(b) + xatol

        # bisection step size
        m = (c - b) / 2

        # If the bracket is small enough or we hit the root
        if (abs(m) < tol || fb == 0 || abs(fb) < ftol)
            return (b, true, i)
        end

        # if the step we took last time `e` is smaller than the tolerance OR if
        # our function really increased just do bisection
        if abs(e) < tol || abs(fa) < abs(fb)
            d = e = m
        else
            # do inverse quadratic interpolation if we can
            # otherwise do secant
            # The use of p & q is to make things more floating point stable

            r = fb / fc
            if a ≠ c
                s = fb / fa
                q = fa / fc
                p = s * (2m * q * (q - r) - (b - a) * (r - 1))
                q = (q - 1) * (r - 1) * (s - 1)
            else
                p = 2m * r
                q = 1 - r
            end

            # either p or q has to flip signs, since we want p to be positive
            # below we prefer to flip the sign of q
            p > 0 ? (q = -q) : (p = -p)

            # As long as the step keeps us in the bracket and the step is not too
            # large we accept the secant / quadratic step
            #
            # The first condition ensures that the step p / q cause the solution
            # to be in the interval [b, b + (3/4)*(a-b)], e.g, we do not get too
            # close to a which is the worst of the two sides.
            #
            # The second condition ensures that the step is at least half as big
            # as the previous step, e.g., p / q < e / 2 (otherwise we are not
            # converging quickly so we should just bisect)
            if 2p < 3m * q - abs(tol * q) && p < abs(e * q / 2)
                e, d = d, p / q
            else
                d = e = m
            end
        end

        # Save the last step
        a, fa = b, fb

        # As long as the step isn't too small accept it, otherwise take a `tol`
        # sizes step in the correct direction
        b = b + (abs(d) > tol ? d : (m > 0 ? tol : -tol))
        fb = f(b)

        # if fb and fc have same sign then really a, b bracket the root
        if (fb > 0) == (fc > 0)
            c, fc = a, fa
            d = e = b - a
        end
    end

    return (b, false, max_iterations)
end

function waveprop!(dq, q, params, t)
    Nqr = params.Nqr
    Nqs = params.Nqs
    rhsops = params.rhsops
    EToF = params.EToF
    EToO = params.EToO
    EToS = params.EToS
    FToB = params.FToB
    FToE = params.FToE
    FToLF = params.FToLF
    friction = params.friction

    nelems = length(rhsops)
    Np = Nqr * Nqs

    lenq0 = 2 * (Np + Nqr + Nqs)

    # @show t

    @assert length(q) == nelems * lenq0
    @inbounds Threads.@threads for e in 1:nelems
        τ̂star = (
            zeros(Float64, Nqs),
            zeros(Float64, Nqs),
            zeros(Float64, Nqr),
            zeros(Float64, Nqr),
        )

        vstar = (
            zeros(Float64, Nqs),
            zeros(Float64, Nqs),
            zeros(Float64, Nqr),
            zeros(Float64, Nqr),
        )

        qe = @view q[(lenq0 * (e - 1) + 1):(e * lenq0)]
        dqe = @view dq[(lenq0 * (e - 1) + 1):(e * lenq0)]

        u = @view qe[1:Np]
        v = @view qe[Np .+ (1:Np)]
        û1 = @view qe[(2Np + 1):(2Np + Nqs)]
        û2 = @view qe[(2Np + 1 + Nqs):(2(Np + Nqs))]
        û3 = @view qe[(2(Np + Nqs) + 1):(2(Np + Nqs) + Nqr)]
        û4 = @view qe[(2(Np + Nqs) + Nqr + 1):(2(Np + Nqs + Nqr))]

        ustar = (û1, û2, û3, û4)

        du = @view dqe[1:Np]
        dv = @view dqe[Np .+ (1:Np)]
        dû1 = @view dqe[(2Np + 1):(2Np + Nqs)]
        dû2 = @view dqe[(2Np + 1 + Nqs):(2(Np + Nqs))]
        dû3 = @view dqe[(2(Np + Nqs) + 1):(2(Np + Nqs) + Nqr)]
        dû4 = @view dqe[(2(Np + Nqs) + Nqr + 1):(2(Np + Nqs + Nqr))]

        #vstar = (dû1, dû2, dû3, dû4)

        τ̂ = (
            rhsops[e].nCB[1] * u + rhsops[e].nCnΓ[1] * û1 -
            rhsops[e].nCnΓL[1] * u,
            rhsops[e].nCB[2] * u + rhsops[e].nCnΓ[2] * û2 -
            rhsops[e].nCnΓL[2] * u,
            rhsops[e].nCB[3] * u + rhsops[e].nCnΓ[3] * û3 -
            rhsops[e].nCnΓL[3] * u,
            rhsops[e].nCB[4] * u + rhsops[e].nCnΓ[4] * û4 -
            rhsops[e].nCnΓL[4] * u,
        )

        glb_fcs = EToF[:, e]  #global faces of element e

        for lf in 1:4
            bc_type_face = FToB[glb_fcs[lf]]
            if bc_type_face == BC_DIRICHLET
                τ̂star[lf] .= τ̂[lf]
                vstar[lf] .= rhsops[e].gDdot(t, e)[lf]

            elseif bc_type_face == BC_NEUMANN
                τ̂star[lf] .= rhsops[e].gN(t, e)[lf]
                vstar[lf] .= rhsops[e].L[lf] * v

            elseif bc_type_face == BC_JUMP_INTERFACE
                els_share = FToE[:, glb_fcs[lf]]
                lf_share_glb = FToLF[:, glb_fcs[lf]]
                cel = [1]
                cfc = [1]
                for i in 1:2
                    if els_share[i] != e
                        cel[1] = els_share[i] # other element face is connected to!
                    end
                end

                if e < cel[1]
                    cfc[1] = lf_share_glb[2]
                else
                    cfc[1] = lf_share_glb[1]
                end

                qplus = @view q[(lenq0 * (cel[1] - 1) + 1):(cel[1] * lenq0)]
                uplus = @view qplus[1:Np]
                vplus = @view qplus[Np .+ (1:Np)]
                û1plus = qplus[(2Np + 1):(2Np + Nqs)]
                û2plus = qplus[(2Np + 1 + Nqs):(2(Np + Nqs))]
                û3plus = qplus[(2(Np + Nqs) + 1):(2(Np + Nqs) + Nqr)]
                û4plus = qplus[(2(Np + Nqs) + Nqr + 1):(2(Np + Nqs + Nqr))]

                ustarplus = (û1plus, û2plus, û3plus, û4plus)
                τ̂plus =
                    rhsops[cel[1]].nCB[cfc[1]] * uplus +
                    rhsops[cel[1]].nCnΓ[cfc[1]] * ustarplus[cfc[1]] -
                    rhsops[cel[1]].nCnΓL[cfc[1]] * uplus

                revrse_e = EToO[lf, e]
                revrse_cel = EToO[cfc, cel]

                vplus_fc = rhsops[cel[1]].L[cfc[1]] * vplus
                Zplus = rhsops[cel[1]].Z[cfc[1]]

                gDdotplus = rhsops[cel[1]].gDdot(t, cel[1])[cfc[1]]
                gNplus = rhsops[cel[1]].gN(t, cel[1])[cfc[1]]
                sJplus = rhsops[cel[1]].sJ[cfc[1]]
                sJminus = rhsops[e].sJ[lf]

                if EToO[lf, e] != EToO[cfc[1], cel[1]]
                    τ̂plus_rev = τ̂plus[end:-1:1]
                    vplus_fc_rev = vplus_fc[end:-1:1] #TODO: could change these and above to avoid defining new arrays
                    Z_rev = Zplus[end:-1:1, end:-1:1]
                    gNplus_rev = gNplus[end:-1:1]
                    gDdotplus_rev = gDdotplus[end:-1:1]
                    sJplus_rev = sJplus[end:-1:1]
                else
                    τ̂plus_rev = τ̂plus
                    vplus_fc_rev = vplus_fc
                    Z_rev = Zplus
                    gNplus_rev = gNplus
                    gDdotplus_rev = gDdotplus
                    sJplus_rev = sJplus
                end

                # Z_refI = spdiagm(0 => 1 ./ diag(Z_rev))
                # ZI = rhsops[e].Zinv[lf]
                Z = rhsops[e].Z[lf]

                Vminus = gDdotplus_rev - rhsops[e].gDdot(t, e)[lf]
                sJgminus =
                    rhsops[e].gN(t, e)[lf] -
                    rhsops[e].sJ[lf] .* friction.(Vminus)
                sJgplus = gNplus_rev - sJplus_rev .* friction.(-Vminus)

                #gminus = rhsops[e].Zinv[lf] * sJgminus
                #gplus = Z_refI * sJgplus

                Zvminus = Z * rhsops[e].L[lf] * v
                Zvplus = Z_rev * vplus_fc_rev

                vfric = zeros(Float64, Nqr)

                for i in 1:Nqr
                    f(V) =
                        2 * sJminus[i] * friction(V) + Z[i, i] * V -
                        (Zvplus[i] - Zvminus[i]) + τ̂plus_rev[i] - τ̂[lf][i] -
                        sJgplus[i] + sJgminus[i]

                    Left =
                        (1 / Z[i, i]) * (
                            sJgplus[i] - sJgminus[i] - τ̂plus_rev[i] +
                            τ̂[lf][i] +
                            Zvplus[i] - Zvminus[i]
                        )

                    Right = -Left

                    if Left >= Right
                        tmp = Left
                        Left = Right
                        Right = tmp
                    end

                    v_root = brentdekker(f, Left, Right)#, xatol = 1e-9, xrtol = 1e-9, ftol = 1e-9)
                    vfric[i] = v_root[1]
                end

                vstar[lf] .=
                    rhsops[e].L[lf] * v +
                    rhsops[e].Zinv[lf] * (
                        -1 * τ̂[lf] +
                        rhsops[e].sJ[lf] .* friction.(vfric) +
                        sJgminus
                    )
                τ̂star[lf] .=
                    rhsops[e].Z[lf] * (vstar[lf] - rhsops[e].L[lf] * v) + τ̂[lf]

            elseif bc_type_face == -BC_JUMP_INTERFACE
                els_share = FToE[:, glb_fcs[lf]]
                lf_share_glb = FToLF[:, glb_fcs[lf]]
                cel = [1]
                cfc = [1]
                for i in 1:2
                    if els_share[i] != e
                        cel[1] = els_share[i] # other element face is connected to!
                    end
                end
                gDdotplus = rhsops[cel[1]].gDdot(t, cel[1])[cfc[1]]
                qplus = @view q[(lenq0 * (cel[1] - 1) + 1):(cel[1] * lenq0)]
                vplus = @view qplus[Np .+ (1:Np)]
                vplus = rhsops[cel[1]].L[cfc[1]] * vplus
                if EToO[lf, e] != EToO[cfc[1], cel[1]]
                    gDdotplus .= gDdotplus[end:-1:1]
                    vplus .= vplus[end:-1:1]
                end
                Vdisc = vplus - rhsops[e].L[lf] * v
                Vexact = gDdotplus - rhsops[e].gDdot(t, e)[lf]
                τ̂star[lf] .=
                    rhsops[e].gN(t, e)[lf] -
                    rhsops[e].sJ[lf] .* (friction.(Vexact) - friction.(Vdisc))
                vstar[lf] .= rhsops[e].L[lf] * v

            elseif bc_type_face == BC_LOCKED_INTERFACE
                els_share = FToE[:, glb_fcs[lf]]
                lf_share_glb = FToLF[:, glb_fcs[lf]]
                cel = [1]
                cfc = [1]
                for i in 1:2
                    if els_share[i] != e
                        cel[1] = els_share[i] # other element face is connected to!
                    end
                end

                if e < cel[1]
                    cfc[1] = lf_share_glb[2]
                else
                    cfc[1] = lf_share_glb[1]
                end

                qplus = @view q[(lenq0 * (cel[1] - 1) + 1):(cel[1] * lenq0)]
                uplus = @view qplus[1:Np]
                vplus = @view qplus[Np .+ (1:Np)]
                û1plus = qplus[(2Np + 1):(2Np + Nqs)]
                û2plus = qplus[(2Np + 1 + Nqs):(2(Np + Nqs))]
                û3plus = qplus[(2(Np + Nqs) + 1):(2(Np + Nqs) + Nqr)]
                û4plus = qplus[(2(Np + Nqs) + Nqr + 1):(2(Np + Nqs + Nqr))]

                ustarplus = (û1plus, û2plus, û3plus, û4plus)
                τ̂plus =
                    rhsops[cel[1]].nCB[cfc[1]] * uplus +
                    rhsops[cel[1]].nCnΓ[cfc[1]] * ustarplus[cfc[1]] -
                    rhsops[cel[1]].nCnΓL[cfc[1]] * uplus

                revrse_e = EToO[lf, e]
                revrse_cel = EToO[cfc, cel]
                vplus_fc = rhsops[cel[1]].L[cfc[1]] * vplus

                if EToO[lf, e] != EToO[cfc[1], cel[1]]#e == 1 || e == 2
                    τ̂plus_rev = τ̂plus[end:-1:1]
                    vplus_fc_rev = vplus_fc[end:-1:1]
                else
                    τ̂plus_rev = τ̂plus
                    vplus_fc_rev = vplus_fc
                end

                τ̂star[lf] .= 0.5 * (τ̂[lf] - τ̂plus_rev)
                vstar[lf] .= 0.5 * (rhsops[e].L[lf] * v + vplus_fc_rev)
            else
                print("illegal bc type")
                return
            end
        end

        du .= v
        # dv .= -rhsops[e].Ã * u
        mul!(dv, rhsops[e].Ã, u, -1, 0)

        for lf in 1:4
            dv .+=
                rhsops[e].L[lf]' * (rhsops[e].H[lf] * τ̂star[lf]) -
                rhsops[e].BtnCH[lf] * (ustar[lf] - rhsops[e].L[lf] * u)
        end

        if isnothing(rhsops[e].F)
            dv .= rhsops[e].JIHI * dv
        else
            dv .= rhsops[e].JIHI * dv .+ rhsops[e].F(t, e)
        end

        dû1 .= vstar[1]
        dû2 .= vstar[2]
        dû3 .= vstar[3]
        dû4 .= vstar[4]
    end
end


text:
using Profile
using ProfileView
using GtkObservables
using AbstractTrees
using Test

function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A, dims = 2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B, dims = 1)
        b = rand(100)
        C = B.*b
    end
end

unstable(x) = x > 0.5 ? true : 0.0

function profile_unstable_test(m, n)
    s = s2 = 0
    for i = 1:n
        for k = 1:m
            s += unstable(rand())
        end
        x = collect(1:20)
        s2 += sum(x)
    end
    s, s2
end

function add2(x)
    y = x[1] + x[2]
    return y, backtrace()
end

@testset "ProfileView" begin
    @testset "windows" begin
        profile_test(1)
        @test isa(@profview(profile_test(10)), ProfileView.Window)
        data, lidict = Profile.retrieve()

        Profile.clear()
        @profile profile_test(10)
        @test isa(ProfileView.view(), ProfileView.Window)
        @test isa(ProfileView.view(C=true), ProfileView.Window)
        @test isa(ProfileView.view(fontsize=18), ProfileView.Window)
        @test isa(ProfileView.view(windowname="ProfileWindow"), ProfileView.Window)

        Profile.clear()
        profile_unstable_test(1, 1)
        @profile profile_unstable_test(10, 10^6)
        @test isa(ProfileView.view(), ProfileView.Window)

        @test isa(ProfileView.view(ProfileView.FlameGraphs.flamegraph()), ProfileView.Window)
        @test isa(ProfileView.view(ProfileView.FlameGraphs.FlameColors()), ProfileView.Window)

        data, lidict = Profile.retrieve()
        @test isa(ProfileView.view(data, lidict=lidict), ProfileView.Window)

        @test isa(ProfileView.view(nothing), ProfileView.Window)

        # Interactivity
        stackframe(func, file, line; C=false) = ProfileView.StackFrame(Symbol(func), Symbol(file), line, nothing, C, false, 0)

        backtraces = UInt64[0, 4, 3, 2, 1,   # order: calles then caller
                            0, 6, 5, 1,
                            0, 8, 7,
                            0, 4, 3, 2, 1,
                            0]
        if isdefined(Profile, :add_fake_meta)
            backtraces = Profile.add_fake_meta(backtraces)
        end
        lidict = Dict{UInt64,Vector{ProfileView.StackFrame}}(
            1=>[stackframe(:f1, :file1, 1)],
            2=>[stackframe(:f2, :file1, 5)],
            3=>[stackframe(:f3, :file2, 1)],
            4=>[stackframe(:f2, :file1, 15)],
            5=>[stackframe(:f4, :file1, 20)],
            6=>[stackframe(:f5, :file3, 1)],
            7=>[stackframe(:f1, :file1, 2)],
            8=>[stackframe(:f6, :file3, 10)])
        g = ProfileView.flamegraph(backtraces; lidict=lidict)
        win, c, fdraw, (tb_open, tb_save_as) = ProfileView.viewgui(ProfileView.FlameGraphs.default_colors, g);
        ProfileView.Gtk.showall(win)
        sleep(1.0)
        @test c.widget.is_realized && c.widget.is_sized  # to ensure the motion test really runs
        btn = c.mouse.motion[]
        c.mouse.motion[] = MouseButton(XY{UserUnit}(2.8, 1.4), btn.button, btn.clicktype, btn.modifiers)
        # do it again to check the repair code
        c.mouse.motion[] = MouseButton(XY{UserUnit}(2.9, 1.4), btn.button, btn.clicktype, btn.modifiers)
        mktemp() do path, io
            redirect_stdout(io) do
                c.mouse.buttonpress[] = MouseButton(XY{UserUnit}(2.8, 1.4), 1, GtkObservables.BUTTON_PRESS, btn.modifiers)
            end
            flush(io)
            @test occursin("file3:1, f5 [inlined]", read(path, String))
        end
        fn = tempname()
        try
            ProfileView._save(fn, backtraces, lidict)
            tmp = []
            ProfileView._open(tmp, fn)
            for (gn, sn) in zip(PreOrderDFS(g), PreOrderDFS(tmp[1]))
                @test gn.data == sn.data
            end
        finally
            rm(fn)
        end

        # Also click on real stackframes
        Profile.clear()
        @profile profile_test(10)
        g = ProfileView.flamegraph()
        win, c, fdraw, (tb_open, tb_save_as) = ProfileView.viewgui(ProfileView.FlameGraphs.default_colors, g);
        ProfileView.Gtk.showall(win)
        sleep(1.0)
        sz = size(ProfileView.flamepixels(ProfileView.FlameGraphs.default_colors, g))
        mktemp() do path, io
            redirect_stdout(io) do
                for j in 1:sz[2], i in 1:sz[1]
                    c.mouse.buttonpress[] = MouseButton(XY{UserUnit}(j, i), 1, GtkObservables.BUTTON_PRESS, btn.modifiers)
                end
            end
            flush(io)
            @test occursin("MethodInstance", read(path, String))
        end

        ProfileView.closeall()
    end

    @testset "warntype_last" begin
        # Test `warntype_last`
        ProfileView.clicked[] = nothing
        @test_logs (:warn, "click on a non-inlined bar to see `code_warntype` info") warntype_last() === nothing
        _, bt = add2(Any[1,2])
        st = stacktrace(bt)
        ProfileView.clicked[] = st[1]
        io = IOBuffer()
        warntype_last(io)
        str = String(take!(io))
        @test occursin("Base.getindex(x, 1)::ANY", str)
        idx = findfirst(sf -> sf.inlined, st)
        ProfileView.clicked[] = st[idx]
        @test_logs (:warn, "click on a non-inlined bar to see `code_warntype` info") warntype_last(io) === nothing
    end
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    init_mpi()

Initialize MPI by calling `MPI.Initialized()`. The function will check if MPI is already initialized
and if yes, do nothing, thus it is safe to call it multiple times.
"""
function init_mpi()
  if MPI_INITIALIZED[]
    return nothing
  end

  if MPI.Initialized()
    @assert MPI.Query_thread() >= MPI.THREAD_FUNNELED "MPI already initialized with insufficient threading support"
  else
    # MPI.THREAD_FUNNELED: Only main thread makes MPI calls
    provided = MPI.Init_thread(MPI.THREAD_FUNNELED)
    atexit(finalize_mpi) # register atexit hook
    @assert provided >= MPI.THREAD_FUNNELED "MPI library with insufficient threading support"
  end

  # Initialize global MPI state
  MPI_RANK[] = MPI.Comm_rank(MPI.COMM_WORLD)
  MPI_SIZE[] = MPI.Comm_size(MPI.COMM_WORLD)
  MPI_IS_PARALLEL[] = MPI_SIZE[] > 1
  MPI_IS_SERIAL[] = !MPI_IS_PARALLEL[]
  MPI_IS_ROOT[] = MPI_IS_SERIAL[] || MPI_RANK[] == 0
  MPI_INITIALIZED[] = true

  return nothing
end

function finalize_mpi()
  if MPI.Initialized()
    MPI.Finalize()
  end
end

const MPI_INITIALIZED = Ref(false)
const MPI_RANK = Ref(-1)
const MPI_SIZE = Ref(-1)
const MPI_IS_PARALLEL = Ref(false)
const MPI_IS_SERIAL = Ref(true)
const MPI_IS_ROOT = Ref(true)

@inline mpi_comm() = MPI.COMM_WORLD

@inline mpi_rank() = MPI_RANK[]

@inline mpi_nranks() = MPI_SIZE[]

@inline mpi_isparallel() = MPI_IS_PARALLEL[]

# This is not type-stable but that's okay since we want to get rid of it anyway
# and it's not used in performance-critical parts. The alternative we used before,
# calling something like `eval(:(mpi_parallel() = Val(true)))` in `init_mpi()`,
# causes invalidations and slows down the first call to Trixi.
mpi_parallel()::Union{Val{true}, Val{false}} = Val(mpi_isparallel())

@inline mpi_isroot() = MPI_IS_ROOT[]

@inline mpi_root() = 0

@inline function mpi_println(args...)
  if mpi_isroot()
    println(args...)
  end
  return nothing
end
@inline function mpi_print(args...)
  if mpi_isroot()
    print(args...)
  end
  return nothing
end

end # @muladd


text:
module BenchFilter

using BenchmarkTools
using GeneratorsX
using Transducers

# some arbitrary predicate function that filters things out roughly half the time
@inline predicate(i, j, k) = xor(xor(i, j), k) < 7

@generator function genx(n)
    for i in 1:n, j in i:n, k in 1:n
        if predicate(i, j, k)
            @yield (i, j, k)
        end
    end
end

base(n) = ((i, j, k) for i in 1:n for j in i:n for k in 1:n if predicate(i, j, k))

const SUITE = BenchmarkGroup()

let s1 = SUITE[:reducer=>:collect] = BenchmarkGroup()
    s2 = s1[:impl=>:base] = BenchmarkGroup()
    s2[:src=>:base] = @benchmarkable(collect(base(10)))
    s2[:src=>:genx] = @benchmarkable(collect(genx(10)))

    s2 = s1[:impl=>:xf] = BenchmarkGroup()
    s2[:src=>:base] = @benchmarkable(collect(Map(identity), base(10)))
    s2[:src=>:genx] = @benchmarkable(collect(Map(identity), genx(10)))
end

let s1 = SUITE[:reducer=>:sum] = BenchmarkGroup()
    s2 = s1[:impl=>:base] = BenchmarkGroup()
    s2[:src=>:base] = @benchmarkable(sum(Base.splat(*), base(10)))
    s2[:src=>:genx] = @benchmarkable(sum(Base.splat(*), genx(10)))

    s2 = s1[:impl=>:xf] = BenchmarkGroup()
    s2[:src=>:base] = @benchmarkable(foldl(+, MapSplat(*), base(10)))
    s2[:src=>:genx] = @benchmarkable(foldl(+, MapSplat(*), genx(10)))

    s2 = s1[:impl=>:xf_init] = BenchmarkGroup()
    s2[:src=>:base] = @benchmarkable(foldl(+, MapSplat(*), base(10); init = 0))
    s2[:src=>:genx] = @benchmarkable(foldl(+, MapSplat(*), genx(10); init = 0))
end

end  # module
BenchFilter.SUITE


text:
#=
Authors:
Darius Russell Kish

Created: April 15, 2020
Modified: April 15, 2020

This implementation is based strongly in the Disjoint Set implementation
in JuliaCollections' DataStructures which is under active development. It
represents proper implementation of data structures for Julia and much of the
typing and integration into the language is taken from their implementation,
however the function implementations are my own.
https://github.com/JuliaCollections/DataStructures.jl/blob/master/src/disjoint_set.jl

Implementations for collections of arbitrary inputs are *not* provided here
and it is assumed this low-level data structure will be wrapped for such
functionality. This avoids issues of key-value association time complexity for
conversion between elements and their index value in analyses of this data
structure.
=#

#=
Notation T <: Integer :: T is a subtype of Integer

An example type tree:
abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end
=#
mutable struct DisjointSet{T<:Integer}
    ids::Vector{T}
    sizes::Vector{T}
    nGroups::T
end

# Construction type 1
# Notation: collect(Base.OneTo(n)) == list(range(n)) in Python
DisjointSet(n::T) where {T <: Integer} =
    DisjointSet{T}(collect(Base.OneTo(n)), ones(n), n)

# Construction type 2
DisjointSet{T}(n::Integer) where {T <: Integer} =
    DisjointSet{T}(collect(Base.OneTo(n)), ones(n), n)

# The length or "size" of a UnionFind is the length of either of its vectors
length(uf::DisjointSet) = Base.length(uf.ids)

"""
    num_groups(uf::DisjointSet)
Get a number of groups.
"""
num_groups(uf::DisjointSet) = uf.nGroups

"""
    root(uf::DisjointSet{T}, x::T)
Finds the root node of node x and performs path compression
"""
function root!(uf::DisjointSet{T}, x::T) where {T <: Integer}
    root = uf.ids[x]
    # extracting one iteration of the loop
    # allows ensuring inbounds for the rest
    @inbounds while root != uf.ids[root]
        root = uf.ids[root]
    end

    # path compression
    @inbounds while x != root
        next = uf.ids[x]
        uf.ids[x] = root
        x = next
    end

    return root
end

"""
    find(uf::DisjointSet{T}, x::T, y::T)
Determine if x and y belong to the same group
"""
function find!(uf::DisjointSet{T}, x::T, y::T) where {T <: Integer}
    return root!(uf, x) == root!(uf, y)
end

"""
    unite(uf::DisjointSet{T}, x::T, y::T)
Merges two groups together
"""
function unite!(uf::DisjointSet{T}, x::T, y::T) where {T <: Integer}
    pid = root!(uf, x)
    tid = root!(uf, y)

    if pid == tid
        return uf.nGroups
    end

    if uf.sizes[pid] < uf.sizes[tid]
        uf.ids[pid] = tid
        uf.sizes[tid] += uf.sizes[pid]
    else
        uf.ids[tid] = pid
        uf.sizes[pid] += uf.sizes[tid]
    end

    uf.nGroups -= 1
end


text:
using Test
using Distributions
using Colors
using FixedPointNumbers
using SegmentationTools
using ImageDraw
using Images
using AxisArrays
using Unitful: μm

@testset "Find peak centers" begin
    # even width peak
    xs = Float64.(1:100)
    ys = fill(0.0, 100)
    ys[20:29] .+= 2
    @test SegmentationTools.find_peak_center(xs, ys) == 24.5

    # odd width peak
    xs = Float64.(1:100)
    ys = fill(0.0, 100)
    ys[20:30] .+= 2
    @test SegmentationTools.find_peak_center(xs, ys) == 25.0

    # double peaks
    xs = Float64.(1:100)
    ys = fill(0.0, 100)
    ys[20:30] .+= 2
    ys[50:60] .+= 2
    @test_throws ErrorException("Dependent variable is not unimodal") SegmentationTools.find_peak_center(xs, ys) == 25.0

    # shoulder
    xs = Float64.(1:100)
    ys = fill(0.0, 100)
    ys[20:30] .+= 1
    ys[26:30] .+= 1
    # high cut
    @test SegmentationTools.find_peak_center(xs, ys) == 28.0
    # low cut
    @test SegmentationTools.find_peak_center(xs, ys, h=0.25) == 25.0
end

@testset "colorize" begin
    # Fixed Point Number
    data = rand(Gray{N0f16}, 1024, 1024, 10)
    # make sure the backing type is still N0f16 changed by the mapping
    @test eltype(SegmentationTools.colorize(data)) <: RGB{N0f16}

    # same for floating point numbers
    data = rand(Gray{Float64}, 1024, 1024, 10)
    @test eltype(SegmentationTools.colorize(data)) <: RGB{Float64}
end

@testset "trackpy table construction" begin
    img = zeros(Gray, 100, 100, 5)

    # create three cells
    centers = [Point(20, 20), Point(25, 17), Point(70, 70), Point(70, 70), Point(16, 80)]
    radii = [5, 5, 15, 5, 7]

    bkg = 0.1

    for i in 1:size(img, 3)
        # draw the "cells" and make them grow dimmer over time
        colors = [Gray(1.0/i) for _ in centers]
        # the 3rd circle is concentric with the 3rd and is a high background
        colors[3] = Gray(bkg)
        ImageDraw.draw!(view(img, :, :, i), [ImageDraw.CirclePointRadius(c, r) for (c,r) in zip(centers, radii)], colors)
    end

    wrapped = AxisArray(img, Axis{:y}(1μm:1μm:100μm), Axis{:x}(1μm:1μm:100μm), Axis{:time}(1:5));

    mask = AxisArray(img .> bkg, AxisArrays.axes(wrapped))
    particles = SegmentationTools.build_tp_df(wrapped, mask)

    # three timepoints and three separate particles
    @test size(particles, 1) == 3*5

    # we now assign unique ids per-frame, so we need to convert these back to
    # consistent ids across time
    ids = mod.(particles[!, :id].-1, 3) .+ 1
    # verify that the total fluorescence values in the dataframe grow dimmer as
    # we expect
    @test all(particles[ids .== 1, :tf_slice] ./ 145.0 .≈ [1.0/i for i in 1:5])

    # for the third particle verify that the total fluorescence values have the
    # correct amount of background subtracted off, i.e for a background value of
    # `bkg` the final value should be tf - area * bkg
    @test all(particles[ids .== 3, :medbkg_slice] .≈ [bkg for i in 1:5])

    # test centroids
    @test particles[1, :x] == centers[5].x
    @test particles[1, :y] == centers[5].y
end

@testset "local background calc" begin
    include("localities.jl")
end

@testset "Flatfield correction" begin
    # construct a 2D gaussian
    D = reshape(pdf.(Normal(50, 100), range(0, length=100)), :, 1) .* 100
    # 0-centered gaussian noise
    noise = Gray.(rand(Normal(0, 0.001), 100, 100))

    illum = Gray.(D * D')
    # make the signal = 2*illumination intensity
    foreground = AxisArray(illum .+ illum .+ noise)
    # background
    background = illum .+ noise

    # if the flatfield correction is working properly than the image after
    # correction should equal 2.0 everywhere since it's 2 * illum
    @test all(SegmentationTools.flatfield_correct(foreground, background, noise) .≈ 2.0)
end