text:
"""
    replace_newlines(s::AbstractString, replacement = " ")

Replace newlines with something other string. This is useful when you want to
"wrap" multiple lines into a single line. Hence the default replacement is just
a blank space.
"""
function replace_newlines(s::AbstractString, replacement = " ")
    return replace(s, r"\n" => replacement)
end

"""
    @discord(ex)

Evaluate any Discord.jl function that is expected to return a future.
The response is captured and logged at Debug level before the
underlying value is returned.
"""
macro discord(ex)
    return quote
        future = $(esc(ex))
        response = fetch(future)
        @debug "Response received" response
        response.val
    end
end

"""
    audit(source, user_id, name, discriminator, message; audit_dir, audit_file)

Create an audit record in the system for a given event.

# Arguments
- `source`: unique string representing where the audit log comes from
- `user_id`: Discord user id
- `name`: Discord user name
- `discriminator`: Discord user discriminator
- `audit_dir`: audit log directory, defaults to `data/audit`
- `audit_file`: audit log file, defaults to `\$audit_dir/\$source.log`
"""
function audit(
    source::AbstractString,
    user_id::UInt64,
    name::AbstractString,
    discriminator::AbstractString,
    message::AbstractString;
    audit_dir = joinpath("data", "audit"),
    audit_file = joinpath(audit_dir, "$source.log")
)
    mkpath(audit_dir)
    open(joinpath(audit_file); write = true, append = true) do io
        println(io, "$(now(tz"UTC"))\t$name#$discriminator\t$user_id\t$message")
    end
end

# Convenient Discord functions

"""
    delete_emoji(c::Client, channel_id::UInt64, message_id::UInt64, user_id::UInt64, emoji::Emoji)

Delete the specified user's emoji from a Discord message.
"""
function delete_emoji(c::Client, channel_id::UInt64, message_id::UInt64, user_id::UInt64, emoji::Emoji)
    message = Message(; id = message_id, channel_id = channel_id)
    user = User(; id = user_id)
    return @discord delete(c, emoji, message, user)
end

"""
    update_message(c::Client, channel_id::UInt64, message_id::UInt64, content::AbstractString)

Update a Discord message with new content.
"""
function update_message(c::Client, channel_id::UInt64, message_id::UInt64, content::AbstractString)
    message = Message(; id = message_id, channel_id = channel_id)
    @discord update(c, message; content)
end

"""
    is_real_user(u::User)

Returns true if the user is not a bot.
"""
is_real_user(u::User) = ismissing(u.bot) || u.bot === false

"""
    is_bot(u::User)

Returns true if the user is a bot.
"""
is_bot(u::User) = !is_real_user(u)

# Generic flipping functions

"""
    previous

Flip state to the previous item. If it passes the beginning of the list,
reset the state to the last item.
"""
function previous end

"""
    next

Flip state to the next item. If it passes the end of the list,
reset the state to the first item.
"""
function next end

"""
    flipper(f::Union{typeof(previous), typeof(next)}, T::Type)

Return a flip function, either `previous` or `next`, that expects
argument of type `T`.
"""
flipper(f::Union{typeof(previous), typeof(next)}, T::Type) = TypedFunction{T}(f)

"""
    extract_command(command::AbstractString, s::AbstractString)

Extract command by removing the command prefix. For example:

```
julia> HoJBot.extract_command("j", ",j doc sin")
"doc sin"
```
"""
function extract_command(command::AbstractString, s::AbstractString)
    prefix = COMMAND_PREFIX * command
    return strip(replace(s, Regex("^" * prefix * " *") => ""))
end


text:
# Autogenerated wrapper script for NativeFileDialog_jll for i686-w64-mingw32
export libnfd

JLLWrappers.@generate_wrapper_header("NativeFileDialog")
JLLWrappers.@declare_library_product(libnfd, "libnfd.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnfd,
        "bin\\libnfd.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using JSON

## build_groups
f = jldopen("data/small.ChromNet.jld")
data = read(f, "data")
ids = read(f, "ids")
metadata = read(f, "metadata")
close(f)
C = cov([data' data'])
groups = build_groups(C, [ids; ids])
for i in 1:8:length(groups[end][2])
    @test groups[end][2][i] == groups[end][2][i+4] # make sure duplicates are grouped together
end

## build_groupgm
C = cov([data[1,:]' data[1,:]' data[2,:]'.+0.01*data[1,:]' data[3,:]'] .+ 0.0001*randn(size(data)[2], 4))
groups = build_groups(C, ids[1:4])
I = inv(C)
G,idsG = build_groupgm(I, ids[1:4], groups)
@test abs(I[1,3] + I[2,3] - G[5,3]) < 1e-7

## parse_config
metadataTmp = open(f->parse_config(f, "data"), "data/single_bed.config")
@test metadataTmp["EX001"]["name"] == "K562_tss"
@test metadataTmp["EX001"]["lab"] == "Sample lab"

## build_network
C = cor(data, vardim=2)
groups = build_groups(C, ids)
G,idsG = build_groupgm(inv(C), ids, groups)
network = build_network(G, groups, metadata; threshold=0.03, groupScoreThreshold=0.7, quiet=true)


text:
using PyPlot
using JLD

path = "/Users/milan/cams"
allvars = ["no2","go3","so2","aermr04","aermr05","aermr06","co","no"]
allvarsLN = [L"NO$_2$",L"O$_3$",L"SO$_2$","AERMR04","AERMR05","AERMR06","CO","NO"]

D = load("cams/error/$vari.jld")
EV = D["EV"]
lz4hcc = D["lz4"]

EVs = Array{Float32,3}(undef,length(allvars),size(EV)[1],size(EV)[2])
LZs = Array{Float32,1}(undef,length(allvars))

for (i,vari) in enumerate(allvars)
    D = load("cams/error/$vari.jld")
    EVs[i,:,:] = D["EV"]
    LZs[i] = D["lz4"]
end

## theoretical
theosize = [32/i for i in 32:-1:1]
x0 = 1f0
theoerror = [log10(nextfloat(x0,2^i)) for i in 0:31]

## plotting
pygui(true)

cm1 = ColorMap("viridis")
cm2 = ColorMap("inferno")

function RGB(cv::Tuple)
    r = UInt8(round(cv[1]*255))
    g = UInt8(round(cv[2]*255))
    b = UInt8(round(cv[3]*255))
    return "#"*repr(r)[3:4]*repr(g)[3:4]*repr(b)[3:4]
end

colours = ["0.6","0.25",RGB(cm2(0.8)),RGB(cm2(0.6)),RGB(cm2(0.4)),RGB(cm2(0.9)),
                    RGB(cm1(0.6)),RGB(cm1(0.4)),RGB(cm1(0.2)),RGB(cm1(0.8))]

names = ["LinQuant24","LogQuant16","RoundNearest16+Blosc","RoundNearest16+LZ4HC",
            "RoundNearest16+Zstd","RoundNearest16+Deflate","RoundNearest24+Blosc",
            "RoundNearest24+LZ4HC","RoundNearest24+Zstd","RoundNearest24+Deflate"]

# markers = ["s","s","d","d","d","d","o","o","o","o"]
markers = ["p","^","o","+","x","4","h","v"]

fig,ax = subplots(1,1,figsize=(8,6))
ax.set_xscale("log")
ax2 = ax.twinx()

ax2.plot(theoerror[1],1,"w",markersize=15,marker="*",lw=0,markeredgewidth=1.4,
            markeredgecolor="k",alpha=1,label="Float32",zorder=10)

for j in 1:length(allvars)
    ax.scatter(theoerror[1],LZs[j],120,color="w",marker=markers[j],
            edgecolor="C3",alpha=0.7)
end

for j in [1,8,2,7,3,4,5,6]
    if j in [4,5,6]
        ax.scatter([-1,-1,-1],[-1,-1,-1],120,color="k",edgecolor="k",marker=markers[j],
            alpha=0.55,lw=1.5,label=allvarsLN[j])
        ax.scatter(theoerror[1],LZs[j],120,color="C3",marker=markers[j],alpha=0.7)
    else
        ax.scatter([-1,-1,-1],[-1,-1,-1],120,color="w",edgecolor="k",marker=markers[j],
            alpha=0.55,lw=1.5,label=allvarsLN[j])
        ax.scatter(theoerror[1],LZs[j],120,color="w",marker=markers[j],
            edgecolor="C3",alpha=0.7)
    end
    ax.scatter(EVs[j,:,5],EVs[j,:,1],120,colours,marker=markers[j],
        alpha=0.55)
end

for i in [1,2,6,3,4,5,10,7,8,9]
    ax2.scatter(-1,-1,120,colours[i],marker="s",alpha=0.7,label=names[i])
end

ax2.plot(theoerror,theosize,"k",label="Theoretical rounding",zorder=-1)
ax.set_xlim(theoerror[1]/2,20)
ax.set_ylim(0.7,10)
ax2.set_ylim(0.7,10)
ax2.set_yticks([])
ax.set_yticks([1,2,4,6,8,10])

# ax.set_ylim(0,110)
ax.set_xlabel("max decimal error")
ax.set_ylabel("compression factor")
ax.legend(loc=2,ncol=1,fontsize=9,scatterpoints=1,title="Variable")
ax2.legend(loc=1,ncol=1,fontsize=9,scatterpoints=1,title="Compression method")
ax.set_title("CAMS data compression with various lossy/lossless methods",loc="left")

# ax.text(EV[1,3],EV[1,1],"median  ",rotation=90,ha="center",va="top",fontsize=9)
# ax.text(EV[1,4],EV[1,1],"90%  ",rotation=90,ha="center",va="top",fontsize=9)
# ax.text(EV[1,5],EV[1,1],"max  ",rotation=90,ha="center",va="top",fontsize=9)

tight_layout()


text:
using VisualRegressionTests
using ExamplePlots

import DataFrames, RDatasets

# don't let pyplot use a gui... it'll crash
# note: Agg will set gui -> :none in PyPlot
ENV["MPLBACKEND"] = "Agg"
try
  @eval import PyPlot
  info("Matplotlib version: $(PyPlot.matplotlib[:__version__])")
end

using Plots, FactCheck

default(size=(500,300))

# TODO: use julia's Condition type and the wait() and notify() functions to initialize a Window, then wait() on a condition that
#       is referenced in a button press callback (the button clicked callback will call notify() on that condition)

function image_comparison_tests(pkg::Symbol, idx::Int; debug = false, popup = isinteractive(), sigma = [1,1], eps = 1e-2)

  # first
  Plots._debugMode.on = debug
  example = ExamplePlots._examples[idx]
  info("Testing plot: $pkg:$idx:$(example.header)")
  backend(pkg)
  backend()

  # ensure consistent results
  srand(1234)

  # reference image directory setup
  refdir = joinpath(Pkg.dir("ExamplePlots"), "test", "refimg", string(pkg))

  # test function
  func = (fn, idx) -> begin
    map(eval, example.exprs)
    png(fn)
  end

  try
    run(`mkdir -p $refdir`)
  catch err
    display(err)
  end
  reffn = joinpath(refdir, "ref$idx.png")

  # the test
  vtest = VisualTest(func, reffn, idx)
  test_images(vtest, popup=popup, sigma=sigma, eps=eps)
end

function image_comparison_facts(pkg::Symbol;
                                skip = [],      # skip these examples (int index)
                                only = nothing, # limit to these examples (int index)
                                debug = false,  # print debug information?
                                sigma = [1,1],  # number of pixels to "blur"
                                eps = 1e-2)     # acceptable error (percent)
  for i in 1:length(ExamplePlots._examples)
    i in skip && continue
    if only == nothing || i in only
      @fact image_comparison_tests(pkg, i, debug=debug, sigma=sigma, eps=eps) |> success --> true
    end
  end
end


text:
function Distributions.minimum(d::Distributions.Gamma{T}) where {T<:Real}
    return 0.0
end
function Distributions.maximum(d::Distributions.Gamma{T}) where {T<:Real}
    return Inf
end


text:
# run all 4 processor tests in a single session
args_orig = copy(ARGS)

cd("./advection")
include(joinpath(pwd(), "runtests_parallel4.jl"))

cd("../euler")
resize!(ARGS, length(args_orig))  # reset to default
copy!(ARGS, args_orig)
include(joinpath(pwd(), "runtests_parallel4.jl"))

cd("../")


text:
"""
    fillnodata_upstream(g, toposort, data, nodata)

Fill `nodata` upstream cells with the value from the first downstream valid cell, based on
directed acyclic graph `g`, topological order `toposort`, nodata value `nodata` and `data`
containing missing values. Returns filled data `data_out`.
"""
function fillnodata_upstream(g, toposort, data, nodata)
    data_out = copy(data)
    for v in reverse(toposort)  # down- to upstream
        idx_ds = outneighbors(g, v)
        if !isempty(idx_ds)
            if data_out[v] == nodata && data_out[only(idx_ds)] != nodata
                data_out[v] = data_out[only(idx_ds)]
            end
        end
    end
    return data_out
end

"""
    stream_order(g, toposort)

Return the Strahler streamorder based on directed acyclic graph `g` and topological order
`toposort`.
"""
function stream_order(g, toposort)
    n = length(toposort)
    strord = fill(1, n)
    for v in toposort
        inds_up = inneighbors(g, v)
        if !isempty(inds_up)
            sto_up = strord[inds_up]
            if length(findall(x -> x == maximum(sto_up), sto_up)) > 1
                strord[v] = (maximum(sto_up) + 1)
            else
                strord[v] = maximum(sto_up)
            end
        end
    end
    return strord
end

"""
    subbasins(g, streamorder, toposort, min_sto)

Return subbasins with a unique id starting at 1, based on a minimum streamorder `min_sto`,
directed acyclic graph `g` and topological order `toposort`.
"""
function subbasins(g, streamorder, toposort, min_sto)
    n = length(toposort)
    subbas = fill(0, n)

    i = 1
    for v in toposort
        if streamorder[v] < min_sto
            continue
        end
        ds_nodes = outneighbors(g, v)
        if !isempty(ds_nodes)
            if streamorder[v] != streamorder[only(ds_nodes)]
                subbas[v] = i
                i += 1
            end
        else
            # also set pits (without a downstream node)
            subbas[v] = i
            i += 1
        end
    end
    return subbas
end

"""
    subbasins_order(g, outlet, max_dist)

Group subbasins (down- to upstream), starting at the `outlet` of the basin, per distance (1
to maximum distance `max_dist`) from the `outlet`, based on the directed acyclic graph `g`
of subbasins. Returns grouped subbasins `order`, ordered from `max_dist` (including
subbasins without upstream neighbor and distance < `max_dist`) to distance 0 (`outlet`).
"""
function subbasins_order(g, outlet, max_dist)
    order = Vector{Vector{Int}}(undef, max_dist + 1)
    order[1] = [outlet]
    for i = 1:max_dist
        v = Vector{Int}()
        for n in order[i]
            ups_nodes = inneighbors(g, n)
            if !isempty(ups_nodes)
                append!(v, ups_nodes)
            end
        end
        order[i+1] = v
    end

    # move subbasins without upstream neighbor (headwater) to index [max_dist+1]
    for i = 1:max_dist
        for s in order[i]
            if isempty(inneighbors(g, s))
                append!(order[max_dist+1], s)
                filter!(e -> e ≠ s, order[i])
            end
        end
    end
    return reverse(order)
end

"""
    graph_from_nodes(graph, subbas, subbas_fill)

Extract directed acyclic graph `g` representing the flow network at subbasin level from
`subbas` containing nodes with a unique id representing subbasin outlets, `subbas_fill` with
subbasin ids for the complete domain, and directed acyclic graph `graph` representing the
flow network for each subbasin cell.
"""
function graph_from_nodes(graph, subbas, subbas_fill)
    n = maximum(subbas)
    g = DiGraph(n)
    for i = 1:n
        idx = findall(x -> x == i, subbas)
        ds_idx = outneighbors(graph, only(idx))
        to_node = subbas_fill[ds_idx]
        if !isempty(to_node)
            add_edge!(g, i, only(to_node))
        end
    end
    return g
end

"""
    function kinwave_set_subdomains(config, graph, toposort, index_pit)

Setup subdomains for parallel execution (threading) of the kinematic wave calculation.
Subdomains are subbasins based on a minumum streamorder `min_sto` (see also `subbasins(g,
streamorder, toposort, min_sto)`). Subbasins are extracted for each basin outlet in Vector
`index_pit`.

# Arguments
- `config` Struct that contains parsed TOML configuration
- `graph` directed acyclic graph of the kinematic wave domain
- `toposort` topological order of `graph`
- `index_pit` Vector with basin outlets (pits)

# Output
- `subbas_order` grouped subbasin ids (`Vector{Vector{Int}}`) ordered upstream (first index)
  to downstream (last index)
- `indices_subbas` list of indices per subbasin id stored as `Vector{Vector{Int}}`
- `topo_subbas` topological order per subbasin id stored as `Vector{Vector{Int}}`
"""
function kinwave_set_subdomains(config, graph, toposort, index_pit)

    if nthreads() > 1
        min_sto = get(config.model, "min_streamorder", 4)

        # extract basins (per outlet/pit), assign unique basin id
        n_pits = length(index_pit)
        basin = fill(0, length(toposort))
        basin[index_pit] = [1:n_pits;]
        basin_fill = fillnodata_upstream(graph, toposort, basin, 0)

        # pre-allocate the Vector with indices matching the topological order of the
        # complete domain (upstream neighbors are stored at these indices)
        index_toposort = fill(0, length(toposort))
        for (i, j) in enumerate(toposort)
            index_toposort[j] = i
        end

        order_subbas = Vector{Vector{Int}}()
        indices_subbas = Vector{Vector{Int}}()
        topo_subbas = Vector{Vector{Int}}()
        index = Vector{Int}()
        total_subbas = 0
        for i = 1:n_pits
            # extract subbasins per basin, make a graph at the subbasin level, calculate the
            # maximum distance of this graph, and group and order the subbasin ids from
            # upstream to downstream
            basin = findall(x -> x == i, basin_fill)
            g, _ = induced_subgraph(graph, basin)
            toposort_b = topological_sort_by_dfs(g)
            streamorder = stream_order(g, toposort_b)
            subbas = subbasins(g, streamorder, toposort_b, min_sto)
            subbas_fill = fillnodata_upstream(g, toposort_b, subbas, 0)
            n_subbas = max(length(subbas[subbas.>0]), 1)
            if n_subbas > 1
                graph_subbas = graph_from_nodes(g, subbas, subbas_fill)
                toposort_subbas = topological_sort_by_dfs(graph_subbas)
                dist = Graphs.Experimental.Traversals.distances(
                    Graph(graph_subbas),
                    toposort_subbas[end],
                )
                max_dist = maximum([dist; 1])
                v_subbas = subbasins_order(graph_subbas, toposort_subbas[end], max_dist)
            else
                v_subbas = [[1]]
            end
            # subbasins need a unique id (in case of multiple basins/outlets in the
            # kinematic wave domain)
            for n = 1:length(v_subbas)
                v_subbas[n] .= v_subbas[n] .+ total_subbas
            end
            total_subbas += n_subbas
            append!(order_subbas, v_subbas)
            append!(index, 1:length(v_subbas))
            # in case of multiple subbasins calculate topological order per subbasin
            # (subgraph of the corresponding basin graph g), and the indices that match the
            # subbasin topological order
            if n_subbas > 1
                for s = 1:n_subbas
                    subbas_s = findall(x -> x == s, subbas_fill)
                    sg, _ = induced_subgraph(g, subbas_s)
                    toposort_sg = topological_sort_by_dfs(sg)
                    push!(topo_subbas, basin[subbas_s[toposort_sg]])
                    push!(indices_subbas, index_toposort[basin[subbas_s[toposort_sg]]])
                end
            else
                push!(topo_subbas, basin[toposort_b])
                push!(indices_subbas, index_toposort[basin[toposort_b]])
            end
        end
        # reduce the order of subbasin ids by merging groups of subbasins that have the same
        # index (multiple basins/outlets in the kinematic wave domain)
        subbas_order = Vector{Vector{Int}}(undef, maximum(index))
        for m = 1:maximum(index)
            subbas_order[m] = reduce(vcat, order_subbas[index.==m])
        end
    else
        subbas_order = [[1]]
        indices_subbas = [[1:length(toposort);]]
        topo_subbas = [toposort]
    end

    return subbas_order, indices_subbas, topo_subbas
end


text:
using NodesAndModes
using Test
using LinearAlgebra

@testset "1D tests" begin
    tol = 1e2*eps()

    N = 0
    r,w = gauss_lobatto_quad(0,0,N)
    @test sum(w) ≈ 2

    for N = 1:2
        r,w = gauss_quad(0,0,N)
        @test sum(w) ≈ 2
        @test abs(sum(w.*r)) < tol

        V = vandermonde(Line(),N,r)
        @test V'*diagm(w)*V ≈ I

        Vr = grad_vandermonde(Line(),N,r)
        Dr = Vr/V
        @test norm(sum(Dr,dims=2)) < tol
        @test Dr*r ≈ ones(N+1)

        r,w = gauss_lobatto_quad(0,0,N)
        @test sum(w) ≈ 2
        @test abs(sum(w.*r)) < tol

        # check if Lobatto is exact for 2N-2 polynoms
        V = vandermonde(Line(),N-1,r)
        @test V'*diagm(w)*V ≈ I

        Dr = grad_vandermonde(Line(),N,r)/vandermonde(Line(),N,r)
        @test norm(sum(Dr,dims=2)) < tol
        @test Dr*r ≈ ones(N+1)
    end
end

area(elem::Quad) = 4.0
area(elem::Tri) = 2.0

@testset "2D $elem basis tests" for elem in (Tri(), Quad())

    tol = 1e2*eps()

    N = 3
    rq, sq, wq = quad_nodes(elem, N)
    @test sum(wq) ≈ area(elem)

    Vq = vandermonde(elem, N, rq, sq)
    @test Vq' * diagm(wq) * Vq ≈ I

    r, s = nodes(elem, N)
    V = vandermonde(elem, N, r, s)
    Dr, Ds = (A->A / V).(grad_vandermonde(elem, N, r, s))
    @test norm(sum(Dr, dims=2)) + norm(sum(Ds, dims=2)) < tol
    @test norm(Dr * s) + norm(Ds * r) < tol
    @test Dr * r ≈ ones(length(r))
    @test Ds * s ≈ ones(length(s))

    r, s = equi_nodes(elem, N)
    V = vandermonde(elem, N, r, s)
    Dr,Ds = (A->A / V).(grad_vandermonde(elem, N, r, s))
    @test norm(sum(Dr, dims=2)) + norm(sum(Ds, dims=2)) < tol
    @test norm(Dr * s) + norm(Ds * r) < tol
    @test Dr*r ≈ ones(length(r))
    @test Ds*s ≈ ones(length(s))

    # check for Kronecker structure
    if elem == Quad()
        r, s = nodes(elem, N)
        V = vandermonde(elem, N, r, s)
        V1D = vandermonde(Line(), N, nodes(Line(), N))
        @test V ≈ kron(V1D, V1D)
    end
end

volume(elem::Hex) = 8
volume(elem::Wedge) = 4
volume(elem::Pyr) = 8/3
volume(elem::Tet) = 4/3

@testset "3D $elem basis tests" for elem in (Hex(), Wedge(), Pyr(), Tet())

    tol = 5e2 * eps()

    N = 3
    rq, sq, tq, wq = quad_nodes(elem,N)
    @test sum(wq) ≈ volume(elem)
    if elem==Hex()
        @test abs(sum(rq .* wq)) + abs(sum(sq .* wq)) + abs(sum(tq .* wq)) < tol
    end

    Vq = vandermonde(elem, N, rq, sq, tq)
    @test Vq' * diagm(wq) * Vq ≈ I

    if elem != Pyr()
        r, s, t = nodes(elem, N)
        V = vandermonde(elem, N, r, s, t)
        Dr, Ds, Dt = (A->A / V).(grad_vandermonde(elem, N, r, s, t))
        @test norm(sum(Dr, dims=2)) + norm(sum(Ds, dims=2)) + norm(sum(Dt, dims=2)) < tol
        @test norm(Dr * s) + norm(Dr * t) < tol
        @test norm(Ds * r) + norm(Ds * t) < tol
        @test norm(Dt * r) + norm(Dt * s) < tol
        @test Dr * r ≈ ones(length(r))
        @test Ds * s ≈ ones(length(s))
        @test Dt * t ≈ ones(length(t))

        r, s, t = equi_nodes(elem, N)
        V = vandermonde(elem, N, r, s, t)
        Dr, Ds, Dt = (A->A / V).(grad_vandermonde(elem, N, r, s, t))
        @test norm(sum(Dr, dims=2)) + norm(sum(Ds, dims=2)) + norm(sum(Dt, dims=2)) < tol
        @test norm(Dr * s) + norm(Dr * t) + norm(Ds * r) + norm(Ds * t) + norm(Dt * r) + norm(Dt * s) < tol
        @test Dr * r ≈ one.(r)
        @test Ds * s ≈ one.(s)
        @test Dt * t ≈ one.(t)
    elseif elem == Pyr()
        r, s, t = nodes(elem, N)
        rq, sq, tq, wq = quad_nodes(elem, N)
        Vq, Vr, Vs, Vt = (A -> A / vandermonde(elem, N, r, s, t)).(basis(elem, N, rq, sq, tq))
        M = Vq' * diagm(wq) * Vq

        # construct weak differentiation matrices instead of nodal differentiation matrices since
        # the pyramid basis is singular at the node at the tip of the pyramid.
        Dr, Ds, Dt = (A -> M \ (Vq' * diagm(wq) * A)).((Vr, Vs, Vt))
        @test norm(sum(Dr, dims=2)) < tol
        @test norm(sum(Ds, dims=2)) < tol
        @test norm(sum(Dt, dims=2)) < tol
        @test norm(Dr * s) + norm(Dr * t) < tol
        @test norm(Ds * r) + norm(Ds * t) < tol
        @test norm(Dt * r) + norm(Dt * s) < tol
        @test Dr * r ≈ one.(r)
        @test Ds * s ≈ one.(s)
        @test Dt * t ≈ one.(t)
    end
end

@testset "Test for Kronecker structure in the Hex basis matrix" begin
    tol = 5e2*eps()
    N = 3
    VDM_1D = vandermonde(Line(), N, nodes(Line(), N))
    VDM_quad_1D = vandermonde(Line(), N, quad_nodes(Line(), N)[1])
    @test abs(norm(vandermonde(Hex(), N, nodes(Hex(), N)...) - kron(VDM_1D, VDM_1D, VDM_1D))) < tol
    @test abs(norm(vandermonde(Hex(), N, quad_nodes(Hex(), N)[1:3]...) - kron(VDM_quad_1D, VDM_quad_1D, VDM_quad_1D))) < tol
end

# test high order quadrature
@testset "Simplicial Stroud quadrature" begin
    tol = 1e3*eps()

    N = 14
    rq, sq, wq = quad_nodes(Tri(), N)
    @test sum(wq) ≈ 2
    @test sum(rq .* wq) ≈ -2/3
    @test sum(sq .* wq) ≈ -2/3

    N = 8
    rq, sq, tq, wq = quad_nodes(Tet(), N)
    @test sum(wq) ≈ 4/3
    @test sum(rq .* wq) ≈ -2/3
    @test sum(sq .* wq) ≈ -2/3
    @test sum(tq .* wq) ≈ -2/3
end

if VERSION >= v"1.6" # apparently inference got better after 1.5?
    @testset "Inferrability for $elementType" for elementType in (Line(), Tri(), Quad(), Tet(), Pyr(), Wedge(), Hex())
        N = 1
        if elementType == Line()
            @test_nowarn (@inferred gauss_quad(0, 0, N)) == ([-0.5773502691896257, 0.5773502691896257], [0.9999999999999998, 0.9999999999999998])
            @test_nowarn (@inferred gauss_lobatto_quad(0, 0, N)) == ([-1.0, 1.0], [1.0, 1.0])
        end

        @test_nowarn @inferred nodes(elementType, N)
        @test_nowarn @inferred quad_nodes(elementType, N)
        @test_nowarn @inferred equi_nodes(elementType, N)
        if elementType == Line()
            @test_nowarn @inferred basis(elementType, N, nodes(elementType,N))
        else
            @test_nowarn @inferred basis(elementType, N, nodes(elementType,N)...)
        end
    end
end


text:
module WindowsTimeZoneIDs

using LazyArtifacts
using ...TimeZones: DEPS_DIR, @artifact_str
using Future: copy!

const UNICODE_CLDR_VERSION = "release-39"

# A mapping of Windows timezone names to Olson timezone names.
# Details on the contents of this file can be found at:
# http://cldr.unicode.org/development/development-process/design-proposals/extended-windows-olson-zid-mapping
const WINDOWS_ZONE_URL = "https://raw.githubusercontent.com/unicode-org/cldr/$UNICODE_CLDR_VERSION/common/supplemental/windowsZones.xml"
const WINDOWS_ZONE_FILE = joinpath("cldr-$UNICODE_CLDR_VERSION", "common", "supplemental", "windowsZones.xml")

const WINDOWS_XML_DIR = joinpath(DEPS_DIR, "local")
const WINDOWS_XML_FILE = joinpath(WINDOWS_XML_DIR, "windowsZones.xml")

const WINDOWS_TRANSLATION = Dict{String, String}()

function __init__()
    isdir(WINDOWS_XML_DIR) || mkdir(WINDOWS_XML_DIR)

    if isfile(WINDOWS_XML_FILE)
        copy!(WINDOWS_TRANSLATION, compile(WINDOWS_XML_FILE))
    end
end

function compile(xml_file::AbstractString)

    translation = Dict{String,String}()

    # Get the timezone conversions from the file
    #
    # Note: Since the XML file is simplistic enough that we can parse what we need via a
    # regex we can avoid having an XML package dependency. Additionally, since this XML file
    # is included as part of the this package we can correct any parsing issues before a
    # TimeZones.jl release occurs.
    for line in readlines(xml_file)
        # Territory "001" is the global default
        occursin("territory=\"001\"", line) || continue
        win_name = match(r"other=\"(.*?)\"", line)[1]
        posix_name = match(r"type=\"(.*?)\"", line)[1]
        translation[win_name] = posix_name
    end

    return translation
end

function build(xml_file::AbstractString=WINDOWS_XML_FILE; force::Bool=false)
    if !isfile(xml_file) || force
        @info "Downloading Windows to POSIX timezone ID XML version: $UNICODE_CLDR_VERSION"
        download(WINDOWS_ZONE_URL, xml_file)
    end

    @info "Compiling Windows time zone name translation"
    copy!(WINDOWS_TRANSLATION, compile(xml_file))
end

end


text:
using DiffEqSensitivity
using OrdinaryDiffEq, Calculus, Test
using Zygote

function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -p[1]
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1]
end

function affect!(integrator)
  @show integrator.t
  println("bounced.")
  integrator.u[2] = -integrator.p[2]*integrator.u[2]
end

cb = ContinuousCallback(condition, affect!)
p = [9.8, 0.8]
prob = ODEProblem(f,eltype(p).([1.0,0.0]),eltype(p).((0.0,1.0)),copy(p))

function test_f(p)
  _prob = remake(prob, p=p)
  solve(_prob,Tsit5(),abstol=1e-14,reltol=1e-14,callback=cb,save_everystep=false)[end]
end
findiff = Calculus.finite_difference_jacobian(test_f,p)
findiff

using ForwardDiff
ad = ForwardDiff.jacobian(test_f,p)
ad

@test ad ≈ findiff

function test_f2(p, sensealg=ForwardDiffSensitivity(), controller=nothing, alg=Tsit5())
  _prob = remake(prob, p=p)
  u = solve(_prob,alg,sensealg=sensealg,controller=controller,
    abstol=1e-14,reltol=1e-14,callback=cb,save_everystep=false)
  u[end][end]
end

@test test_f2(p) == test_f(p)[end]

g1 = Zygote.gradient(θ->test_f2(θ,ForwardDiffSensitivity()), p)
g2 = Zygote.gradient(θ->test_f2(θ,ReverseDiffAdjoint()), p)
g3 = Zygote.gradient(θ->test_f2(θ,ReverseDiffAdjoint(), IController()), p)
g4 = Zygote.gradient(θ->test_f2(θ,ReverseDiffAdjoint(), PIController(7//50, 2//25)), p)
@test_broken g5 = Zygote.gradient(θ->test_f2(θ,ReverseDiffAdjoint(), PIDController(1/18. , 1/9., 1/18.)), p)
g6 = Zygote.gradient(θ->test_f2(θ,ForwardDiffSensitivity(),
  OrdinaryDiffEq.PredictiveController(), TRBDF2()), p)
@test_broken g7 = Zygote.gradient(θ->test_f2(θ,ReverseDiffAdjoint(),
    OrdinaryDiffEq.PredictiveController(), TRBDF2()), p)

@test g1[1] ≈ findiff[2,1:2]
@test g2[1] ≈ findiff[2,1:2]
@test g3[1] ≈ findiff[2,1:2]
@test g4[1] ≈ findiff[2,1:2]
@test_broken g5[1] ≈ findiff[2,1:2]
@test g6[1] ≈ findiff[2,1:2]
@test_broken g7[1] ≈ findiff[2,1:2]


text:
function visualize!(vis, p::Rocket, q;
	Δt=0.1, mesh=false, T_off=length(q))
	default_background!(vis)

	if mesh
		obj_rocket = joinpath(pwd(), "models/starship/Starship.obj")
		mtl_rocket = joinpath(pwd(), "models/starship/Starship.mtl")
		ctm = ModifiedMeshFileObject(obj_rocket, mtl_rocket, scale=1.0)
		setobject!(vis["rocket"]["starship"], ctm)

		settransform!(vis["rocket"]["starship"],
			compose(Translation(0.0, 0.0, -p.length),
				LinearMap(0.25 * RotY(0.0) * RotZ(0.5 * π) * RotX(0.5 * π))))

        body = Cylinder(Point3f0(0.0, 0.0, -1.25),
          Point3f0(0.0, 0.0, 0.5),
          convert(Float32, 0.125))

        setobject!(vis["rocket"]["body"], body,
          MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, 1.0)))
	else
		body = Cylinder(Point3f0(0.0, 0.0, -1.0 * p.length),
			Point3f0(0.0, 0.0, 1.0 * p.length),
			convert(Float32, 0.15))

		setobject!(vis["rocket"], body,
			MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
			anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))
	end

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for t = 1:length(q)
	    MeshCat.atframe(anim, t) do
            if t >= T_off
				setvisible!(vis["rocket"]["body"], false)
			else
				setvisible!(vis["rocket"]["body"], true)
			end
	        settransform!(vis["rocket"],
	              compose(Translation(q[t][1:3]),
	                    LinearMap(MRP(q[t][4:6]...) * RotX(0.0))))
	    end
	end

	MeshCat.setanimation!(vis, anim)
end


text:
using FactCheck
using Distributions

srand(12347)

facts("Checking HMM distribution") do
    K = 4
    T = 50
    pars = Dirichlet(K, 1)
    A = Array{Float64}(K, K)
    π0 = rand(pars)
    for k in 1:K
        A[:, k] = rand(pars)
    end
    ψ = rand(K, T)

    x = HMM(ψ, π0, A)

    context("Constructors") do
        @fact isa(x, HMM) --> true
        @fact isa(x, DiscreteMatrixDistribution) --> true
        @fact isa(x, Distribution) --> true
        @fact x.ψ --> ψ
        @fact x.A --> A
        @fact x.π0 --> π0
        @fact size(x.ξ) --> size(ψ)
        @fact size(x.Ξ) --> (K, K, T - 1)
    end

    context("Bad constructor inputs throw errors") do
        @fact_throws ErrorException HMM(ψ[2:end, :], π0, A)
        @fact_throws ErrorException HMM(ψ, π0[2:end], A)
        @fact_throws ErrorException HMM(ψ, π0, A[2:end, :])
        @fact_throws ErrorException HMM(ψ, π0, A[:, 2:end])
    end

    context("Conversions") do
        ψ32 = convert(Matrix{Float32}, ψ)
        π32 = convert(Vector{Float32}, π0)
        A32 = convert(Matrix{Float32}, A)
        x32 = HMM(ψ32, π32, A32)
        @fact isa(convert(HMM{Float64}, ψ32, π32, A32), HMM{Float64}) --> true
        @fact isa(convert(HMM{Float64}, x32), HMM{Float64}) --> true
    end

    context("Check basic interface") do
        @fact nstates(x) --> K
        @fact size(x) --> (K, T)
        @fact length(x) --> K * T
    end

    context("Check mean") do
        ξ = mean(x)
        @fact ξ --> x.ξ
    end

    context("Check two-slice") do
        Ξ = cov(x)
        @fact Ξ --> x.Ξ
    end

    context("Check sampling") do
        z = rand(x)
        @fact size(z) --> (K, T)
        @fact length(find(z)) --> T
        @fact sum(z .!= 0, 1) --> ones(1, T)
    end

    context("Check logpdf") do
        # draw a chain of states
        z = zeros(K, T)
        init_state = rand(Categorical(π0))
        z[init_state, 1] = 1
        for t in 2:T
            newstate = rand(Categorical(A * z[:, t - 1]))
            z[newstate, t] = 1
        end
        lpdf = logpdf(x, z)

        @fact lpdf --> less_than_or_equal(0)
    end

    context("Check entropy") do
        @fact entropy(x) --> greater_than_or_equal(0)
    end

    context("Check natural parameters") do
        nats = naturals(x)
        pars = naturals_to_params(nats, typeof(x))

        @fact length(nats) --> 3
        @fact size(nats[1]) --> size(x)
        @fact length(pars) --> 3
        @fact size(pars[1]) --> size(x)
        @fact nats[1] --> roughly(log(ψ))
        @fact nats[2] --> roughly(log(π0))
        @fact nats[3] --> roughly(log(A))
        @fact pars[1] --> roughly(ψ)
        @fact pars[2] --> roughly(π0)
        @fact pars[3] --> roughly(A)
    end

    context("Check params") do
        @fact params(x) --> (x.ψ, x.π0, x.A)
    end
end

facts("Checking MarkovChain distribution") do
    K = 4
    T = 50
    pars = Dirichlet(K, 1)
    A = Array{Float64}(K, K)
    π0 = rand(pars)
    for k in 1:K
        A[:, k] = rand(pars)
    end

    x = MarkovChain(π0, A, T)

    context("Constructors") do
        @fact isa(x, MarkovChain) --> true
        @fact isa(x, DiscreteMatrixDistribution) --> true
        @fact isa(x, Distribution) --> true
        @fact x.A --> A
        @fact x.π0 --> π0
    end

    context("Bad constructor inputs throw errors") do
        @fact_throws ErrorException MarkovChain(π0, A, 5.5)
        @fact_throws ErrorException MarkovChain(π0[2:end], A, T)
        @fact_throws ErrorException MarkovChain(π0, A[2:end, :], T)
        @fact_throws ErrorException MarkovChain(π0, A[:, 2:end], T)
    end

    context("Conversions") do
        π32 = convert(Vector{Float32}, π0)
        A32 = convert(Matrix{Float32}, A)
        x32 = MarkovChain(π32, A32, T)
        @fact isa(convert(MarkovChain{Float64}, π32, A32, T), MarkovChain{Float64}) --> true
        @fact isa(convert(MarkovChain{Float64}, x32), MarkovChain{Float64}) --> true
    end

    context("Check basic interface") do
        @fact nstates(x) --> K
        @fact size(x) --> (K, T)
        @fact length(x) --> K * T
    end

    context("Check mean") do
        ξ = mean(x)
        @fact ξ[:, 1] --> π0
        @fact ξ[:, 2] --> A * π0
    end

    context("Check sampling") do
        z = rand(x)
        @fact size(z) --> (K, T)
        @fact length(find(z)) --> T
        @fact sum(z .!= 0, 1) --> ones(1, T)
    end

    context("Check logpdf") do
        # draw a chain of states
        z = zeros(K, T)
        init_state = rand(Categorical(π0))
        z[init_state, 1] = 1
        for t in 2:T
            newstate = rand(Categorical(A * z[:, t - 1]))
            z[newstate, t] = 1
        end
        lpdf = logpdf(x, z)

        @fact lpdf --> less_than_or_equal(0)
    end

    context("Check entropy") do
        @fact entropy(x) --> greater_than_or_equal(0)
    end

    context("Check natural parameters") do
        nats = naturals(x)
        pars = naturals_to_params(nats, typeof(x))

        @fact length(nats) --> 2
        @fact length(pars) --> 2
        @fact nats[1] --> roughly(log(π0))
        @fact nats[2] --> roughly(log(A))
        @fact pars[1] --> roughly(π0)
        @fact pars[2] --> roughly(A)
    end

    context("Check params") do
        @fact params(x) --> (x.π0, x.A)
    end
end

facts("Checking MarkovMatrix distribution") do
    d = 5

    pars = Dirichlet{Float64}[Dirichlet(rand(d)) for i in 1:d]
    parmat = hcat([rand(Dirichlet(d, 1)) for i in 1:d]...)
    parvec = [rand(Dirichlet(d, 1)) for i in 1:d]

    x = MarkovMatrix(pars)

    context("Inner constructor") do
        @fact isa(x, MarkovMatrix) --> true
        @fact isa(x, ContinuousMatrixDistribution) --> true
        @fact isa(x, Distribution) --> true
        @fact x.cols --> pars
    end

    context("Inner constructor consistency checks") do
        @fact_throws ErrorException MarkovMatrix(pars[2:end])
    end

    context("Outer constructor") do
        x = MarkovMatrix(parmat)
        @fact isa(x, MarkovMatrix) --> true

        y = MarkovMatrix(parvec)
        @fact isa(y, MarkovMatrix) --> true
    end

    context("Outer constructor consistency checks") do
        @fact_throws ErrorException MarkovMatrix(parmat[2:end, :])
    end

    context("Conversions") do
        pars32 = Dirichlet{Float32}[Dirichlet(rand(5)) for _ in 1:d]
        parvec32 = [rand(Float32, d) for _ in 1:d]
        parmat32 = rand(Float32, d, d)
        x32 = MarkovMatrix(pars32)
        @fact isa(convert(MarkovMatrix{Float64}, pars32), MarkovMatrix{Float64}) --> true
        @fact isa(convert(MarkovMatrix{Float64}, parvec32), MarkovMatrix{Float64}) --> true
        @fact isa(convert(MarkovMatrix{Float64}, parmat32), MarkovMatrix{Float64}) --> true
        @fact isa(convert(MarkovMatrix{Float64}, x32), MarkovMatrix{Float64}) --> true
    end

    context("Check basic interface") do
        @fact nstates(x) --> d
        @fact size(x) --> (d, d)
        @fact length(x) --> d^2
    end

    context("Check mean") do
        m = mean(x)
        @fact size(m) --> (d, d)
        @fact sum(m, 1) --> roughly(ones(1, d))
    end

    context("Check Elog") do
        El = Elog(x)
        @fact size(El) --> (d, d)
    end

    context("Check sampling") do
        A = rand(x)
        @fact size(A) --> (d, d)
        @fact sum(A, 1) --> roughly(ones(1, d))
    end

    context("Check logpdf") do
        # draw a matrix
        A = rand(x)
        lpdf = logpdf(x, A)

        ll = sum([logpdf(x.cols[i], A[:, i]) for i in 1:nstates(x)])
        @fact lpdf --> roughly(ll)
    end

    context("Check entropy") do
        H = sum([entropy(c) for c in x.cols])
        @fact entropy(x) --> roughly(H)
    end

    context("Check natural parameters") do
        nats = naturals(x)
        pars = naturals_to_params(nats, typeof(x))

        params = hcat([c.alpha for c in x.cols]...)

        @fact length(nats) --> 1
        @fact size(nats[1]) --> size(x)
        @fact length(pars) --> 1
        @fact size(pars[1]) --> size(x)
        @fact pars[1] --> roughly(params)
    end

    context("Check params") do
        pars = hcat([c.alpha for c in x.cols]...)
        @fact params(x) --> (pars,)
    end
end

facts("Checking forward-backward algorithm.") do
    # set up some contants for test dataset
    T = 500  # time points
    U = 50  # number of observations at each time
    K = 5  # states
    dt = 1  # time step

    # make transition probabilities
    # we assume that the transition matrix A is left stochastic
    pars = Dirichlet(K, 1)
    A = Array{Float64}(K, K)
    π0 = rand(pars)
    for k in 1:K
        A[:, k] = rand(pars)
    end

    # draw a chain of states
    z = zeros(K, T)
    z_ints = Array{Int}(T)
    init_state = rand(Categorical(π0))
    z[init_state, 1] = 1
    z_ints[1] = init_state
    for t in 2:T
        newstate = rand(Categorical(A * z[:, t - 1]))
        z[newstate, t] = 1
        z_ints[t] = newstate
    end

    # observation model: Poisson
    λ = 10 * collect(1:K)
    N = Array{Int}(U, T)
    for t in 1:T
        N[:, t] = rand(Poisson(λ[z_ints[t]]), U)
    end

    # logpdf of observations, conditioned on state
    logψ = Array{Float64}(K, T)
    for t in 1:T
        for k in 1:K
            logψ[k, t] = 0
            for u in 1:U
                logψ[k, t] += logpdf(Poisson(λ[k]), N[u, t])
            end
        end
    end

    # run forward-backward algorithm
    γ, logZ, Ξ = forwardbackward(π0, A, map(exp, logψ))

    context("Checking validity of test data") do
        @fact sum(π0) --> roughly(1)
        @fact sum(A, 1) --> roughly(ones(1, K))
        @fact length(find(z)) --> T
        @fact sum(z .!= 0, 1) --> ones(1, T)
        @fact logψ .≤ 0 --> all
    end

    context("Check return shapes.") do
        @fact size(γ) --> (K, T)
        @fact sum(γ, 1) --> roughly(ones(1, T))
        @fact size(logZ) --> ()
        @fact size(Ξ) --> (K, K, T - 1)
        @fact sum(Ξ, [1, 2]) --> roughly(ones(1, 1, T - 1))
    end

    context("Rescaling logψ compensated by logZ") do
        offset = randn(1, T)
        logψ_r = logψ .+ offset
        γ_r, logZ_r, Ξ_r = forwardbackward(π0, A, map(exp, logψ_r))

        @fact γ --> roughly(γ_r)
        @fact logZ_r --> roughly(logZ + sum(offset))
        @fact Ξ_r --> roughly(Ξ)
    end

    context("Rescaling π0 compensated by logZ") do
        offset = maximum(log(π0))
        π0_r = exp(log(π0) + offset)
        γ_r, logZ_r, Ξ_r = forwardbackward(π0_r, A, map(exp, logψ))

        @fact γ --> roughly(γ_r)
        @fact logZ_r --> roughly(logZ + offset)
        @fact Ξ_r --> roughly(Ξ)
    end

    context("Rescaling A compensated by logZ") do
        offset = maximum(log(A))
        A_r = exp(log(A) + offset)
        γ_r, logZ_r, Ξ_r = forwardbackward(π0, A_r, map(exp, logψ))

        @fact γ --> roughly(γ_r)
        @fact logZ_r --> roughly(logZ + (T - 1) * offset)
        @fact Ξ_r --> roughly(Ξ)
    end

    context("Test entropy positive") do
        # H = E_q[-log q] logZ - LL ≥ 0
        emission_piece = sum(γ .* logψ)
        initial_piece = sum(γ[:, 1] .* log(π0))
        transition_piece = 0.
        for t in 1:(T - 1)
            transition_piece += sum(Ξ[:, :, t] .* log(A))
        end

        LL = emission_piece + initial_piece + transition_piece

        @fact (logZ ≥ LL || logZ ≈ LL) --> true
    end

    context("Test entropy positive (subadditive parameters)") do
        # same test for positive entropy, but with subadditive A and π0
        logπ0_r = log(π0) - rand(K)
        logA_r = log(A) - rand(size(A))

        # run forward-backward algorithm
        γ, logZ, Ξ = forwardbackward(exp(logπ0_r), exp(logA_r), map(exp, logψ))

        # calculate log likelihood
        emission_piece = sum(γ .* logψ)
        initial_piece = sum(γ[:, 1] .* logπ0_r)
        transition_piece = 0.
        for t in 1:(T - 1)
            transition_piece += sum(Ξ[:, :, t] .* logA_r)
        end

        LL = emission_piece + initial_piece + transition_piece

        @fact (logZ ≥ LL || logZ ≈ LL) --> true
    end

    context("Test accuracy of inference") do
        tmin = 10
        @fact γ[:, tmin:end] --> roughly(z[:, tmin:end], atol=1e-3)
    end
end


text:
################################################################################
#
#	CONCRETE TYPE
#
#   Bandstructure{path} <: AbstractBandstructure{path}
#   --> path is the path in k-space
#
#   FILE CONTAINS
#       - concrete struct definition
#        - interface implementation
#
################################################################################

mutable struct Bandstructure{P,H} <: AbstractBandstructure{P,H}

    # the path along which the band structure is calcualted
    path :: P
    # the hamiltonian
    h :: H

    # bands for each segment
    # bands[i] gives all bands of segment i
    # bands[i][j] gives all energy values for band j of segment i
    # bands[i][j][k] gives the energy value at kpoint index k of band j in segment i
    bands :: Array{Array{Array{Float64, 1}, 1}, 1}

end

# export the type
export Bandstructure

################################################################################
#
#	INTERFACING / ACCESSING BANDSTRUCTURES
#	(functions have to be overwritten by concrete types)
#
################################################################################

# obtaining the hamiltonian
function hamiltonian(
            bs :: Bandstructure{P,H}
        ) :: H where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # return the field
    return bs.h
end

# export the function
export hamiltonian

# obtaining the path
function path(
            bs :: Bandstructure{P,H}
        ) :: P where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # return the field
    return bs.path
end

# export the function
export path

# obtaining the energy values
function energies(
            bs :: Bandstructure{P,H}
        ) :: Array{Array{Array{Float64, 1}, 1}, 1} where {RP, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # return the energy array
    return bs.bands
end

# recalculate the band structure (energy values)
function recalculate!(
            bs :: Bandstructure{P,H}
            ;
            resolution :: Union{Int64, Vector{Int64}} = 100
        ) where {D,RP<:AbstractReciprocalPoint{D}, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # calculate the number of segments
    N_segments = length(path(bs)) - 1
    # parse the segment resolution
    seg_resolution = (zeros(Int64,N_segments).+1).*resolution
    if typeof(resolution)==Vector{Int64}
        if length(resolution) == N_segments
            seg_resolution = resolution
        else
            error("given resolution \"$(resolution)\" has not enough components ($(N_segments) needed)")
        end
    end

    # calculate number of bands
    N_bands = dim(hamiltonian(bs))

    # make a list of all k values
    k_values = Vector{Vector{Vector{Float64}}}(undef, N_segments)
    # make a list for all energy values
    e_values = Vector{Vector{Vector{Float64}}}(undef, N_segments)

    # fill all k values
    @inbounds for i in 1:N_segments
        # fill the k values with points on the path segment
        k_values[i] = getPointsOnLine(path(bs)[i], path(bs)[i+1], seg_resolution[i])
        # fill energy values with zeros
        e_values[i] = Vector{Float64}[
            zeros(Float64, seg_resolution[i]) for j in 1:N_bands
        ]
    end

    # calculate the energy values
    @inbounds for i in 1:N_segments
        for j in 1:seg_resolution[i]
            # use the k_value j of segment i
            k = k_values[i][j]
            # obtain the matrix
            m = matrixAtK(hamiltonian(bs), k)
            # diagonalize the matrix and save the eigenvalues
            m_eigenvalues = sort(real.(eigvals(m)))
            # put the eigenvalues into the bands
            @simd for b in 1:N_bands
                e_values[i][b][j] = m_eigenvalues[b]
            end
        end
    end

    # set the energy values within band structure
    bs.bands = e_values

    # return nothing
    return nothing
end

# export the function
export recalculate!

# Convinience function if path and hamiltonian are known
function getBandstructure(
            path :: P,
            h    :: H
            ;
            recalculate :: Bool = true
        ) :: Bandstructure{P,H} where {D,RP<:AbstractReciprocalPoint{D}, P<:AbstractReciprocalPath{RP}, L,UC,HB,H<:AbstractHamiltonian{L,UC,HB}}

    # create a new object
    bs = Bandstructure{P,H}(path, h, Vector{Vector{Vector{Float64}}}(undef, length(path)-1))

    # recalculate the numerical energy values
    if recalculate
        recalculate!(bs)
    end

    # return the new object
    return bs
end

# Convinience function if only path, unitcell and the bond hamiltonian is known (uses concrete Hamiltonian)
function getBandstructure(
            uc   :: UC,
            path :: P,
            h    :: HB
            ;
            recalculate :: Bool = true
        ) :: Bandstructure{P,Hamiltonian{L,UC,HB}} where {D,RP<:AbstractReciprocalPoint{D}, P<:AbstractReciprocalPath{RP}, L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B},NS,HB<:AbstractBondHamiltonian{L,NS}}

    # create a new object
    bs = Bandstructure{P,Hamiltonian{L,UC,HB}}(path, Hamiltonian{L,UC,HB}(uc, h), Vector{Vector{Vector{Float64}}}(undef, length(path)-1))

    # recalculate the numerical energy values
    if recalculate
        recalculate!(bs)
    end

    # return the new object
    return bs
end

# Convinience function if only path and unitcell but not even the bond hamiltonian is known (uses concrete Hamiltonian and simple hopping hamiltonian)
function getBandstructure(
            uc   :: UC,
            path :: P
            ;
            recalculate :: Bool = true
        ) :: Bandstructure{P,Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}} where {D,RP<:AbstractReciprocalPoint{D}, P<:AbstractReciprocalPath{RP}, L,S,N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B}}

    # create a new object
    bs = Bandstructure{P,Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}}(path, Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}(uc, getHoppingHamiltonianSimple(uc)), Vector{Vector{Vector{Float64}}}(undef, length(path)-1))

    # recalculate the numerical energy values
    if recalculate
        recalculate!(bs)
    end

    # return the new object
    return bs
end

# Convinience function if only unitcell and points are known
function getBandstructure(
            uc     :: UC,
            points :: Union{Symbol, Tuple{Symbol, Int64}} ...
            ;
            recalculate :: Bool = true
        ) :: Bandstructure{ReciprocalPath{ReciprocalPoint{D}},Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}} where {LS,D,L,S<:AbstractSite{LS,D},N,B<:AbstractBond{L,N},UC<:AbstractUnitcell{S,B}}

    # create a new object
    bs = Bandstructure{ReciprocalPath{ReciprocalPoint{D}},Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}}(getReciprocalPath(uc, points...), Hamiltonian{L,UC,BondHoppingHamiltonianSimple{L}}(uc, getHoppingHamiltonianSimple(uc)), Vector{Vector{Vector{Float64}}}(undef, length(points)-1))

    # recalculate the numerical energy values
    if recalculate
        recalculate!(bs)
    end

    # return the new object
    return bs
end

# export the convenience function
export getBandstructure


text:
"""
    BoseFS{N,M,S} <: SingleComponentFockAddress

Address type that represents a Fock state of `N` spinless bosons in `M` modes
by wrapping a bitstring of type `S <: BitString`.

# Constructors

* `BoseFS{N,M}(bs::BitString)`: Unsafe constructor. Does not check whether the number of
  particles in `bs` is equal to `N`.

* `BoseFS(::BitString)`: Automatically determine `N` and `M`. This constructor is not type
  stable!

* `BoseFS{[N,M,S]}(onr)`: Create `BoseFS{N,M}` from [`onr`](@ref) representation. This is
  efficient as long as at least `N` is provided.

See also: [`SingleComponentFockAddress`](@ref), [`FermiFS`](@ref), [`BitString`](@ref).
"""
struct BoseFS{N,M,S<:BitString} <: SingleComponentFockAddress{N,M}
    bs::S
end

function BoseFS{N,M}(bs::BitString{B}) where {N,M,B}
    # Check for consistency between parameter, but NOT for the correct number of bits.
    N + M - 1 == B || throw(ArgumentError("type parameter mismatch"))
    return BoseFS{N,M,typeof(bs)}(bs)
end

function BoseFS(bs::BitString{B}) where B
    N = count_ones(bs)
    M = B - N + 1
    return BoseFS{N,M}(bs)
end

@inline function BoseFS{N,M,S}(
    onr::Union{SVector{M},NTuple{M}}
) where {N,M,S<:BitString{<:Any,1}}
    @boundscheck sum(onr) == N || error("invalid ONR")
    T = chunk_type(S)
    result = zero(T)
    for i in M:-1:1
        curr_occnum = T(onr[i])
        result <<= curr_occnum + T(1)
        result |= one(T) << curr_occnum - T(1)
    end
    return BoseFS{N,M,S}(S(SVector(result)))
end

@inline function BoseFS{N,M,S}(onr::Union{SVector{M},NTuple{M}}) where {N,M,S<:BitString}
    @boundscheck sum(onr) == N || error("invalid ONR")
    K = num_chunks(S)
    result = zeros(MVector{K,UInt64})
    offset = 0
    bits_left = chunk_bits(S, K)
    i = 1
    j = K
    while true
        # Write number to result
        curr_occnum = onr[i]
        while curr_occnum > 0
            x = min(curr_occnum, bits_left)
            mask = (one(UInt64) << x - 1) << offset
            @inbounds result[j] |= mask
            bits_left -= x
            offset += x
            curr_occnum -= x

            if bits_left == 0
                j -= 1
                offset = 0
                bits_left = chunk_bits(S, j)
            end
        end
        offset += 1
        bits_left -= 1

        if bits_left == 0
            j -= 1
            offset = 0
            bits_left = chunk_bits(S, j)
        end
        i += 1
        i > M && break
    end
    return BoseFS{N,M}(S(SVector(result)))
end
function BoseFS{N,M}(onr::Union{AbstractVector,Tuple}) where {N,M}
    S = typeof(BitString{N + M - 1}(0))
    return BoseFS{N,M,S}(SVector{M}(onr))
end
function BoseFS{N}(onr::Union{SVector{M},NTuple{M}}) where {N,M}
    return BoseFS{N,M}(onr)
end
function BoseFS(onr::Union{AbstractVector,Tuple})
    M = length(onr)
    N = sum(onr)
    return BoseFS{N,M}(onr)
end

function Base.show(io::IO, b::BoseFS{N,M,S}) where {N,M,S}
    print(io, "BoseFS{$N,$M}(", tuple(onr(b)...), ")")
end
Base.bitstring(b::BoseFS) = bitstring(b.bs)

Base.isless(a::BoseFS, b::BoseFS) = isless(a.bs, b.bs)
Base.hash(bba::BoseFS,  h::UInt) = hash(bba.bs, h)
Base.:(==)(a::BoseFS, b::BoseFS) = a.bs == b.bs

"""
    near_uniform_onr(N, M) -> onr::SVector{M,Int}

Create occupation number representation `onr` distributing `N` particles in `M`
modes in a close-to-uniform fashion with each mode filled with at least
`N ÷ M` particles and at most with `N ÷ M + 1` particles.
"""
function near_uniform_onr(n::Number, m::Number)
    return near_uniform_onr(Val(n),Val(m))
end
function near_uniform_onr(::Val{N}, ::Val{M}) where {N, M}
    fillingfactor, extras = divrem(N, M)
    # startonr = fill(fillingfactor,M)
    startonr = fillingfactor * @MVector ones(Int,M)
    startonr[1:extras] .+= 1
    return SVector{M}(startonr)
end

"""
    near_uniform(BoseFS{N,M})
    near_uniform(BoseFS{N,M,S}) -> bfs::BoseFS{N,M,S}

Create bosonic Fock state with near uniform occupation number of `M` modes with
a total of `N` particles. Specifying the bit address type `S` is optional.

# Examples
```jldoctest
julia> near_uniform(BoseFS{7,5,BitString{14}})
BoseFS{7,5}((2, 2, 1, 1, 1))

julia> near_uniform(BoseFS{7,5})
BoseFS{7,5}((2, 2, 1, 1, 1))
```
"""
function near_uniform(::Type{<:BoseFS{N,M}}) where {N,M}
    return BoseFS{N,M}(near_uniform_onr(Val(N),Val(M)))
end
near_uniform(b::AbstractFockAddress) = near_uniform(typeof(b))

@deprecate nearUniform near_uniform

"""
    onr(bs)

Compute and return the occupation number representation of the bit string
address `bs` as an `SVector{M,Int32}`, where `M` is the number of modes.
"""
onr(bba::BoseFS) = SVector(m_onr(bba))

"""
    m_onr(bs)

Compute and return the occupation number representation of the bit string
address `bs` as an `MVector{M,Int32}`, where `M` is the number of modes.
"""
@inline m_onr(bba::BoseFS) = m_onr(Val(num_chunks(bba.bs)), bba)

# Version specialized for single-chunk addresses.
@inline function m_onr(::Val{1}, bba::BoseFS{N,M}) where {N,M}
    result = zeros(MVector{M,Int32})
    address = bba.bs
    for mode in 1:M
        bosons = Int32(trailing_ones(address))
        @inbounds result[mode] = bosons
        address >>>= (bosons + 1) % UInt
        iszero(address) && break
    end
    return result
end

# Version specialized for multi-chunk addresses. This is quite a bit faster for large
# addresses.
@inline function m_onr(::Val{K}, bba::BoseFS{N,M}) where {K,N,M}
    B = num_bits(bba.bs)
    result = zeros(MVector{M,Int32})
    address = bba.bs
    mode = 1
    i = K
    while true
        chunk = chunks(address)[i]
        bits_left = chunk_bits(address, i)
        while !iszero(chunk)
            bosons = trailing_ones(chunk)
            @inbounds result[mode] += unsafe_trunc(Int32, bosons)
            chunk >>>= bosons % UInt
            empty_modes = trailing_zeros(chunk)
            mode += empty_modes
            chunk >>>= empty_modes % UInt
            bits_left -= bosons + empty_modes
        end
        i == 1 && break
        i -= 1
        mode += bits_left
    end
    return result
end

function num_occupied_modes(b::BoseFS{<:Any,<:Any,S}) where S
    return num_occupied_modes(Val(num_chunks(S)), b)
end
# For vacuum state
function num_occupied_modes(b::BoseFS{0})
    return 0
end

@inline function num_occupied_modes(::Val{1}, b::BoseFS)
    chunk = b.bs.chunks[1]
    result = 0
    while true
        chunk >>= (trailing_zeros(chunk) % UInt)
        chunk >>= (trailing_ones(chunk) % UInt)
        result += 1
        iszero(chunk) && break
    end
    return result
end

@inline function num_occupied_modes(_, b::BoseFS)
    # This version is faster than using the occupied_mode iterator
    address = b.bs
    result = 0
    K = num_chunks(address)
    last_mask = UInt64(1) << 63 # = 0b100000...
    prev_top_bit = false
    # This loop compiles away for address<:BSAdd*
    for i in K:-1:1
        chunk = chunks(address)[i]
        # This part handles modes that span across chunk boundaries.
        # If the previous top bit and the current bottom bit are both 1, we have to subtract
        # 1 from the result or the mode will be counted twice.
        result -= (chunk & prev_top_bit) % Int
        prev_top_bit = (chunk & last_mask) > 0
        while !iszero(chunk)
            chunk >>>= trailing_zeros(chunk)
            chunk >>>= trailing_ones(chunk)
            result += 1
        end
    end
    return result
end

"""
    BoseFSIndex

Struct used for indexing and performing [`excitation`](@ref)s on a [`BoseFS`](@ref).

## Fields:

* `occnum`: the occupation number.
* `mode`: the index of the mode.
* `offset`: the bit offset of the mode.
"""
struct BoseFSIndex<:FieldVector{3,Int}
    occnum::Int
    mode::Int
    offset::Int
end

function Base.show(io::IO, i::BoseFSIndex)
    @unpack occnum, mode, offset = i
    print(io, "BoseFSIndex(occnum=$occnum, mode=$mode, offset=$offset)")
end
Base.show(io::IO, ::MIME"text/plain", i::BoseFSIndex) = show(io, i)

"""
    BoseOccupiedModes{C,S<:BoseFS}
Iterator for occupied modes. `C` is the number of chunks. See [`occupied_modes`](@ref).
"""
struct BoseOccupiedModes{C,S<:BoseFS}
    address::S
end

function occupied_modes(b::BoseFS{<:Any,<:Any,S}) where {S}
    return BoseOccupiedModes{num_chunks(S),typeof(b)}(b)
end

function is_occupied(::BoseFS, i::BoseFSIndex)
    return i.occnum > 0
end

Base.length(o::BoseOccupiedModes) = num_occupied_modes(o.address)
Base.eltype(::BoseOccupiedModes) = BoseFSIndex

# Single chunk versions are simpler.
@inline function Base.iterate(osi::BoseOccupiedModes{1})
    chunk = osi.address.bs.chunks[1]
    empty_modes = trailing_zeros(chunk)
    return iterate(
        osi, (chunk >> (empty_modes % UInt), empty_modes, 1 + empty_modes)
    )
end
@inline function Base.iterate(osi::BoseOccupiedModes{1}, (chunk, bit, mode))
    if iszero(chunk)
        return nothing
    else
        bosons = trailing_ones(chunk)
        chunk >>>= (bosons % UInt)
        empty_modes = trailing_zeros(chunk)
        chunk >>>= (empty_modes % UInt)
        next_bit = bit + bosons + empty_modes
        next_mode = mode + empty_modes
        return BoseFSIndex(bosons, mode, bit), (chunk, next_bit, next_mode)
    end
end

# Multi-chunk version
@inline function Base.iterate(osi::BoseOccupiedModes)
    bitstring = osi.address.bs
    i = num_chunks(bitstring)
    chunk = chunks(bitstring)[i]
    bits_left = chunk_bits(bitstring, i)
    mode = 1
    return iterate(osi, (i, chunk, bits_left, mode))
end
@inline function Base.iterate(osi::BoseOccupiedModes, (i, chunk, bits_left, mode))
    i < 1 && return nothing
    bitstring = osi.address.bs
    S = typeof(bitstring)
    bit_position = 0

    # Remove and count trailing zeros.
    empty_modes = min(trailing_zeros(chunk), bits_left)
    chunk >>>= empty_modes % UInt
    bits_left -= empty_modes
    mode += empty_modes
    while bits_left < 1
        i -= 1
        i < 1 && return nothing
        @inbounds chunk = chunks(bitstring)[i]
        bits_left = chunk_bits(S, i)
        empty_modes = min(bits_left, trailing_zeros(chunk))
        mode += empty_modes
        bits_left -= empty_modes
        chunk >>>= empty_modes % UInt
    end

    bit_position = chunk_bits(S, i) - bits_left + 64 * (num_chunks(bitstring) - i)

    # Remove and count trailing ones.
    result = 0
    bosons = trailing_ones(chunk)
    bits_left -= bosons
    chunk >>>= bosons % UInt
    result += bosons
    while bits_left < 1
        i -= 1
        i < 1 && break
        @inbounds chunk = chunks(bitstring)[i]
        bits_left = chunk_bits(S, i)

        bosons = trailing_ones(chunk)
        bits_left -= bosons
        result += bosons
        chunk >>>= bosons % UInt
    end
    return BoseFSIndex(result, mode, bit_position), (i, chunk, bits_left, mode)
end

function find_mode(b::BoseFS, index)
    last_occnum = last_mode = last_offset = 0
    for (occnum, mode, offset) in occupied_modes(b)
        dist = index - mode
        if dist == 0
            return BoseFSIndex(occnum, index, offset)
        elseif dist < 0
            return BoseFSIndex(0, index, offset + dist)
        end
        last_occnum = occnum
        last_mode = mode
        last_offset = offset
    end
    offset = last_offset + last_occnum + index - last_mode
    return BoseFSIndex(0, index, offset)
end
# Multiple in a single pass
function find_mode(b::BoseFS, indices::NTuple{N}) where {N}
    # Idea: find permutation, then use the permutation to find indices in order even though
    # they are not sorted.
    perm = sortperm(SVector(indices))
    # perm_i is the index in permutation and goes from 1:N.
    perm_i = 1
    # curr_i points to indices and result
    curr_i = perm[1]
    # index is the current index we are looking for.
    index = indices[curr_i]

    result = ntuple(_ -> BoseFSIndex(0, 0, 0), Val(N))
    last_occnum = last_mode = last_offset = 0
    @inbounds for (occnum, mode, offset) in occupied_modes(b)
        dist = index - mode
        # While loop handles duplicate entries in indices.
        while dist ≤ 0
            if dist == 0
                @set! result[curr_i] = BoseFSIndex(occnum, mode, offset)
            else
                @set! result[curr_i] = BoseFSIndex(0, index, offset + dist)
            end
            perm_i += 1
            perm_i > N && return result
            curr_i = perm[perm_i]
            index = indices[curr_i]
            dist = index - mode
        end
        last_occnum = occnum
        last_mode = mode
        last_offset = offset
    end
    # Now we have to find all indices that appear after the last occupied site.
    # While true because we break out of the loop early anyway.
    @inbounds while true
        offset = last_offset + last_occnum + index - last_mode
        @set! result[curr_i] = BoseFSIndex(0, index, offset)
        perm_i += 1
        perm_i > N && return result
        curr_i = perm[perm_i]
        index = indices[curr_i]
    end
    return result # not reached
end

function find_occupied_mode(b::BoseFS, index::Integer, n=1)
    for (occnum, mode, offset) in occupied_modes(b)
        index -= occnum ≥ n
        if index == 0
            return BoseFSIndex(occnum, mode, offset)
        end
    end
    return BoseFSIndex(0, 0, 0)
end

function _move_particle(b::BoseFS, from, to)
    if to == from
        return b
    elseif to < from
        return typeof(b)(partial_left_shift(b.bs, to, from))
    else
        return typeof(b)(partial_right_shift(b.bs, from, to - 1))
    end
end

###
### Multiple excitation stuff
###
# Fix offsets that changed after performing a move.
@inline function _fix_offset(pair, index::BoseFSIndex)
    fst, snd = pair[1], pair[2]
    if fst.offset < snd.offset
        return @set index.offset += fst.offset < index.offset ≤ snd.offset
    else
        return @set index.offset -= fst.offset > index.offset > snd.offset
    end
end
_fix_offset(pair) = Base.Fix1(_fix_offset, pair)

# Move multiple particles. This does not care about values, so it performs moves in an
# arbitrary order (from left to right in pairs).
@inline function _move_particles(b::BoseFS, (c,)::NTuple{1}, (d,)::NTuple{1})
    return _move_particle(b, d.offset, c.offset)
end
@inline function _move_particles(b::BoseFS, (c, cs...), (d, ds...))
    b = _move_particle(b, d.offset, c.offset)
    fix = _fix_offset(c => d)
    b = _move_particles(b, map(fix, cs), map(fix, ds))
    return b
end

# Apply destruction operator to BoseFSIndex.
@inline _destroy(d, index) = @set index.occnum -= (d.mode == index.mode)
@inline _destroy(d) = Base.Fix1(_destroy, d)
# Apply creation operator to BoseFSIndex.
@inline _create(c, index) = @set index.occnum += (c.mode == index.mode)
@inline _create(c) = Base.Fix1(_create, c)

# Compute the value of an excitation. Starts by applying all destruction operators, and
# then applying all creation operators. The operators must be given in reverse order.
# Will return 0 if move is illegal.
@inline _compute_value(::Tuple{}, ::Tuple{}) = 1
@inline function _compute_value((c, cs...), ::Tuple{})
    return _compute_value(map(_create(c), cs), ()) * (c.occnum + 1)
end
@inline function _compute_value(creations, (d, ds...))
    return _compute_value(map(_destroy(d), creations), map(_destroy(d), ds)) * d.occnum
end

function excitation(b::BoseFS, creations::NTuple{N}, destructions::NTuple{N}) where N
    # We start by computing the value. This is where the check if the move is even legal
    # is done.
    creations_rev = reverse(creations)
    value = _compute_value(creations_rev, reverse(destructions))
    if iszero(value)
        return b, 0.0
    else
        # Now that we know the value and that the move is legal, we can apply the moves
        # without worrying about doing something weird.
        return _move_particles(b, creations_rev, destructions), √value
    end
end


text:
#
# Some small helpfers to hande (the product manifold of) manifold-valued images,
# namely \(\mathcal M^{n,m}\) of an manifold-valued array M
#
export constructImageGraph

"""
    constructImageGraph(img,type)

construct the graph that can be used within the algorithms to model forward
differences
# Arguments
* `img` – the manifold-valued (MPoint-entries) image
* `type` – the type of graph, here `firstOrderDifference` and `secondOrderDifference`.
 """
function constructImageGraph(data::Array{T}, graphtype::String)::Array{Tuple} where {T<:MPoint}
  if graphtype == "firstOrderDifference"
    numdims = ndims(data)
    dims = size(data)
    newnumel = 0
    for i=1:numdims
      ithDiffs = 1
      for j=1:numdims
        if i==j
          ithDiffs = ithDiffs*(dims[j]-1)
        else
          ithDiffs = ithDiffs*dims[j]
        end
      end
      newnumel = newnumel + ithDiffs
    end
    edges = Array{Tuple{Int,Int}}(undef, newnumel)
    # fill with differences
    n=1
    for i=1:numdims
      # generate forward differences along ith dimension
      for index in eachindex(data) #for each index (linear)
        ind = CartesianIndices(dims)[index] # get index (sub)
        if ind[i] < dims[i]     # forward difference in ith direction possible
          ind2 = [Tuple(ind)...]       # extract into array
          ind2[i] = ind2[i] + 1 # generate neighbor
          edges[n] = ( (LinearIndices(dims))[Tuple(ind)...],(LinearIndices(dims))[Tuple(ind2)...] )
          n=n+1
        end
      end
    end
  elseif graphtype == "secondOrderDifference"
    numdims = ndims(data)
    dims = size(data)
    newnumel = 0
    for i=1:numdims
      ithDiffs = 1
      for j=1:numdims
        if i==j
          ithDiffs = ithDiffs*(dims[j]-2)
        else
          ithDiffs = ithDiffs*dims[j]
        end
      end
      newnumel = newnumel + ithDiffs
    end
    edges = Array{Tuple{Int,Int,Int}}(undef,newnumel)
    # fill with differences
    n=1
    for i=1:numdims
      # generate forward differences along ith dimension
      for index in eachindex(data) #for each index (linear)
        ind = CartesianIndices(dims)[index] # get index (sub)
        if (ind[i] < dims[i]) && (1<ind[i])     # left and right neighbot in ith direction within array
          ind2 = [Tuple(ind)...]       # extract into array
          ind2[i] = ind2[i] + 1 # generate neighbor
          ind3 = [Tuple(ind)...]       # extract into array
          ind3[i] = ind3[i] - 1 # generate neighbor
          edges[n] = ( (LinearIndices(dims))[Tuple(ind3)...], (LinearIndices(dims))[Tuple(ind)...], (LinearIndices(dims))[Tuple(ind2)...] )
          n=n+1
        end
      end
    end
  else # end of type checks
      throw( ErrorException(" constructImageGraph – unknown graph type: $graphtype " ) )
  end
  return edges
end


text:
mutable struct NeuralNetwork
    layers::AbstractArray{AbstractLayer,1}
    cost_function::AbstractCostFunction
    n_layers::Int64 # cache the number of layers
    function NeuralNetwork(
        layers_::AbstractArray{AbstractLayer,1},
        cost_function_::AbstractCostFunction,
    )
        new(copy(layers_), cost_function_, size(layers_)[1])
    end
end

function fitting(
    nn::NeuralNetwork,
    X_::AbstractArray{Float64,2},
    t::AbstractArray{Float64,2},
    learning_rate::Float64,
    policy = "Adam",
)
    # X = [x_1 x_2 ,,, x_N]
    # t = [t_1 t_2 ,,, t_N]
    X = copy(X_)
    for layer in nn.layers
        X = forward_propagation(layer, X)
    end

    typeof(X) == typeof(t) || error("fitting: size does not match")

    diff = X - t

    for layer in reverse(nn.layers)
        diff = backward_propagation(layer, diff, learning_rate, policy)
    end

end

function predict(nn::NeuralNetwork, X_::AbstractArray{Float64,2})
    # X = [x_1 x_2 ,,, x_N]
    X = copy(X_)
    for layer in nn.layers
        X = forward_propagation(layer, X)
    end

    return X
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
# TODO should drop_metadata be delete_metadata?

const modules = Module[]
const GLOBAL_METADATA    = gensym(:metadata)

"""
    GlobalMetadata <: AbstractDict{UInt,Dict{Symbol,Any}}

Stores metadata for instances of types at the module level. It has restricted support
for dictionary methods to ensure that references aren't unintentionally
"""
struct GlobalMetadata <: AbstractDict{UInt,MDict}
    data::IdDict{UInt,MDict}

    function GlobalMetadata(m::Module)
        if !isdefined(m, GLOBAL_METADATA)
            Core.eval(m, :(const $GLOBAL_METADATA = $(new(IdDict{UInt,Dict{Symbol,Any}}()))))
            push!(modules, m)
        end
        return metadata(m)
    end
end

data(m::GlobalMetadata) = getfield(m, :data)

Base.get(m::GlobalMetadata, k::UInt, @nospecialize(default)) = get(data(m), k, default)

Base.get!(m::GlobalMetadata, k::UInt, @nospecialize(default)) = get!(data(m), k, default)

Base.isempty(m::GlobalMetadata) = isempty(data(m))

Base.getindex(x::GlobalMetadata, k::UInt) = getindex(data(x), k)

Base.setindex!(x::GlobalMetadata, v::MDict, k::UInt) = setindex!(data(x), v, k)

Base.length(m::GlobalMetadata) = length(data(m))

Base.keys(m::GlobalMetadata) = keys(data(m))

Base.delete!(m::GlobalMetadata, k::UInt) = delete!(data(m), k)

Base.iterate(m::GlobalMetadata) = iterate(data(m))
Base.iterate(m::GlobalMetadata, state) = iterate(data(m), state)

"""
    NoMetadata

Internal type for the `Metadata` package that indicates the absence of any metadata.
_DO NOT_ store metadata with the value `NoMetadata()`.
"""
struct NoMetadata end

const no_metadata = NoMetadata()

Base.show(io::IO, ::NoMetadata) = print(io, "no_metadata")

Base.haskey(::NoMetadata, @nospecialize(k)) = false
Base.get(::NoMetadata, @nospecialize(k), default) = default
Base.getindex(::NoMetadata, @nospecialize(k)) = no_metadata
Base.getproperty(::NoMetadata, ::Symbol) = no_metadata

"""
    MetaStruct(p, m)

Binds a parent instance (`p`) to some metadata (`m`). `MetaStruct` is the generic type
constructed when `attach_metadata(p, m)` is called.

See also: [`attach_metadata`](@ref), [`attach_eachmeta`](@ref)
"""
struct MetaStruct{P,M}
    parent::P
    metadata::M
end

Base.parent(m::MetaStruct) = getfield(m, :parent)
ArrayInterface.parent_type(::Type{MetaStruct{P,M}}) where {P,M} = P

Base.eltype(::Type{T}) where {T<:MetaStruct} = eltype(parent_type(T))

function Base.show(io::IO, ::MIME"text/plain", x::MetaStruct)
    print(io, "attach_metadata($(parent(x)), ::$(metadata_type(x)))\n")
    Metadata.metadata_summary(io, x)
end

Base.copy(x::MetaStruct) = propagate_metadata(x, deepcopy(parent(x)))

@_define_function_no_prop(Base,  ==, MetaStruct, MetaStruct)
@_define_function_no_prop_first(Base,  ==, MetaStruct, Any)
@_define_function_no_prop_last(Base,  ==, Any, MetaStruct)
@_define_function_no_prop_first(Base,  ==, MetaStruct, Missing)
@_define_function_no_prop_last(Base,  ==, Missing, MetaStruct)
@_define_function_no_prop_first(Base,  ==, MetaStruct, WeakRef)
@_define_function_no_prop_last(Base,  ==, WeakRef, MetaStruct)

"""
    MetaUnitRange(x::AbstractUnitRange, meta)

Type for storing metadata alongside a anything that is subtype of `AbstractUnitRange`.

## Examples

```jldoctest
julia> using Metadata

julia> Metadata.MetaUnitRange(1:2, (m1 =1, m2=[1, 2]))
attach_metadata(1:2, ::NamedTuple{(:m1, :m2), Tuple{Int64, Vector{Int64}}})
  • metadata:
     m1 = 1
     m2 = [1, 2]

```
"""
struct MetaUnitRange{T,P<:AbstractRange{T},M} <: AbstractUnitRange{T}
    parent::P
    metadata::M

    function MetaUnitRange{T,P,M}(p::P, m::M) where {T,P,M}
        if known_step(P) == oneunit(T)
            return new{T,P,M}(p, m)
        else
            throw(ArgumentError("step must be 1, got $(step(p))"))
        end
    end

    function MetaUnitRange{T}(p::AbstractRange, m) where {T}
        if eltype(p) <: T
            return MetaUnitRange{T,typeof(p),typeof(m)}(p, m)
        else
            return MetaUnitRange{T}(AbstractUnitRange{T}(p), m)
        end
    end

    MetaUnitRange(p::AbstractRange, m) = MetaUnitRange{eltype(p)}(p, m)
end

Base.parent(r::MetaUnitRange) = getfield(r, :parent)
ArrayInterface.parent_type(::Type{<:MetaUnitRange{<:Any,P,<:Any}}) where {P} = P

@_define_function_no_prop(Base, first, MetaUnitRange)
@_define_function_no_prop(Base, step, MetaUnitRange)
@_define_function_no_prop(Base, last, MetaUnitRange)
@_define_function_no_prop(Base, length, MetaUnitRange)

Base.@propagate_inbounds Base.getindex(r::MetaUnitRange, i::Integer) = parent(r)[i]

ArrayInterface.known_first(::Type{T}) where {T<:MetaUnitRange} = known_first(parent_type(T))

ArrayInterface.known_last(::Type{T}) where {T<:MetaUnitRange} = known_last(parent_type(T))

@propagate_inbounds Base.getindex(r::MetaUnitRange, i) = propagate_metadata(r, parent(r)[i])

function Base.show(io::IO, m::MIME"text/plain", x::MetaUnitRange)
    if haskey(io, :compact)
        show(io, parent(x))
    else
        print(io, "attach_metadata(")
        print(io, parent(x))
        print(io, ", ", Metadata.showarg_metadata(x), ")\n")
        Metadata.metadata_summary(io, x)
    end
end

@propagate_inbounds function Base.getindex(r::MetaUnitRange, s::StepRange{T}) where T<:Integer
    return propagate_metadata(r, getindex(parent(r), s))
end
@propagate_inbounds function Base.getindex(r::MetaUnitRange, s::AbstractUnitRange{T}) where {T<:Integer}
    return propagate_metadata(r, getindex(parent(r), s))
end

Base.getindex(r::MetaUnitRange, ::Colon) = copy(r)

Base.copy(x::MetaUnitRange) = copy_metadata(x, copy(parent(x)))

"""
    MetaIO(io, meta)

Type for storing metadata alongside subtypes of `IO`.
"""
struct MetaIO{T<:IO,M} <: IO
    parent::T
    metadata::M
end

Base.parent(x::MetaIO) = getfield(x, :parent)
ArrayInterface.parent_type(::Type{T}) where {IOType,T<:MetaIO{IOType}} = IOType

@_define_function_no_prop(Base, isreadonly, MetaIO)
@_define_function_no_prop(Base, isreadable, MetaIO)
@_define_function_no_prop(Base, iswritable, MetaIO)
@_define_function_no_prop(Base, stat, MetaIO)
@_define_function_no_prop(Base, eof, MetaIO)
@_define_function_no_prop(Base, position, MetaIO)
@_define_function_no_prop(Base, close, MetaIO)
@_define_function_no_prop(Base, isopen, MetaIO)
@_define_function_no_prop(Base, ismarked, MetaIO)
@_define_function_no_prop(Base, mark, MetaIO)
@_define_function_no_prop(Base, unmark, MetaIO)
@_define_function_no_prop(Base, reset, MetaIO)
@_define_function_no_prop(Base, seekend, MetaIO)

@_define_function_no_prop_first(Base, skip, MetaIO, Integer)
@_define_function_no_prop_first(Base, seek, MetaIO, Integer)
@_define_function_no_prop_first(Base, read, MetaIO, Integer)
@_define_function_no_prop_first(Base, read!, MetaIO, Ref)
@_define_function_no_prop_first(Base, read!, MetaIO, AbstractArray)
@_define_function_no_prop_first(Base, read!,  MetaIO, Array{UInt8})
@_define_function_no_prop_first(Base, read!,  MetaIO, BitArray)

@_define_function_no_prop_first(Base, write, MetaIO, Array)
@_define_function_no_prop_first(Base, write, MetaIO, AbstractArray)
@_define_function_no_prop_first(Base, write, MetaIO, BitArray)
@_define_function_no_prop_first(Base, write, MetaIO, Base.CodeUnits)
@_define_function_no_prop_first(Base, write, MetaIO, Union{Float16, Float32, Float64, Int128, Int16, Int32, Int64, UInt128, UInt16, UInt32, UInt64})

#Base.read(@nospecialize(s::MetaIO), n::Int) = read(parent(s), n)

function Base.write(@nospecialize(s::MetaIO), x::SubArray{T,N,P,I,L} where L where I where P<:Array) where {T, N}
    return write(parent(s), x)
end

@inline metadata_type(::Type{T}) where {IOT,M,T<:MetaIO{IOT,M}} = M

"""
    MetaArray(parent::AbstractArray, metadata)

Custom `AbstractArray` object to store an `AbstractArray` `parent` as well as
some `metadata`.

## Examples

```jldoctest
julia> using Metadata

julia> Metadata.MetaArray(ones(2,2), metadata=(m1 =1, m2=[1, 2]))
2×2 attach_metadata(::Matrix{Float64}, ::NamedTuple{(:m1, :m2), Tuple{Int64, Vector{Int64}}}
  • metadata:
     m1 = 1
     m2 = [1, 2]
)
 1.0  1.0
 1.0  1.0

```
"""
struct MetaArray{T, N, M, A<:AbstractArray} <: ArrayInterface.AbstractArray2{T, N}
    parent::A
    metadata::M

    MetaArray{T,N,M,A}(a::A, m::M) where {T,N,M,A} = new{T,N,M,A}(a, m)
    MetaArray{T,N,M,A}(a::A, m) where {T,N,M,A} = new{T,N,M,A}(a, M(m))
    MetaArray{T,N,M,A}(a, m::M) where {T,N,M,A} = MetaArray{T,N,M,A}(A(a), m)
    MetaArray{T,N,M,A}(a, m) where {T,N,M,A} = MetaArray{T,N,M,A}(A(a), M(m))
    function MetaArray{T,N,M,A}(a::AbstractArray; metadata=Dict{Symbol,Any}(), kwargs...) where {T,N,M,A}
        return MetaArray{T,N,M,A}(a, _construct_meta(metadata, values(kwargs)))
    end

    function MetaArray{T,N,M,A}(args...; metadata=Dict{Symbol,Any}(), kwargs...) where {T,N,M,A}
        return MetaArray{T,N,M,A}(A(args...); metadata=metadata, kwargs...)
    end

    ###
    ### MetaArray{T,N,M}
    ###
    function MetaArray{T,N,M}(x::AbstractArray, m::M) where {T,N,M}
        if eltype(x) <: T
            return MetaArray{T,N,M,typeof(x)}(x, m)
        else
            return MetaArray{T,N,M}(convert(AbstractArray{T}, x), m)
        end
    end

    ###
    ### MetArray{T,N}
    ###
    MetaArray{T,N}(a::AbstractArray, m::M) where {T,N,M} = MetaArray{T,N,M}(a, m)
    function MetaArray{T,N}(a::AbstractArray; metadata=Dict{Symbol,Any}(), kwargs...) where {T,N}
        return MetaArray{T,N}(a, _construct_meta(metadata, values(kwargs)))
    end
    function MetaArray{T,N}(args...; metadata=Dict{Symbol,Any}(), kwargs...) where {T,N}
        return MetaArray{T,N}(Array{T,N}(args...); metadata=metadata, kwargs...)
    end

    ###
    ### MetArray{T}
    ###
    function MetaArray{T}(args...; metadata=Dict{Symbol,Any}(), kwargs...) where {T}
        return MetaArray{T}(Array{T}(args...); metadata=metadata, kwargs...)
    end
    MetaArray{T}(a::AbstractArray, m::M) where {T,M} = MetaArray{T,ndims(a)}(a, m)
    function MetaArray{T}(a::AbstractArray; metadata=Dict{Symbol,Any}(), kwargs...) where {T}
        return MetaArray{T,ndims(a)}(a; metadata=metadata, kwargs...)
    end

    ###
    ### MetaArray
    ###
    MetaArray(v::AbstractArray{T,N}, m::M) where {T,N,M} = new{T,N,M,typeof(v)}(v, m)
    function MetaArray(a::AbstractArray; metadata=Dict{Symbol,Any}(), kwargs...)
        return MetaArray{eltype(a)}(a; metadata=metadata, kwargs...)
    end
end

ArrayInterface.parent_type(::Type{MetaArray{T,M,N,A}}) where {T,M,N,A} = A

Base.parent(A::MetaArray) = getfield(A, :parent)

Base.copy(A::MetaArray) = copy_metadata(A, copy(parent(A)))

function Base.show(io::IO, ::MIME"text/plain", X::MetaArray)
    summary(io, X)
    isempty(X) && return
    Base.show_circular(io, X) && return

    if !haskey(io, :compact) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if get(io, :limit, false) && eltype(X) === Method
        io = IOContext(io, :limit => false)
    end

    if get(io, :limit, false) && Base.displaysize(io)[1] - 4 <= 0
        return print(io, " …")
    else
        println(io)
    end

    io = IOContext(io, :typeinfo => eltype(X))

    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, parent(X))
end

function Base.similar(x::MetaArray, ::Type{T}, dims::NTuple{N,Int}) where {T,N}
    return Metadata.share_metadata(x, similar(parent(x), T, dims))
end

function Base.similar(
    x::MetaArray,
    ::Type{T},
    dims::Tuple{Union{Integer,OneTo},Vararg{Union{Integer,OneTo}}}
) where {T}

    return Metadata.propagate_metadata(x, similar(parent(x), T, dims))
end
function Base.similar(x::MetaArray, ::Type{T}, dims::Tuple{Integer, Vararg{Integer}}) where {T}
    return Metadata.propagate_metadata(x, similar(parent(x), T, dims))
end

Base.summary(io::IO, x::MetaArray) = Base.showarg(io, x, true)
function Base.showarg(io::IO, x::MetaArray, toplevel)
    if toplevel
        print(io, Base.dims2string(length.(axes(x))), " ")
    end
    print(io, "attach_metadata(")
    Base.showarg(io, parent(x), false)
    print(io, ", ", showarg_metadata(x))
    println(io)
    metadata_summary(io, x)
    print(io, "\n)")
end

function ArrayInterface.defines_strides(::Type{T}) where {T<:MetaArray}
    return ArrayInterface.defines_strides(parent_type(T))
end

@propagate_inbounds function Base.getindex(A::MetaArray{T}, args...) where {T}
    return _getindex(A, getindex(parent(A), args...))
end

_getindex(A::MetaArray{T}, val::T) where {T} = val
_getindex(A::MetaArray{T}, val) where {T} = propagate_metadata(A, val)