text:
module KPP_TKE

using OceanTurb

import ..OceanTurb: oncell, onface
import .KPP: ∂B∂z, isunstable, ωτ, ωb
import .ModularKPP: AbstractModularKPPModel

const nsol = 5
@solution U V T S e

Base.@kwdef struct TKEParameters{T} <: AbstractParameters
        CLz :: T = 0.4    # Dissipation parameter
        CLb :: T = 0.7    # Dissipation parameter
         Cτ :: T = 400.0  # Dissipation parameter

        CDe :: T = 2.0    # Dissipation parameter
        Ke₀ :: T = 1e-6   # Interior diffusivity for salinity

       CK_U :: T = 0.1    # Diffusivity parameter
       CK_T :: T = 0.1    # Diffusivity parameter
       CK_e :: T = 0.1    # Diffusivity parameter

        KU₀ :: T = 1e-6   # Interior viscosity for velocity
        KT₀ :: T = 1e-7   # Interior diffusivity for temperature
        KS₀ :: T = 1e-9   # Interior diffusivity for salinity
end

mutable struct State{T} <: FieldVector{6, T}
    Fu :: T
    Fv :: T
    Fθ :: T
    Fs :: T
    Fb :: T
     h :: T
    function State(T=Float64)
        new{T}(0, 0, 0, 0, 0, 0)
    end
end

"""
    update_state!(model)

Update the top flux conditions and mixing depth for `model`
and store in `model.state`.
"""
function update_state!(m)
    m.state.Fu = getbc(m, m.bcs.U.top)
    m.state.Fv = getbc(m, m.bcs.V.top)
    m.state.Fθ = getbc(m, m.bcs.T.top)
    m.state.Fs = getbc(m, m.bcs.S.top)
    m.state.Fb = m.constants.g * (m.constants.α * m.state.Fθ - m.constants.β * m.state.Fs)
    m.state.h  = ModularKPP.mixing_depth(m)
    return nothing
end

struct Model{KP, HP, SP, S, BC, TS, G, T} <: AbstractModularKPPModel{KP, HP, Nothing, TS, G, T}
      @add_clock_grid_timestepper
        solution :: S
             bcs :: BC
             tke :: TKEParameters{T}
     diffusivity :: KP
     mixingdepth :: HP
        kprofile :: SP
       constants :: Constants{T}
           state :: State{T}
end

function Model(; N=10, L=1.0,
            grid = UniformGrid(N, L),
       constants = Constants(),
             tke = TKEParameters(),
     diffusivity = ModularKPP.LMDDiffusivity(),
     mixingdepth = ModularKPP.LMDMixingDepth(),
        kprofile = ModularKPP.DiffusivityShape(),
         stepper = :ForwardEuler,
    )

    solution = Solution((CellField(grid) for i=1:nsol)...)

    bcs = (
        U = DefaultBoundaryConditions(eltype(grid)),
        V = DefaultBoundaryConditions(eltype(grid)),
        T = DefaultBoundaryConditions(eltype(grid)),
        S = DefaultBoundaryConditions(eltype(grid)),
        e = DefaultBoundaryConditions(eltype(grid))
    )

    Kϕ = (U=KU, V=KV, T=KT, S=KS, e=Ke)
    Rϕ = (U=RU, V=RV, T=RT, S=RS, e=Re)
    eq = Equation(K=Kϕ, R=Rϕ, update=update_state!)
    lhs = OceanTurb.build_lhs(solution)

    timestepper = Timestepper(stepper, eq, solution, lhs)

    return Model(Clock(), grid, timestepper, solution, bcs, tke,
                    diffusivity, mixingdepth, kprofile, constants, State())
end

# Note: to increase readability, we use 'm' to refer to 'model' in function
# definitions below.
#

@inline function mixing_time(m, i)
    return 400.0
    #N = sqrt( max(0, ∂B∂z(m, i)) )
    #@inbounds ũ = sqrt( max(0, m.solution.e[i]) ) + 1e-16
    #return @inbounds min(m.tke.Cτ, m.tke.CLz * abs(m.grid.zf[i]) / ũ, m.tke.CLb / N)
    #return @inbounds min(m.tke.Cτ, m.tke.CLb / N)
end

@inline K_mixing_time(m::AbstractModel{TS, G, T}, i) where {TS, G, T} =
    max(zero(T), mixing_time(m, i) * onface(m.solution.e, i)) + m.tke.KU₀

@inline oncell(f::Function, m, i) = 0.5 * (f(m, i) + f(m, i+1))
@inline onface(f::Function, m, i) = 0.5 * (f(m, i) + f(m, i-1))

@inline ∂B∂z(m, i) = ∂B∂z(m.solution.T, m.solution.S, m.constants.g, m.constants.α,
                            m.constants.β, i)

@inline production(m, i) = KU(m, i) * (∂z(m.solution.U, i)^2 + ∂z(m.solution.V, i)^2)

@inline buoyancy_flux(m, i) = - m.constants.g * (
      m.constants.α * KT(m, i) * ∂z(m.solution.T, i)
    - m.constants.β * KS(m, i) * ∂z(m.solution.S, i)
    )

@inline dissipation(m, i) = @inbounds m.tke.CDe * m.solution.e[i] / mixing_time(m, i)

#
# Equation entry
#

@inline K(m, i) = @inbounds mixing_time(m, i) * onface(m.solution.e, i)

#=
@inline KU(m, i) = m.tke.KU₀ #+ m.tke.CK_U * K(m, i)
@inline KT(m, i) = m.tke.KT₀ #+ m.tke.CK_T * K(m, i)
@inline KS(m, i) = m.tke.KS₀ #+ m.tke.CK_T * K(m, i)
@inline Ke(m, i) = m.tke.Ke₀ #+ m.tke.CK_e * K(m, i)
=#

@inline KU(m, i) = ModularKPP.KU(m, i)
@inline KT(m, i) = ModularKPP.KT(m, i)
@inline KS(m, i) = ModularKPP.KS(m, i)
@inline Ke(m, i) = ModularKPP.KU(m, i)

const KV = KU

@inline RU(m, i) = @inbounds   m.constants.f * m.solution.V[i]
@inline RV(m, i) = @inbounds - m.constants.f * m.solution.U[i]
@inline RT(m, i) = 0
@inline RS(m, i) = 0

@inline Re(m, i) = oncell(production, m, i) + oncell(buoyancy_flux, m, i) - dissipation(m, i)

end # module


text:
#### Simple example of a new Variable in Caesar.
#### This variable is a duplication of Pose2Pose2 in RoME.

"""
    $(SIGNATURES)

A duplication of the Pose2Pose2 factor structure for demonstration purposes.
Note that this a FunctorPairwise factor.
"""
mutable struct ExamplePose2Pose2{T} <: IncrementalInference.FunctorPairwise where {T <: IIF.SamplableBelief}
 z::T
 ExamplePose2Pose2{T}() where {T <: IIF.SamplableBelief} = new{T}()
 ExamplePose2Pose2{T}(z1::T) where {T <: IIF.SamplableBelief} = new{T}(z1)
end

"""
    $(SIGNATURES)

Constructor for the ExamplePose2Pose2 factor structure.
"""
ExamplePose2Pose2(z::T) where {T <: IIF.SamplableBelief} = ExamplePose2Pose2{T}(z)

"""
    $(SIGNATURES)

The sampler for the ExamplePose2Pose2 structure.
This provides N samples from the factor s.
"""
getSample(s::Pose2Pose2{<:IIF.SamplableBelief}, N::Int=1) = (rand(s.z,N), )

"""
    $(SIGNATURES)

???
"""
function (s::Pose2Pose2{<:IIF.SamplableBelief})(res::Array{Float64},
           userdata,
           idx::Int,
           meas::Tuple,
           wxi::Array{Float64,2},
           wxj::Array{Float64,2}  )
 #
 wXjhat = SE2(wxi[1:3,idx])*SE2(meas[1][1:3,idx])
 jXjhat = SE2(wxj[1:3,idx]) \ wXjhat
 se2vee!(res, jXjhat)
 nothing
end

ext: jl


text:
abstract type aggregator end

struct BordaCount <: aggregator end
struct BordaCount <: aggregator end

struct agentPreference
function BordaCount(A)
    nCandid = max(A[1])
    Votes = zeros(nCandid)
    for x in A

    end
end
using Random
A = [randperm(10) for i=1:15]


text:
#=
Benchmarking routines for characters and strings

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md

THIS IS STILL VERY WIP AND HARDCODED!

The variable userdir controls where to save different files downloaded from the internet,
and store the renamed copies (stripped of ASCII lines) for later processing
(I may change that part so that the it is done automatically, and the files are just renamed,
or the files to be processed are picked up from a configuration file)

books = load_books()     # This downloads and processes the books in the gutenbergbooks list
save_books(books)        # This saves them in an output directory (samples) for later benchmarking
tst = checktests()       # Runs a bunch of tests, compares results to the results from String
res = benchdir()         # Runs the benchmark functions on all the documents found in the directory
save_results(fname, res) # Saves the results to the given file

res = load_results(fname) # Loads the results from the given file
dispbench(res)            # Displays the results in a pretty format
=#

const V6_COMPAT = VERSION < v"0.7.0-DEV"

using BenchmarkTools

const test_legacy = false

@static if test_legacy
    using LegacyStrings
end

using StringLiterals
using Strs
import Strs: LineCounts, CharTypes, CharStat, calcstats
import Strs: _LatinStr, _UCS2Str, _UTF32Str, _LatinChr

uninit(T, len) = @static V6_COMPAT ? T(len) : T(undef, len)
create_vector(T, len)  = uninit(Vector{T}, len)

import Base: show

const inppath = "textsamples"
const gutpath = "gutenberg"
const smppath = "samples"

const gutenbergbooks =
    (("files/2600/2600-0",        "English"), # War & Peace, some other languages in quotes
     ("files/1400/1400-0",        "English"), # Great Expectations, uses Unicode quotes
     ("files/42286/42286-0",      "Hungarian"),
     #("files/8119/8119-0",        "Polish"), # couldn't get this to load correctly
     ("files/32941/32941-0",      "Japanese"),
     ("files/24264/24264-0",      "Chinese"),
     ("files/40687/40687-0",      "Telugu"), # Third most spoken in India, official
     ("files/50513/50513-0",      "French"),
#     ("files/43007/43007-0",      "Arabic"),
     ("cache/epub/38496/pg38496", "Portuguese"),
     ("cache/epub/2000/pg2000",   "Spanish"), # Don Quijote
     ("cache/epub/48750/pg48750", "Swedish"),
     ("cache/epub/48322/pg48322", "German"),
     )

const downloadedbooks =
    (("LYSAIa GORA DIeVICh'Ia - SIeRGIeI GOLOVAChIoV.txt", "Russian"),
     )

getdefdir(dir)::String = dir === nothing ? homedir() : dir

function filter_lines(lines)
    out = Vector{String}()
    sizehint!(out, length(lines))
    # Eliminate initial lines, empty lines, trailing lines
    checkbeg = true
    for l in lines
        if sizeof(l) > 41 && starts_with(l, "***") && ends_with(l, "***") &&
            occurs_in(" PROJECT GUTENBERG EBOOK", l) &&
            occurs_in(checkbeg ? "START OF TH" : "END OF TH", l)
            checkbeg || break # Found "end of"
            checkbeg = false
        else
            push!(out, l)
        end
    end
    out
end

"""
Load books from Project Gutenberg site, removing lines added at beginning and end that
are not part of the book, as much as possible
"""
function load_gutenberg!(books, list, dict, gutenbergdir)
    for (nam, lang) in list
        cnt = get(dict, lang, 0)
        dict[lang] = cnt + 1
        outnam = cnt == 0 ? "$lang.txt" : "$lang-$cnt.txt"
        lname = joinpath(gutenbergdir, outnam)
        download(joinpath("http://www.gutenberg.org/", nam * ".txt"), lname)
        println("Saved to: ", lname)
        push!(books, (outnam, filter_lines(readlines(lname))))
    end
    books
end

function remove_empty(lines)
    # Eliminate empty lines
    len = length(lines)
    out = Vector{String}()
    sizehint!(out, len)
    for l in lines
        is_empty(l) || push!(out, l)
    end
    out
end

"""
load_books(; dir=nothing)

Loads a set of books from a local directory, and downloads a set of books from Project Gutenberg
Returns them as a dictionary with names -> vectors of strings
"""
function load_books(; dir::Any=nothing)
    defdir = getdefdir(dir)
    inputdir = joinpath(defdir, inppath)
    dict = Dict{String,Int}()
    books = Vector{Tuple{String, Vector{String}}}()
    for (nam, lang) in downloadedbooks
        cnt = get(dict, lang, 0)
        dict[lang] = cnt + 1
        outnam = cnt == 0 ? "$lang.txt" : "$lang-$cnt.txt"
        push!(books, (outnam, readlines(joinpath(inputdir, nam))))
    end
    load_gutenberg!(books, gutenbergbooks, dict, joinpath(defdir, gutpath))
end

"""
save_books(books; dir=nothing)

Saves the collection of downloaded books into the given directory, in a "samples" subdirectory.
If the directory is not set, it will default to the user's home directory
"""
function save_books(books; dir::Any=nothing)
    sampledir = joinpath(getdefdir(dir), smppath)
    for (nam, book) in books
        outnam = joinpath(sampledir, nam)
        open(outnam, "w") do io
            for lin in book
                println(io, lin)
            end
            println("Saved $nam")
        end
    end
end

show(io::IO, cnt::LineCounts) =
    pr"\(io)\%10d(cnt.bytes)\%12.3f(cnt.bytes/cnt.chars)"

function show(io::IO, s::CharTypes)
    pr"\(io)\%10d(s.ascii)\%10d(s.latin)\%10d(s.utf2byte)\%10d(s.ucs2)"
    pr"\(io)\%10d(s.utf32)\%10d(s.surr)\%10d(s.invalid)"
end

function show(io::IO, v::Tuple{String,CharStat})
    s = v[2]
    pr_ul(io, "File name              Lines     Chars   Avg C/L     Empty")
    pr_ul(io, "       Min       Max   MaxType\n")
    pr"\(io)\%-18s(v[1])\%10d(s.num)\%10d(s.len)\%10.3f(s.len/s.num)\%10d(s.empty)"
    pr"\(io)\%10d(s.minlen)\%10d(s.maxlen)\%10s(string(enctyp(s.maxtyp)))\n\n"
    pr_ul(io, "Character Types:       ")
    pr_ul(io, "ASCII    Latin1    2-Byte      UCS2     UTF32 Surrogate   Invalid\n")
    pr"\(io)Total characters: \(s.chars)\n"
    pr"\(io)Lines with > 0:   \(s.lines)\n"
end

_stdout() = @static V6_COMPAT ? STDOUT : stdout
@static if V6_COMPAT
    const pwc = print_with_color
else
    pwc(c, io, str) = printstyled(io, str; color = c)
end
pwc(c, l) = pwc(c, _stdout(), l)

pr_ul(io, l) = pwc(:underline, io, l)
pr_ul(l)     = pwc(:underline, l)

print_size_ratio(io, val) =
    pwc(val < .95 ? :green : val > 1.05 ? :red : :normal, io, f"\%6.3f(val)")
print_time_ratio(io, val) =
    pwc(val < .95 ? :red : val > 1.05 ? :green : :normal, io, f"\%10.3f(val)")

print_size_ratio(val) = print_size_ratio(_stdout(), val)
print_time_ratio(val) = print_time_ratio(_stdout(), val)

print_ratio(val) = pwc(val < .95 ? :red : val > 1.05 ? :green : :normal, f"\%12.3f(val)")
print_ratio_rev(val) = pwc(val < .95 ? :green : val > 1.05 ? :red : :normal, f"\%12.3f(val)")

function display_results(io, xres)
    # (fname, stats, sizes, res)
    (fname, stats, sizes, selstat, selsiz, res) = xres
    show(io, (fname, stats))
    show(io, (fname, selstat))
    t1 = res[1][3]
    numchars = selstat.len
    pr"\(io)\n\n\%-12.12s(res[1][1])\%6.3f(sizes[1]/stats.len)"
    for rn in res[i]
        pr"\(io)\n\%-12.12s(rn[1])\%6.3f(sizes[i]/stats.len)"
    end
    println(io)
    for i = 1:length(res)
        pr"\(io)\%10.3f(t1[i][3]/numchars)"
        for rn in res[i]
            print_time_ratio(io, t1[i][3]/rn[3])
        end
    end
    println(io)
    # pwc(:yellow, io, f"\%10.3f(rn[3]/numchars)")
end

function dispres(io, xres)
    # (fname, stats, sizes, res)
    (fname, stats, sizes, selstat, selsiz, res) = xres
    show(io, (fname, stats))
    show(io, (fname, selstat))
    maxlen = 0
    pos = 0
    for i = 1:length(res)
        len = length(res[i][3])
        len > maxlen && (maxlen = len; pos = i)
    end
    rs = res[pos][3]
    nam1 = []
    nam2 = []
    for result in rs
        names = split(result[1], '\n')
        push!(nam1, names[1])
        push!(nam2, length(names) > 1 ? names[2] : "")
    end
    pr"\(io)\nType        B/Char"
    for nam in nam1 ; pr"\(io)\%10s(nam)" ; end
    pr_ul(io, "\n                  ")
    for nam in nam2 ; pr_ul(io, f"\%10s(nam)") ; end
    pr"\nResults:          "
    for results in rs
        pr"\(io)\%10d(results[2])"
    end
    r1 = res[1]
    t1 = r1[3]
    numchars = selstat.len
    pr"\(io)\n\n\%-12.12s(r1[1])\%6.3f(sizes[1]/stats.len)"
    for tim in t1
        pr"\(io)\%10.3f(tim[3]/numchars)"
    end
    for i = 2:length(res)
        rn = res[i]
        pr"\(io)\n\%-12.12s(rn[1])\%6.3f(sizes[i]/stats.len)"
        tn = rn[3]
        minres = min(length(t1), length(tn))
        for i = 1:minres
            print_time_ratio(io, t1[i][3]/tn[i][3])
        end
        for i = minres+1:length(tn)
            pwc(:yellow, io, f"\%10.3f(tn[i][3]/numchars)")
        end
    end
    println(io)
end

const divline = string(repeat('#', 100),'\n','\f')

function dispbench(io, totres)
    for res in totres[1]
        dispres(io, res)
        print(io, divline)
    end
end

dispbench(totres) = dispbench(_stdout(), totres)

function display_benchmark(io, totres)
    for res in totres[1]
        display_results(io, res)
        print(io, divline)
    end
end

display_benchmark(totres) = displaybench(_stdout(), totres)

function calchash(lines::Vector{<:AbstractString})
    hashval = 0%UInt
    for text in lines
        hashval = hash(text, hashval)
    end
    hashval%UInt32 # Just show bottom 32-bits, for displaying better
end

function countlength(lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += length(text)
    end
    cnt
end

function dolowercase(lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += lowercase(text) !== text
    end
    cnt
end

function douppercase(lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += uppercase(text) !== text
    end
    cnt
end

@inline function iteratechars(text::AbstractString)
    cnt = 0
    for ch in text
        cnt += is_digit(ch)
    end
    cnt
end

@inline function iteratecps(text::AbstractString)
    cnt = 0
    for ch in codepoints(text)
        cnt += is_digit(ch)
    end
    cnt
end

@inline function iteratecus(text::AbstractString)
    cnt = 0
    for ch in codeunits(text)
        cnt += is_digit(ch)
    end
    cnt
end

function countchars(lines::Vector{T}) where {T<:AbstractString}
    cnt = 0
    for text in lines
        cnt += iteratechars(text)
    end
    cnt
end

function countcps(lines::Vector{T}) where {T<:AbstractString}
    cnt = 0
    for text in lines
        cnt += iteratecps(text)
    end
    cnt
end

function countsize(lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += sizeof(text)
    end
    cnt
end

@inline function countcodeunits(text::AbstractString)
    cnt = 0
    @inbounds for cu in codeunits(text)
        cnt += cu
    end
    cnt
end

function countsize2(lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += countcodeunits(text)
    end
    cnt
end

function checkstr(fun, lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += fun(text)
    end
    cnt
end

function checkuni(lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += sizeof(UniStr(text))
    end
    cnt
end

function checkcvt(::Type{T}, lines) where {T}
    cnt = 0
    for text in lines
        cnt += sizeof(convert(T, text))
    end
    cnt
end

function checktrans(::Type{T}, lines) where {T}
    #println("checktrans($T, $(eltype(lines))")
    cnt = 0
    for text in lines
        cnt += sizeof(transcode(T, text))
    end
    cnt
end

function checktext(fun, lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        is_empty(text) || fun(text)
        cnt += 1
    end
    cnt
end

@inline function iteratenextind(text)
    cnt = 0
    len = ncodeunits(text)
    pos = 1
    while pos <= len
        pos = nextind(text, pos)
        cnt += 1
    end
    cnt
end

@inline function iteratefunchars(fun, text)
    cnt = 0
    for ch in text
        cnt += fun(ch)
    end
    cnt
end

@inline function iteratefuncps(fun, text)
    cnt = 0
    for ch in codepoints(text)
        cnt += fun(ch)
    end
    cnt
end

@inline function iteratefuncus(fun, text)
    cnt = 0
    for ch in codeunits(text)
        cnt += fun(ch)
    end
    cnt
end

function checkchars(fun, lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += iteratefunchars(fun, text)
    end
    cnt
end

function checkcp(fun, lines::Vector{<:AbstractString})
    cnt = 0
    for text in lines
        cnt += iteratefuncps(fun, text)
    end
    cnt
end

checkjoin(lines) = sizeof(join(lines))

@inline function sumchars(text)
    t = 0
    for ch in text
        t += UInt32(ch)
    end
    t
end
@inline function sumcp(text)
    t = 0
    for ch in codepoints(text)
        t += UInt32(ch)
    end
    t
end
@inline function sumcu(text)
    t = 0
    for ch in codeunits(text)
        t += UInt32(ch)
    end
    t
end

function sumcharvals(lines::Vector{<:AbstractString})
    t = 0
    for text in lines
        t += sumchars(text)
    end
    t
end
function sumcodepnts(lines::Vector{<:AbstractString})
    t = 0
    for text in lines
        t += sumcp(text)
    end
    t
end
function sumcodeunits(lines::Vector{<:AbstractString})
    t = 0
    for text in lines
        t += sumcu(text)
    end
    t
end

searchres(x::UnitRange) = x.start
searchres(x::Int) = x

function searchlines(lines, v, rev=false)
    t = 0
    if rev
    for text in lines
        t += searchres(fnd(Last, v, text))
    end
    else
    for text in lines
        t += searchres(fnd(First, v, text))
    end
    end
    t
end

function searchlinescvt(lines, v, rev=false)
    t = 0
    if rev
    for text in lines
        t += searchres(fnd(Last, v, utf8(text)))
    end
    else
    for text in lines
        t += searchres(fnd(First, v, utf8(text)))
    end
    end
    t
end

function searchlines(lines::Vector{UniStr}, v, rev=false)
    t = 0
    if rev
    for text in lines
        if typeof(text) === ASCIIStr
            t += searchres(fnd(Last, v, text::ASCIIStr))
        elseif typeof(text) === _LatinStr
            t += searchres(fnd(Last, v, text::_LatinStr))
        elseif typeof(text) === _UCS2Str
            t += searchres(fnd(Last, v, text::_UCS2Str))
        else
            t += searchres(fnd(Last, v, text::_UTF32Str))
        end
    end
    else
    for text in lines
        if typeof(text) === ASCIIStr
            t += searchres(fnd(First, v, text::ASCIIStr))
        elseif typeof(text) === _LatinStr
            t += searchres(fnd(First, v, text::_LatinStr))
        elseif typeof(text) === _UCS2Str
            t += searchres(fnd(First, v, text::_UCS2Str))
        else
            t += searchres(fnd(First, v, text::_UTF32Str))
        end
    end
    end
    t
end

# maybe change all this to be table driven and use @eval!

searchchar(lines::Vector{String}, d=false) =
    searchlines(lines, '\ufffd', d)
searchchar(lines::Vector{T}, d=false) where {T} =
    searchlines(lines, eltype(T)(0xfffd), d)
searchchar(lines::Vector{T}, d=false) where {T<:Union{ASCIIStr,LatinStr,_LatinStr}} =
    searchlines(lines, eltype(T)(0x1a), d)

searchstr(lines::Vector{String}, d=false) = searchlines(lines, "this is a test", d)
searchstr(lines::Vector{T}, d=false) where {T} = searchlines(lines, T("this is a test"), d)

searchreg(lines::Vector{<:AbstractString}, d=false) =
    searchlinescvt(lines, r"this is a test", d)
searchreg(lines::Vector{T}, d=false) where {T<:Str{<:Strs.Regex_CSEs}} =
    searchlines(lines, r"this is a test", d)
searchreg(lines::Vector{String}, d=false) = searchlines(lines, r"this is a test", d)

rsearchchar(lines) = searchchar(lines, true)
rsearchstr(lines) = searchstr(lines, true)

# normalize
# isassigned

sklength(str) = length(str)
@inline function _length(s, i::Int, n::Int, c::Int)
    i < n || return c
    @inbounds b = codeunit(s, i)
    @inbounds while true
        while true
            (i += 1) ≤ n || return c
            0xc0 ≤ b ≤ 0xf7 && break
            b = codeunit(s, i)
        end
        l = b
        b = codeunit(s, i) # cont byte 1
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xe0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 2
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xf0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 3
        c -= (b & 0xc0 == 0x80)
    end
end
sklength(s::T) where {T<:AbstractString} = _length(s, 1, ncodeunits(s), ncodeunits(s))

_is_surrogate_trail(c::Unsigned) = ((c & ~0x003ff) == 0xdc00)
_is_valid_continuation(c) = ((c & 0xc0) == 0x80)

function oldlength(str::T) where {T<:AbstractString}
    dat = Vector{UInt8}(str)
    cnt = 0
    for i = 1:sizeof(dat)
        @inbounds cnt += !_is_valid_continuation(dat[i])
    end
    cnt
end
function oldlength(s::UTF16Str)
    dat = reinterpret(UInt16, s.data)
    len = length(dat)
    len == 0 && return 0
    cnt = 0
    for i = 1:len
        @inbounds cnt += !_is_surrogate_trail(dat[i])
    end
    cnt
end

repeat1(str)  = repeat(str, 1)
repeat10(str) = repeat(str, 10)
repeat1c(str)  = @inbounds repeat(str[1], 1)
repeat10c(str) = @inbounds repeat(str[1], 10)

countsklength(l)  = checkstr(sklength, l)
countoldlength(l) = checkstr(oldlength, l)
checktextwidth(l) = checkstr(text_width, l)

checkrepeat1(l)   = checktext(repeat1, l)
checkrepeat10(l)  = checktext(repeat10, l)
checkrepeat1c(l)  = checktext(repeat1c, l)
checkrepeat10c(l) = checktext(repeat10c, l)
checkreverse(l)   = checktext(reverse, l)

checknextind(l) = checkstr(iteratenextind, l)
countchars(l)   = checkstr(iteratechars, l)
countcps(l)     = checkstr(iteratecps, l)
countcus(l)     = checkstr(iteratecus, l)

validstr(l)     = checkstr(is_valid, l)
asciistr(l)     = checkstr(is_ascii, l)

validchars(l)   = checkchars(is_valid, l)
asciichars(l)   = checkchars(is_ascii, l)

checkvalid(l)   = checkcp(is_valid,     l)
checkascii(l)   = checkcp(is_ascii,     l)
checkcntrl(l)   = checkcp(is_control,   l)
checkdigit(l)   = checkcp(is_digit,     l)
checkxdigit(l)  = checkcp(is_hex_digit, l)
checklower(l)   = checkcp(is_lowercase, l)
checkupper(l)   = checkcp(is_uppercase, l)
checkalpha(l)   = checkcp(is_alpha,     l)
checknumeric(l) = checkcp(is_numeric,   l)
checkspace(l)   = checkcp(is_space,     l)
checkalnum(l)   = checkcp(is_alphanumeric, l)
checkprint(l)   = checkcp(is_printable,    l)
checkpunct(l)   = checkcp(is_punctuation,  l)
checkgraph(l)   = checkcp(is_graphic,      l)

function mintime(f, lines)
    m = typemax(UInt)
    for i=1:10
        ns = time_ns()
        f(lines)
        t = time_ns() - ns
        t < m && (m = t)
    end
    m
end

function wrap(f, lines, io, cnts::LineCounts, t, msg, fast=true, basetime=0%UInt)
    tst = ""
    try
        res = f(lines)
        tim = fast ? mintime(f, lines) : round(UInt64, 1e9 * @belapsed ($f)($lines))
        push!(t, (msg, res, tim))
        pr"\(io)\%-22s(replace(msg, '\n' => ' ')*':') \%12d(res)"
        pr"\(io) \%12.3f(tim/1000000)"
        pr"\(io) \%12.3f(tim/cnts.lines)"
        pr"\(io) \%12.3f(tim/cnts.chars)"
        pr"\(io) \%12.3f(tim/cnts.bytes)"
        basetime != 0 && print_ratio(basetime/tim)
        pr"\(io)\n"
    catch ex
        typeof(ex) == InterruptException && rethrow()
        push!(t, (msg, 0, 0%UInt))
        println(io, sprint(showerror, ex, catch_backtrace()))
    end
end

function mintime(::Type{T}, f, lines) where {T}
    m = typemax(UInt)
    for i=1:10
        ns = time_ns()
        f(T, lines)
        t = time_ns() - ns
        t < m && (m = t)
    end
    m
end

function wrap(::Type{T}, f, lines, io, cnts::LineCounts, t, msg,
              fast=true, basetime=0%UInt) where {T}
    tst = ""
    try
        res = f(T, lines)
        tim = fast ? mintime(T, f, lines) : round(UInt64, 1e9 * @belapsed ($f)($lines))
        push!(t, (msg, res, tim))
        pr"\(io)\%-22s(replace(msg, '\n' => ' ')*':') \%12d(res)"
        pr"\(io) \%12.3f(tim/1000000)"
        pr"\(io) \%12.3f(tim/cnts.lines)"
        pr"\(io) \%12.3f(tim/cnts.chars)"
        pr"\(io) \%12.3f(tim/cnts.bytes)"
        basetime != 0 && print_ratio(basetime/tim)
        pr"\(io)\n"
    catch ex
        typeof(ex) == InterruptException && rethrow()
        push!(t, (msg, 0, 0%UInt))
        println(io, sprint(showerror, ex, catch_backtrace()))
    end
end

# For each set of lines: save number of characters (should be the same for all)
# For each test (currently just counting via Char iteration, length)

const tests =
    (
     (countsize,   "sizeof"),
     (countlength, "length"),
     (calchash,     "hash\nstring"),
     (checknextind, "nextind\nchars"),
     (checkjoin,    "join\nlines"),
     (checkreverse, "reverse"),
     (checkrepeat1,  "repeat 1\nstring"),
     (checkrepeat10,  "repeat 10\nstring"),
     (searchstr,    "search\nstring"),
     (searchchar,    "search\nchar"),
     (searchreg,     "search\nregex"),
#     (rsearchstr,    "rsearch\nstring"),
#     (rsearchchar,    "rsearch\nchar"),
#     (checkrepeat1c,  "repeat 1\nchar"),
#     (checkrepeat10c,  "repeat 10\nchar"),
#    (countsklength,  "length\nSK"),
#    (countoldlength, "length\nOld"),
     (countchars,   "iteration\nChar"),
     (sumcharvals,  "sum\nchar vals"),
     (asciistr,     "isascii\nstring"),
     (validstr,     "isvalid\nstring"),
     (checkascii,   "isascii\nchars"),
     (checkvalid,   "isvalid\nchars"),
     (checklower,   "islower\nchars"),
     (checkalnum,   "isalnum\nchars"),
     #=
     (checkcntrl,   "iscntrl\nchars"),
     (checkupper,   "isupper\nchars"),
     (checkalpha,   "isalpha\nchars"),
     (checkspace,   "isspace\nchars"),
     (checkprint,   "isprint\nchars"),
     (checkpunct,   "ispunct\nchars"),
     (checkgraph,   "isgraph\nchars"),
     (checkdigit,   "isdigit\nchars"),
     (checkxdigit,  "isxdigit\nchars"),
     (dolowercase,  "lowercase\nstring"),
     =#
     (douppercase,  "uppercase\nstring"),
     (checktextwidth, "textwidth\nstring"),
     #    (dotitlecase,  "titlecase\nstring"),
    )

function select_lines(lines::Vector{<:AbstractString}; num=1000)
    out = Vector{Int}()
    sizehint!(out, num)
    len = length(lines)
    i = len>>1
    j = len-1
    while length(out) < num
        is_empty(lines[i]) || push!(out, i)
        is_empty(lines[j]) || push!(out, j)
        i += 1 > len && break
        j -= 1 < 1   && break
    end
    out
end

function testperf(lines::Vector{T}, io, cnts, docnam, basetime, fast) where {T<:AbstractString}
    # Test performance
    pr_ul(io, f"""\%-22s(docnam) \%12s("Result") \%12s("ms total") """)
    pr_ul(io, f"""\%12s("ns/line") \%12s("ns/char") \%12s("ns/byte")\n""")
    t = []
    # Just run everything once
    for (tst, nam) in tests
        try
            tst(lines)
        catch ex
            typeof(ex) == InterruptException && rethrow()
            str = f"Failed test \(nam): \(sprint(showerror, ex, catch_backtrace()))"
            println("File: ", docnam, ": ", str)
            return docnam, 0, t, str
        end
    end
    pos = 0
    for (tst, nam) in tests
        wrap(tst, lines, io, cnts, t, nam, fast,
             basetime == nothing ? 0%UInt : basetime[pos += 1][3])
    end
    pr"\n"
    docnam, cnts.bytes, t, ""
end

function testcvt(::Type{S}, vec, io, numlines, numchars, siz, docnam,
                 basetime) where {S}
    tr = checktrans
    cv = checkcvt
    #println(typeof(vec))

    # Test performance
    pr_ul(io, f"""\%-22s(docnam) \%12s("Result") \%12s("ms total") """)
    pr_ul(io, f"""\%12s("ns/line") \%12s("ns/char") \%12s("ns/byte")\n""")
    t = []
    # Just run everything once
    try
        if S === String
            tr(UInt16, vec[1])
            tr(UInt32, vec[1])
            tr(String, vec[2])
            tr(String, vec[3])
        else
            cv(S, vec[1])
            cv(S, vec[4])
            cv(S, vec[5])
        end
    catch ex
        typeof(ex) == InterruptException && rethrow()
        str = f"Failed test: \(sprint(showerror, ex, catch_backtrace()))"
        println("File: ", docnam, ": ", str)
        return docnam, 0, t, str
    end
    if S === String
        wrap(UInt16, tr, vec[1], io, LineCounts(numlines, numchars, siz[1]), t,
             "transcode to UTF16", true, 0%UInt)
        wrap(UInt32, tr, vec[1], io, LineCounts(numlines, numchars, siz[1]), t,
             "transcode to UTF32", true, 0%UInt)
        wrap(String, tr, vec[2], io, LineCounts(numlines, numchars, siz[2]), t,
             "transcode UInt16", true, 0%UInt)
        wrap(String, tr, vec[3], io, LineCounts(numlines, numchars, siz[3]), t,
             "transcode UInt32", true, 0%UInt)
    else
        wrap(S, cv, vec[1], io, LineCounts(numlines, numchars, siz[1]), t,
             "convert String", true, basetime[1 + (codeunit(S) === UInt32)][3])
        wrap(S, cv, vec[4], io, LineCounts(numlines, numchars, siz[2]), t,
             "convert UInt16", true, basetime[3][3])
        wrap(S, cv, vec[5], io, LineCounts(numlines, numchars, siz[3]), t,
             "convert UInt32", true, basetime[4][3])
    end
    pr"\n"
    docnam, siz, t, ""
end

function sumsizes1(lines)
    cnt = 0
    for text in lines
        cnt += sizeof(text)
    end
    cnt
end

sumsizes2(lines) = sum(sizeof, lines)

function benchsizes(lines)
    f = sumsizes1
    t1 = @belapsed sumsizes1($lines)
    t2 = @belapsed sumsizes2($lines)
    t3 = @belapsed ($f)($lines)
    t4 = @belapsed sum(sizeof, $lines)
    t1, t2, t3, t4
end

enctyp(t) = t == 1 ? ASCIIStr : t == 2 ? LatinStr : t < 5 ? UCS2Str : UTF32Str

function encode_lines(list, lines)
    numlines = length(lines)
    enc = []
    push!(enc, lines)
    for i = 2:length(list)
        T = list[i]
        vec = create_vector(T, numlines)
        @inbounds for j = 1:numlines
            vec[j] = convert(T, lines[j])
        end
        push!(enc, vec)
    end
    enc
end

function checkboolchar(fun, lines)
    res = []
    for text in lines
        len = length(text)
        bv = uninit(BitVector, len)
        for (i, ch) in enumerate(text)
            bv[i] = fun(ch)
        end
        push!(res, bv)
    end
    res
end

function checkline(::Type{Bool}, fun, lines)
    len = length(lines)
    bv = uninit(BitVector, len)
    for (i, text) in enumerate(lines)
        bv[i] = fun(text)
    end
    bv
end

function checkline(::Type{T}, fun, lines) where {T}
    res = T[]
    for text in lines
        push!(res, fun(text))
    end
    res
end

function runcheckline(::Type{T}, lines, list) where {T}
    totresults = []
    for fun in list
        push!(totresults, checkline(T === Integer ? typeof(fun(lines[1])) : T, fun, lines))
    end
    totresults
end

function runcheckchar(lines, list)
    totresults = []
    for fun in list
        push!(totresults, checkboolchar(fun, lines))
    end
    totresults
end

function runcheckcu(lines, list)
    totresults = []
    for fun in list
        res = []
        for text in lines
            push!(res, [fun(cu) for cu in codeunits(text)])
        end
        push!(totresults, res)
    end
    totresults
end

striplist(list) = replace(string(list), "Base.Unicode." => "")
function comparetestline(lines, results, list, displist)
    pr"  Lines:     \(displist)"
    diff = []
    for (i, fun) in enumerate(list)
        fundiff = []
        funres = results[i]
        try
            for (j, text) in enumerate(lines)
                res = fun(text)
                res == funres[j] ||
                    push!(fundiff, typeof(res) == Bool ? (j, text) : (j, text, res, funres[j]))
            end
            is_empty(fundiff) || push!(diff, (i, fun, fundiff))
        catch ex
            typeof(ex) == InterruptException && rethrow()
            pr"Failed test \(fun): \(sprint(showerror, ex, catch_backtrace()))"
        end
    end
    if is_empty(diff)
        pwc(:green, "\r\u2714\n")
    else
        pwc(:red, "\e[s\rX\e[u")
        io = IOBuffer()
        try
            print(io, " => ", diff)
        catch ex
            typeof(ex) == InterruptException && rethrow()
            pr"Can't display diffs \(diff[end][2]): \(sprint(showerror, ex, catch_backtrace()))"
        end
        str = String(take!(io))
        println(str[1:200])
    end
    diff
end

function comparetestchar(lines, results, list, displist)
    pr"  Chars:     \(displist)"
    diff = []
    for (i, fun) in enumerate(list)
        fundiff = []
        lineres = results[i]
        try
            for (j, text) in enumerate(lines)
                chrdiff = []
                chrres = lineres[j]
                for (k, ch) in enumerate(text)
                    res = fun(ch)
                    res == chrres[k] ||
                        push!(chrdiff, typeof(res) == Bool ? (k, ch) : (k, ch, res, chrres[k]))
                end
                is_empty(chrdiff) || push!(fundiff, (j, text, chrdiff))
            end
            is_empty(fundiff) || push!(diff, (i, fun, fundiff))
        catch ex
            typeof(ex) == InterruptException && rethrow()
            pr"Failed test \(fun): \(sprint(showerror, ex, catch_backtrace()))"
        end
    end
    if is_empty(diff)
        pwc(:green, "\r\u2714\n")
    else
        pwc(:red, "\e[s\rX\e[u")
        io = IOBuffer()
        try
            print(io, " => ", diff)
        catch ex
            typeof(ex) == InterruptException && rethrow()
            pr"Can't display diffs \(diff[end][2]): \(sprint(showerror, ex, catch_backtrace()))"
        end
        str = String(take!(io))
        println(str[1:200])
    end
    diff
end

function comparetestcu(lines, results, list, displist)
    pr"  CodeUnits: \(displist)"
    diff = []
    for (i, fun) in enumerate(list)
        fundiff = []
        lineres = results[i]
        try
            for (j, text) in enumerate(lines)
                chrdiff = []
                chrres = lineres[j]
                for (k, ch) in enumerate(codeunits(text))
                    res = fun(ch)
                    res == chrres[k] ||
                        push!(chrdiff, typeof(res) == Bool ? (k, ch) : (k, ch, res, chrres[k]))
                end
                is_empty(chrdiff) || push!(fundiff, (j, text, chrdiff))
            end
            is_empty(fundiff) || push!(diff, (i, fun, fundiff))
        catch ex
            typeof(ex) == InterruptException && rethrow()
            pr"Failed test \(fun): \(sprint(showerror, ex, catch_backtrace()))"
        end
    end
    if is_empty(diff)
        pwc(:green, "\r\u2714\n")
    else
        pwc(:red, "\e[s\rX\e[u")
        io = IOBuffer()
        try
            print(io, " => ", diff)
        catch ex
            typeof(ex) == InterruptException && rethrow()
            pr"Can't display diffs \(diff[end][2]): \(sprint(showerror, ex, catch_backtrace()))"
        end
        str = String(take!(io))
        println(str[1:200])
    end
    diff
end

const testlist =
    (((length, text_width, hash), "length, text_width, hash"),
     ((is_ascii, is_valid), "is_ascii, is_valid"),
     ((lowercase, uppercase, reverse), "lowercase, uppercase, reverse"),
     ((is_ascii, is_valid, is_control, is_lowercase, is_uppercase, is_alpha,
       is_alphanumeric, is_space, is_printable, is_punctuation, is_graphic,
       is_digit, is_hex_digit),
      "is_(ascii|valid|control|lowercase|uppercase|alpha|alphanumeric|space|printable|" *
      "punctuation|graphic|digit|hex_digit)"),
     ((UInt32, ), "UInt32"),
     ((sizeof, ), "sizeof"))

# Use String to calculate a baseline (note, sizeof needs to be checked separately, as only
# UTF8Str should be the same).

function compareall(io, lines, res)
    pr"\(io)\(eltype(lines))\n"
    (comparetestline(lines, res[1], testlist[1]...),
     comparetestline(lines, res[2], testlist[2]...),
     comparetestline(lines, res[3], testlist[3]...),
     comparetestchar(lines, res[4], testlist[4]...),
     eltype(lines) == UTF8Str ? comparetestcu(lines, res[5], testlist[5]...) : [],
     eltype(lines) == UTF8Str ? comparetestline(lines, res[6], testlist[6]...) : [])
end

function checktests(io = _stdout(); dir::Any=nothing)
    totres = []
    totcmp = []
    sampledir = joinpath(getdefdir(dir), smppath)
    for fname in readdir(sampledir)
        lines = readlines(joinpath(sampledir, fname))
        stats = calcstats(lines)
        list = [String, UTF8Str, UTF16Str, UTF32Str, UniStr]
        MT = enctyp(stats.maxtyp)
        MT != UTF32Str && push!(list, MT)
        isdefined(Main, :UTF8String) && push!(list, UTF8String, UTF16String, UTF32String)
        enc = encode_lines(list, lines)
        res = (runcheckline(Integer, lines, testlist[1][1]),
               runcheckline(Bool, lines, testlist[2][1]),
               runcheckline(Any, lines, testlist[3][1]),
               runcheckchar(lines, testlist[4][1]),
               runcheckcu(lines,   testlist[5][1]),
               runcheckline(Int, lines, testlist[6][1]))
        push!(totres, (fname, res))
        cmp = []
        pr"\(io)Checking \(fname):\n"
        for i = 2:length(list)
            push!(cmp, compareall(io, enc[i], res))
        end
        push!(totres, (fname, res))
        push!(totcmp, (fname, cmp))
    end
    totres, totcmp
end

function benchdir(io = _stdout();  dir::Any=nothing, fast::Bool=true)
    totres = []
    totlines = []
    totnames = []
    totsizes = []
    sampledir = joinpath(getdefdir(dir), smppath)
    for fname in readdir(sampledir)
        lines = readlines(joinpath(sampledir, fname))
        stats = calcstats(lines)
        show(io, (fname, stats))
        list = [String, UTF8Str, UTF16Str, UTF32Str, UniStr]
        MT = enctyp(stats.maxtyp)
        MT != UTF32Str && push!(list, MT)
        isdefined(Main, :UTF8String) && push!(list, UTF8String, UTF16String, UTF32String)
        enc = encode_lines(list, lines)
        # Now calculate and display the size statistics for each
        sizes = [sum(sizeof, enclines) for enclines in enc]
        basesize = sizes[1]
        names = String[string(typ) for typ in list]

        pr_ul(io, "\nTypes:      ")
        for nam in names ; pr_ul(io, f"\%12s(nam)") ; end
        pr"\(io)\nBytes:      "
        for siz in sizes ; pr"\(io)\%12d(siz)" ; end
        pr"\(io)\nBytes/Char: "
        for siz in sizes ; pr"\(io)\%12.3f(siz/stats.len)" ; end
        pr"\(io)\nRelative:               "
        for siz in sizes[2:end] ; print_ratio_rev(siz/basesize) ; end
        pr"\n\n"

        # Select the middle 100 non-empty lines of the file for benchmarking
        sellines = select_lines(lines)
        sel = [enclines[sellines] for enclines in enc]
        selsiz = [sum(sizeof, enclines) for enclines in sel]
        selstat = calcstats(sel[1])
        numchars = selstat.len
        numlines = selstat.num

        # Now test the performance for each
        res = create_vector(Any, length(list))
        res[1] = testperf(sel[1], io, LineCounts(numlines, numchars, selsiz[1]),
                          names[1], nothing, fast)
        basetime = res[1][3]
        for i = 2:length(list)
            res[i] = testperf(sel[i], io, LineCounts(numlines, numchars, selsiz[i]),
                              names[i], basetime, fast)
        end

        # Now test the performance of conversions for each
        v16 = [transcode(UInt16, s) for s in sel[1]]
        v32 = [transcode(UInt32, s) for s in sel[1]]
        cvtsel = (sel[1], v16, v32,
                  [convert(Text2Str, v) for v in v16],
                  [convert(Text4Str, v) for v in v32])
        cvtsiz = [sum(sizeof, enclines) for enclines in cvtsel]

        cvtres = create_vector(Any, length(list))
        cvtres[1] = testcvt(list[1], cvtsel, io, numlines, numchars, cvtsiz, names[1], nothing)
        basetime = cvtres[1][3]
        for i = 2:length(list)
            cvtres[i] = testcvt(list[i], cvtsel, io, numlines, numchars, cvtsiz, names[i], basetime)
        end
        push!(totres, (fname, stats, sizes, selstat, selsiz, res, cvtres))

        push!(totnames, names)
        push!(totlines, enc)

        print(divline)
    end
    totres, totlines, totnames, totsizes
end

function save_results(fname, res)
    open(fname, "w") do io
        serialize(io, res)
    end
end

function load_results(fname)
    open(fname) do io
        deserialize(io)
    end
end


text:
struct Publish <: Packet
    header::UInt8
    id::UInt16
    message::Message
end
Publish(message::Message) = Publish(convert(UInt8, PUBLISH) | ((message.dup & 0x1) << 3) | (convert(UInt8, message.qos) << 1) | message.retain, 0x0000, message)
Publish(packet::Publish, id::UInt16) = Publish(packet.header, id, packet.message)

function read(s::IO, flags::UInt8, ::Type{Publish})
    dup = (flags & 0x08) >> 3
    qos = (flags & 0x06) >> 1
    retain = (flags & 0x01)
    topic = mqtt_read(s, String)
    id = 0x0000
    if qos != 0x00
        id = mqtt_read(s, UInt16)
    end
    payload = take!(s)
    return Publish(PUBLISH, id, Message(dup, qos, retain, topic, payload))
end

function write(s::IO, packet::Publish)
    mqtt_write(s, packet.message.topic)
    if packet.message.qos != AT_MOST_ONCE
        mqtt_write(s, packet.id)
    end
    write(s, packet.message.payload)
end

function needs_id(packet::Publish)
    return packet.message.qos != AT_MOST_ONCE
end

function has_id(packet::Publish)
    return packet.message.qos != AT_MOST_ONCE
end

Base.show(io::IO, x::Publish) = print(io, "publish[", ((x.message.qos == AT_MOST_ONCE) ? "" : "id: $(x.id) "), "message: ", x.message, "]")

struct Puback <: Ack
    header::UInt8
    id::UInt16
end
Puback(id) = Puback(PUBACK, id)
Base.show(io::IO, x::Puback) = print(io, "puback[id: ", x.id, "]")

struct Pubrec <: Ack
    header::UInt8
    id::UInt16
end
Pubrec(id) = Pubrec(PUBREC, id)
has_id(packet::Pubrec) = false
Base.show(io::IO, x::Pubrec) = print(io, "pubrec[id: ", x.id, "]")

struct Pubrel <: Ack
    header::UInt8
    id::UInt16
end
Pubrel(id) = Pubrel(convert(UInt8, PUBREL) | 0x02, id)
has_id(packet::Pubrel) = false
Base.show(io::IO, x::Pubrel) = print(io, "pubrel[id: ", x.id, "]")

struct Pubcomp <: Ack
    header::UInt8
    id::UInt16
end
Pubcomp(id) = Pubcomp(PUBCOMP, id)
Base.show(io::IO, x::Pubcomp) = print(io, "pubcomp[id: ", x.id, "]")

max_issues_repo_name: kivaari/jlMQTT

max_forks_repo_name: kivaari/jlMQTT


text:
using Test, DISSEQT.BamReader

# TODO: add Seq tests

# usage
# @test @typestable(expr) == Int
# @test @typestable(expr) == subtypeof(Integer)
macro typestable(expr)
	esc(:(typeof(@inferred($expr))))
end

# a little hack to be able to use --> subtypeof(ParentType)
import Base: ==, show
struct SubTypeOf
	t::Union{DataType,UnionAll}
end
==(x::Union{DataType,UnionAll},y::SubTypeOf) = x <: y.t
==(x::SubTypeOf,y::Union{DataType,UnionAll}) = y==x
show(io::IO, s::SubTypeOf) = print(io,"<:",s.t)

subtypeof(t::Union{DataType,UnionAll}) = SubTypeOf(t)

@testset "CIGAR" begin
	strs = ["1M", "10S268435455M21S", "20S10M1I9M2D10M3X7=10M4P24M10H"]
	cigArrays = Array[[CIGARElement(1<<4 | CIGAR_M)],
	                 [CIGARElement(10<<4 | CIGAR_S), CIGARElement(268435455<<4 | CIGAR_M), CIGARElement(21<<4 | CIGAR_S)],
				     [CIGARElement(20<<4 | CIGAR_S), CIGARElement(10<<4 | CIGAR_M), CIGARElement(1<<4 | CIGAR_I), CIGARElement(9<<4 | CIGAR_M), CIGARElement(2<<4 | CIGAR_D), CIGARElement(10<<4 | CIGAR_M), CIGARElement(3<<4 | CIGAR_X), CIGARElement(7<<4 | CIGAR_Eq), CIGARElement(10<<4 | CIGAR_M), CIGARElement(4<<4 | CIGAR_P), CIGARElement(24<<4 | CIGAR_M), CIGARElement(10<<4 | CIGAR_H)]]
	for (str,cigArray) in zip(strs,cigArrays)
		cig = CIGAR(str)
		@test string(cig) == str
		@test collect(cig) == cigArray
	end
end

function verifymyread(r::BamRead)
	rname = "TheReadName";

	@test refID(r) == 0
	@test @typestable(refID(r)) == Int

	@test pos(r) == 1234
	@test @typestable(pos(r)) == Int

	@test bin_mq_nl(r) == UInt((Int(0x1249)<<16)|(14<<8)|(length(rname)+1))
	@test @typestable(bin_mq_nl(r)) == UInt

	@test flag_nc(r) == UInt(0x600000c)
	@test @typestable(flag_nc(r)) == UInt

	@test l_seq(r) == 94
	@test @typestable(l_seq(r)) == Int

	@test next_refID(r) == -1
	@test @typestable(next_refID(r)) == Int

	@test next_pos(r) == 0
	@test @typestable(next_pos(r)) == Int

	@test tlen(r) == 0
	@test @typestable(tlen(r)) == Int

	@test BamReader.bin(r) == 0x1249
	@test @typestable(BamReader.bin(r)) == Int

	@test mapq(r) == UInt8(14)
	@test @typestable(mapq(r)) == Int

	@test l_read_name(r) == UInt8(length(rname)+1)
	@test @typestable(l_read_name(r)) == Int

	@test flag(r) == 0x600
	@test @typestable(flag(r)) == Int

	@test n_cigar_op(r) == 12
	@test @typestable(n_cigar_op(r)) == Int

	@test read_name(r) == rname
	@test @typestable(read_name(r)) == subtypeof(AbstractString)

	@test ref_name(r) == "TheChrName"
	@test @typestable(ref_name(r)) == subtypeof(AbstractString)

	cig = "20S10M1I9M2D10M3X7=10M4P24M10H"
	@test string(cigar(r)) == cig
	@test @typestable(string(cigar(r))) == subtypeof(AbstractString)

	@test string(copy(cigar(r))) == cig
	@test @typestable(string(copy(cigar(r)))) == subtypeof(AbstractString)

	@test cigar_str(r) == cig
	@test @typestable(cigar_str(r)) == subtypeof(AbstractString)

	cigArray = [CIGARElement(20<<4 | CIGAR_S), CIGARElement(10<<4 | CIGAR_M), CIGARElement(1<<4 | CIGAR_I), CIGARElement(9<<4 | CIGAR_M), CIGARElement(2<<4 | CIGAR_D), CIGARElement(10<<4 | CIGAR_M), CIGARElement(3<<4 | CIGAR_X), CIGARElement(7<<4 | CIGAR_Eq), CIGARElement(10<<4 | CIGAR_M), CIGARElement(4<<4 | CIGAR_P), CIGARElement(24<<4 | CIGAR_M), CIGARElement(10<<4 | CIGAR_H)]
	@test collect(cigar(r)) == cigArray
	@test @typestable(collect(cigar(r))) == Array{CIGARElement,1}

	s = "ACGTAACCGGTTAAACCCGGGTTTAAAACCCCGGGGTTTTAAAAACCCCCGGGGGTTTTTAAAAAACCCCCCGGGGGGTTTTTTGTCACATTAT"
	@test string(seq(r)) == s
	@test @typestable(string(seq(r))) == subtypeof(AbstractString)

	@test string(copy(seq(r))) == s
	@test @typestable(string(copy(seq(r)))) == subtypeof(AbstractString)

	@test seq_str(r) == s
	@test @typestable(seq_str(r)) == subtypeof(AbstractString)

	@test join(map(x->Char(base_char(x)), seq(r))) == s

	rngs = UnitRange{Int}[ 1:length(s), 2:length(s)-1, 1:length(s)-1, 2:length(s)-2, 3:3, 4:4, 1:2, 2:3, 4:87, 1:0, 2:1, length(s):length(s)-1 ]
	for rng in rngs
		@test string(view(seq(r),rng)) == s[rng]
		@test @typestable(string(view(seq(r),rng))) == subtypeof(AbstractString)
	end

	@test qual(r) == collect(UInt8('!'):UInt8('~')).-33
	@test @typestable(qual(r)) == subtypeof(AbstractArray)

	@test hastag(r,"CC") == true
	@test hastag(r,"AA") == false
	@test tagint(r,"NH") == 1
	@test tagint(r,"CP") == 6789
	@test tagstring(r,"CC") == "AnotherChr"

	@test @typestable(hastag(r,"CC")) == Bool
	@test @typestable(tagint(r,"NH")) == Int
	@test @typestable(tagstring(r,"CC")) == subtypeof(AbstractString)

	@test hastag(r,tag"CC") == true
	@test hastag(r,tag"AA") == false
	@test tagint(r,tag"NH") == 1
	@test tagint(r,tag"CP") == 6789
	@test tagstring(r,tag"CC") == "AnotherChr"

	@test @typestable(hastag(r,tag"CC")) == Bool
	@test @typestable(tagint(r,tag"NH")) == Int
	@test @typestable(tagstring(r,tag"CC")) == subtypeof(AbstractString)

	# TODO: test _all_ tag types

	# compare printout with SAM file
	samfromfile = read("BamReader/data/singleread.sam", String)
	samfromfile = split(samfromfile,"\n")[3] # hardcoded that there are two header rows
	@test sam(r) == samfromfile
end

@testset "BamRead" begin
	buf = UInt8[0,1,0,0,0,0,0,0,209,4,0,0,12,14,73,18,12,0,0,6,94,0,0,0,255,255,255,255,255,255,255,255,0,0,0,0,84,104,101,82,101,97,100,78,97,109,101,0,68,1,0,0,160,0,0,0,17,0,0,0,144,0,0,0,34,0,0,0,160,0,0,0,56,0,0,0,119,0,0,0,160,0,0,0,70,0,0,0,128,1,0,0,165,0,0,0,18,72,17,34,68,136,17,18,34,68,72,136,17,17,34,34,68,68,136,136,17,17,18,34,34,68,68,72,136,136,17,17,17,34,34,34,68,68,68,136,136,136,72,33,33,136,24,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,78,72,67,1,67,67,90,65,110,111,116,104,101,114,67,104,114,0,67,80,83,133,26]

	bamFile = BamFile("dummyname", "@HD\tVN1.3\n@SQ\tSN:TheChrName\tLN:1000000\n", [("TheChrName",1000000)], BamIndex())

	# verify header
	@test filename(bamFile) == "dummyname"
	@test header(bamFile) == "@HD\tVN1.3\n@SQ\tSN:TheChrName\tLN:1000000\n"
	@test sequences(bamFile) == [("TheChrName",1000000)]

	rd = BamRead(bamFile)
	BamReader.nextread!(IOBuffer(buf),rd)
	BamReader.compute_offsets!(rd)

	for r in [rd,copy(rd)] # test on both read and copy of read
		verifymyread(r)
	end
end

# test with small BAM
@testset "BAM" begin
	#bamData = UInt8[0,0,0,0,209,4,0,0,12,14,73,18,12,0,71,202,94,0,0,0,255,255,255,255,255,255,255,255,0,0,0,0,84,104,101,82,101,97,100,78,97,109,101,0,68,1,0,0,160,0,0,0,17,0,0,0,144,0,0,0,34,0,0,0,160,0,0,0,56,0,0,0,119,0,0,0,160,0,0,0,70,0,0,0,128,1,0,0,165,0,0,0,18,72,17,34,68,136,17,18,34,68,72,136,17,17,34,34,68,68,136,136,17,17,18,34,34,68,68,72,136,136,17,17,17,34,34,34,68,68,68,136,136,136,72,33,33,136,24,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,78,72,67,1,67,67,90,65,110,111,116,104,101,114,67,104,114,0,67,80,83,133,26,96,119,5,132,0,0,0,0,192,118,5,132,0,0,0,0,0,120,5,132,0,0,0,0,16,120,5,132,0,0,0,0,48,120,5,132,0,0,0,0,64,120,5,132,0,0,0,0,80,120,5,132,0,0,0,0,240,120,5,132,0,0,0,0,32,121,5,132,0,0,0,0,64,121,5,132,0,0,0,0,80,121,5,132,0,0,0,0,32,123,5,132,0,0,0,0,48,123,5,132,0,0,0,0,48,124,5,132,0,0,0,0,16,123,5,132,0,0,0,0,0,132,6,132,0,0,0,0,32,132,6,132,0,0,0,0,48,134,6,132,0,0,0,0,144,134,6,132,0,0,0,0,160,134,6,132,0,0,0,0,176,134,6,132,0,0,0,0,32,135,6,132,0,0,0,0,160,135,6,132,0,0,0,0,96,137,6,132,0,0,0,0,112,137,6,132,0,0,0,0,48,138,6,132,0,0,0,0,128,138,6,132,0,0,0,0,240,138,6,132,0,0,0,0,0,139,6,132,0,0,0,0,240,139,6,132,0,0,0,0,16,140,6,132,0,0,0,0,112,141,6,132,0,0,0,0,192,141,6,132,0,0,0,0,64,142,6,132,0,0,0,0,48,141,6,132,0,0,0,0,224,154,6,132,0,0,0,0,32,155,6,132,0,0,0,0,96,142,6,132,0,0,0,0,112,142,6,132,0,0,0,0,128,142,6,132,0,0,0,0,160,143,6,132,0,0,0,0,192,143,6,132,0,0,0,0,128,143,6,132,0,0,0,0,192,144,6,132,0,0,0,0,0,145,6,132,0,0,0,0,64,137,6,132,0,0,0,0,16,147,6,132,0,0,0,0,32,147,6,132,0,0,0,0,80,147,6,132,0,0,0,0,112,124,5,132,0,0,0,0,144,124,5,132,0,0,0,0,192,124,5,132,0,0,0,0,32,125,5,132,0,0,0,0,80,125,5,132,0,0,0,0,48,121,5,132,0,0,0,0,96,157,6,132,0,0,0,0,176,157,6,132,0,0,0,0,0,158,6,132,0,0,0,0,48,158,6,132,0,0,0,0,224,158,6,132,0,0,0,0,128,125,5,132,0,0,0,0,176,125,5,132,0,0,0,0,64,126,5,132,0,0,0,0,16,121,5,132,0,0,0,0,240,161,6,132,0,0,0,0,80,162,6,132,0,0,0,0,80,126,5,132,0,0,0,0,144,126,5,132,0,0,0,0,176,126,5,132,0,0,0,0,192,126,5,132,0,0,0,0,224,126,5,132,0,0,0,0,32,127,5,132,0,0,0,0,96,127,5,132,0,0,0,0,16,127,5,132,0,0,0,0,16,128,6,132,0,0,0,0,96,128,6,132,0,0,0,0,64,128,6,132,0,0,0,0,64,130,6,132,0,0,0,0,96,114,5,132,0,0,0,0,80,164,6,132,0,0,0,0,96,164,6,132,0,0,0,0,176,164,6,132,0,0,0,0,192,164,6,132,0,0,0,0,80,165,6,132,0,0,0,0,96,165,6,132,0,0,0,0,208,165,6,132,0,0,0,0,240,165,6,132,0,0,0,0,224,166,6,132,0,0,0,0,112,167,6,132,0,0,0,0,208,168,6,132,0,0,0,0,240,168,6,132,0,0,0,0,160,169,6,132,0,0,0,0,32,170,6,132,0,0,0,0,224,170,6,132,0,0,0,0,0,171,6,132,0,0,0,0,192,175,6,132,0,0,0,0]
	# TODO: work with BAM from memory instead of from file?

	fn = "BamReader/data/singleread.bam"
	bamFile = BamFile(fn)

	# verify header
	@test filename(bamFile) == fn
	@test header(bamFile) == "@HD\tVN1.3\n@SQ\tSN:TheChrName\tLN:1000000\n"
	@test sequences(bamFile) == [("TheChrName",1000000)]

	for (i,rd) in enumerate(bamFile)
		@test i == 1 # there is only one read in this file

		for r in [rd,copy(rd)] # test on both read and copy of read
			verifymyread(r)
		end
	end

end


text:
function checkpointgeneric(etype::String, config::Dict, jobid::Integer, checkpointcallback!::Function)
    if jobid % 50 == 0
        @info ("Main Thread: Processing job $(jobid)/$(config["nruns"]), jobID $(config["EXPERIMENT_ID"])")
        lock(config["checkpointlock"])
        try
            checkpoint = Dict()
            try
                checkpoint = loadcheckpoint(etype)  # Fails if file not found -> Create new checkpoint dict
            catch
            end
            checkpointcallback!(checkpoint, config)
            savecheckpoint(etype, checkpoint)
        catch e
            # throw(e)
        finally
            unlock(config["checkpointlock"])
        end
    end
end

function parallelrun(etype::String, experiment::Function, config::Dict, checkpointcallback!::Function, resultcallback!::Function, finalcallback!::Function)
    Threads.@threads for jobid in config["nprocessed"]+1:config["nruns"]
        # Increment nprocessed (threadsafe)
        lock(config["paramlock"])
        try
            @assert config["nprocessed"] < config["nruns"]
            config["nprocessed"] += 1
        catch
            continue
        finally
            unlock(config["paramlock"])
        end
        # Main experiment
        result = nothing  # Define in outer scope to avoid undefined variable in case of experiment failure
        try
            result = experiment(config, jobid)  # GaussianProcesses.optimize! spams exceptions
        catch e
            display(e)
            # throw(e)
        end
        lock(config["resultlock"])
        # Writing the results
        try
            resultcallback!(result, config)
        catch e
            display(e)
            # throw(e)
        finally
            unlock(config["resultlock"])
        end
        checkpointgeneric(etype, config, jobid, checkpointcallback!)  # Threadsafe
    end  # End of threaded program
    results = Dict()
    try
        results = loadcheckpoint(etype*"_final")  # Fails if file not found -> Create new results dict
    catch
    end
    finalcallback!(results, config)
    savecheckpoint(etype*"_final", results)
    @info ("Main Thread: Parallel run $(config["EXPERIMENT_ID"]) finished successfully.")
end

function parallelsim(experiment::Function, config::Dict; idmod::String = "")
    etype = "noisy" * idmod  # Mean dynamics, sin need different ID
    function checkpointcallback!(checkpoint::Dict, config::Dict)
        checkpoint[config["EXPERIMENT_ID"]] = Dict("nprocessed" => config["nprocessed"], "kstep_mse" => config["kstep_mse"])  # Append to results if any
    end

    function resultcallback!(result, config::Dict)
        result[1] !== nothing ? kstep_mse = simulationerror(result[2], result[1]) : kstep_mse = nothing
        (length(result) == 3 && result[3] !== nothing) ? projectionerror = result[3] : projectionerror = 0
        push!(config["kstep_mse"], kstep_mse)
        push!(config["projectionerror"], projectionerror)
    end

    function finalcallback!(results, config::Dict)
        results[config["EXPERIMENT_ID"]] = Dict("nprocessed" => config["nprocessed"], "kstep_mse" => config["kstep_mse"], "projectionerror" => config["projectionerror"])
    end
    parallelrun(etype, experiment, config, checkpointcallback!, resultcallback!, finalcallback!)
end

function parallelsearch(experiment::Function, config::Dict)
    etype = "params"

    function checkpointcallback!(checkpoint::Dict, config::Dict)
        checkpoint[config["EXPERIMENT_ID"]] = Dict("nprocessed" => config["nprocessed"], "params"=> config["params"], "kstep_mse" => config["kstep_mse"])  # Append to results if any
    end

    function resultcallback!(result, config::Dict)
        result[1] !== nothing ? kstep_mse = simulationerror(result[2], result[1]) : kstep_mse = Inf
        push!(config["kstep_mse"], kstep_mse)
        push!(config["params"], result[3])
    end

    function finalcallback!(results, config::Dict)
        results[config["EXPERIMENT_ID"]] = Dict("nprocessed" => config["nprocessed"], "params" => config["params"], "kstep_mse" => config["kstep_mse"])
    end

    parallelrun(etype, experiment, config, checkpointcallback!, resultcallback!, finalcallback!)
end


text:
import RowEchelon

#=
Homework 5
Section 2.1
Problem 2
=#

# Initialize all Matrices
A = [2 0 -1 ; 4 -5 2]
B = [7 -5 1 ; 1 -4 -3]
C = [1 2 ;-2 1]
D = [3 5 ; -1 4]
E = [-5 ; 3]

println("A + 3B: ")
println(A+(3*B))

# After this point, products may not exist
# Use try/catch to find nonexistent products.
println("2C - 3E:")
try
	println((2*C) - (3*E))
catch
	println("2C - 3E cannot be computed.")
end

println("DB:")
try
	println(D*B)
catch
	println("DB cannot be computed.")
end

println("EC:")
try
	println(E*C)
catch
	println("EC cannot be computed.")
end


text:
const MAX_SHOW_COMPONENTS = 10
const MAX_SHOW_FORECASTS = 10
const MAX_SHOW_FORECAST_INITIAL_TIMES = 1

function Base.summary(components::Components)
    return "$(typeof(components)): $(get_num_components(components))"
end

function Base.show(io::IO, components::Components)
    i = 1
    for component in iterate_components(components)
        if i <= MAX_SHOW_COMPONENTS
            show(io, component)
            println(io)
        end
        i += 1
    end

    if i > MAX_SHOW_COMPONENTS
        num = i - MAX_SHOW_COMPONENTS
        println(io, "\n***Omitted $num components***\n")
    end
end

function Base.show(io::IO, ::MIME"text/plain", components::Components)
    df = create_components_df(components)
    println(io, "Components")
    println(io, "==========")
    println(io, "Num components: $(get_num_components(components))\n")
    show(io, df)
end

function Base.show(io::IO, ::MIME"text/html", components::Components)
    df = create_components_df(components)
    println(io, "
Components
")
    println(io, "
Num components: $(get_num_components(components))
")
    withenv("LINES" => 100, "COLUMNS" => 200) do
        show(io, MIME"text/html"(), df)
    end
end

function Base.summary(forecast::Forecast)
    return "$(typeof(forecast)) forecast"
end

function Base.summary(forecast::ForecastInternal)
    return "$(typeof(forecast)) forecast"
end

function Base.show(io::IO, data::SystemData)
    show(io, data.components)
    println(io, "\n")
    show(io, data.forecast_metadata)
end

function Base.show(io::IO, ::MIME"text/plain", data::SystemData)
    show(io, MIME"text/plain"(), data.components)
    println(io, "\n")

    println(io, "Forecasts")
    println(io, "=========")
    println(io, "Resolution: $(Dates.Minute(get_forecasts_resolution(data)))")
    println(io, "Horizon: $(get_forecasts_horizon(data))")
    initial_times = [string(x) for x in get_forecast_initial_times(data)]
    println(io, "Initial Times: $(join(initial_times, ", "))")
    println(io, "Interval: $(get_forecasts_interval(data))")
    component_count, forecast_count = get_forecast_counts(data)
    println(io, "Components with Forecasts: $component_count")
    println(io, "Total Forecasts: $forecast_count")
end

function Base.show(io::IO, ::MIME"text/html", data::SystemData)
    show(io, MIME"text/html"(), data.components)
    println(io, "\n")

    println(io, "
Forecasts
")
    println(io, "
Resolution: $(Dates.Minute(get_forecasts_resolution(data)))
")
    println(io, "
Horizon: $(get_forecasts_horizon(data))
")
    initial_times = [string(x) for x in get_forecast_initial_times(data)]
    println(io, "
Initial Times: $(join(initial_times, ", "))
")
    println(io, "
Interval: $(get_forecasts_interval(data))
")
    component_count, forecast_count = get_forecast_counts(data)
    println(io, "
Components with Forecasts: $component_count
")
    println(io, "
Total Forecasts: $forecast_count
")
end

function Base.summary(ist::InfrastructureSystemsType)
    # All InfrastructureSystemsType subtypes are supposed to implement get_name.
    # Some don't.  They need to override this function.
    return "$(get_name(ist)) ($(typeof(ist)))"
end

function Base.show(io::IO, ::MIME"text/plain", ist::InfrastructureSystemsType)
    print(io, summary(ist), ":")
    for name in fieldnames(typeof(ist))
        name == :internal && continue
        val = getfield(ist, name)
        # Not allowed to print `nothing`
        if isnothing(val)
            val = "nothing"
        end
        print(io, "\n   ", name, ": ", val)
    end
end

function create_components_df(components::Components)
    counts = Dict{String, Int}()
    rows = []

    for (subtype, values) in components.data
        type_str = strip_module_name(string(subtype))
        counts[type_str] = length(values)
        parents = [strip_module_name(string(x)) for x in supertypes(subtype)]
        row = (ConcreteType=type_str,
               SuperTypes=join(parents, " <: "),
               Count=length(values))
        push!(rows, row)
    end

    sort!(rows, by = x -> x.ConcreteType)

    return DataFrames.DataFrame(rows)
end

function Base.show(io::IO, ::MIME"text/plain", period::Union{Dates.TimePeriod, Dates.DatePeriod})
    total = convert_compound_period(period)
    println(io, "$total")
end

function Base.show(io::IO, ::MIME"text/html", period::Union{Dates.TimePeriod, Dates.DatePeriod})
    total = convert_compound_period(period)
    println(io, "
$total
")
end

## This function takes in a time period or date period and returns a compound period
function convert_compound_period(period::Union{Dates.TimePeriod, Dates.DatePeriod})
    period = time_period_conversion(period)

    milli_weeks = period - (period % Dates.Millisecond(604800000))
    weeks = convert(Dates.Week, milli_weeks)
    period -= milli_weeks

    milli_days = period - (period % Dates.Millisecond(86400000))
    days = convert(Dates.Day, milli_days)
    period -= milli_days

    milli_hours = period - (period % Dates.Millisecond(3600000))
    hours = convert(Dates.Hour, milli_hours)
    period -= milli_hours

    milli_minutes = period - (period % Dates.Millisecond(60000))
    minutes = convert(Dates.Minute, milli_minutes)
    period -= milli_minutes

    seconds = period - (period % Dates.Millisecond(1000)) # finding the seconds
    seconds = convert(Dates.Second, seconds)
    remainder = period % Dates.Millisecond(1000) #finding the remainding milliseconds
    total = weeks + days + hours + minutes + seconds + remainder
    return total
end

function create_forecasts_df(forecasts::Forecasts)
    initial_times = _get_forecast_initial_times(forecasts.data)
    dfs = Vector{DataFrames.DataFrame}()

    for (i, initial_time) in enumerate(initial_times)
        if i > MAX_SHOW_FORECAST_INITIAL_TIMES
            break
        end
        counts = Dict{String, Int}()
        rows = []

        for (key, values) in forecasts.data
            if key.initial_time != initial_time
                continue
            end

            type_str = strip_module_name(string(key.forecast_type))
            counts[type_str] = length(values)
            parents = [strip_module_name(string(x)) for x in supertypes(key.forecast_type)]
            row = (ConcreteType=type_str,
                   SuperTypes=join(parents, " <: "),
                   Count=length(values))
            push!(rows, row)
        end

        sort!(rows, by = x -> x.ConcreteType)

        df = DataFrames.DataFrame(rows)
        push!(dfs, df)
    end

    return initial_times, dfs
end


text:
# common functions and constants

using DiffModels, Optim

# model base class
abstract type ModelBase end

# constants
const pmin = 1e-100
const sqrttwo = √(2)
const sqrttwoπ = √(2 * π)
const logtwoπ = log(2 * π)
const minwnorm = 1e-60  # smallest size of w (to avoid zero drift)

# returns if currently in interactive environment
inrepl() = isinteractive() || isdefined(Main, :IJulia) && Main.IJulia.inited

# robust mean, cutting top and bottom 10% (p) of values before performing mena
function robustmean(x::AbstractVector, p::Real=0.1)
    i = floor(Integer, p * length(x))
    mean(sort(x)[i:end-i])
end

# cumulative Gaussians utility functions
invΦ(p) = sqrttwo * erfinv(2p - 1)
Φ(x) = (1 + erf(x / sqrttwo)) / 2

# computes N(x|0,1) / Φ(x), more stable for asymptotes at small x
NoverΦ(x) = (x > -6 ? sqrt(2 / π) * exp(-0.5abs2(x)) / (1 + erf(x / sqrttwo)) :
                      -x / (1 - 1/abs2(x) + 3/x^4))

# returns moving average of window size w
function movingavg(x::AbstractVector, w)
    n = length(x)
    if w > n
        return [mean(x)]
    end
    y = Array(typeof(mean(x[1])), n - w + 1)
    # this could be made faster by taking a running moving average, at the cost
    # of numerical stability
    for i = 1:length(y)
        y[i] = mean(slice(x, i:(i+w-1)))
    end
    y
end

# function to find probit fit to psychometric curve
# takes xn vector and associated pn=p(y=1|xn) vector and returns
# μ and σ that maximize likelihood of this data
function psychprobitfit(xn::AbstractVector, pn::AbstractVector)
    n = length(xn)
    @assert n == length(pn)
    # (neg) log-likelihood function, x = [μ, σ]
    function f(x::AbstractVector)
        llh = 0.0
        for i = 1:n
            xi = (xn[i] - x[1]) / x[2]
            llh += pn[i] * log(max(pmin,Φ(xi))) + (1 - pn[i]) * log(max(pmin,Φ(-xi)))
        end
        return -llh
    end
    # fit function using Optim.jl
    res = optimize(DifferentiableFunction(f),
                   [-Inf, 0.0], [Inf, Inf], [0.0, 3.0], Fminbox(LBFGS()))
    return Optim.minimizer(res)
end

# returns orthogonal matrix U with U[1,:] = u, and all other rows normalized
function orth_matrix(u::Vector{T}) where T <: Real
    n = length(u)
    unorm2 = dot(u, u)
    # first row
    U = Matrix{T}(undef, n, n)
    U[1,:] = u
    # following the modified Gram Schmidt procedure, assuming that
    # vk has single 1 at location i[k-1], where i[k] is the index of the
    # kth smallest item in u.
    i = sortperm(abs.(u))
    for k = 2:n
        # perform orthogonalisation with respect to unnormalized u1 separately
        ik = i[k-1]
        uk = -(u[ik] / unorm2) * u
        uk[ik] += 1
        for j = 2:(k-1)
            uj = view(U, j, :)      # normalized for all j >= 2
            uk -= dot(uk, uj) * uj
        end
        U[k,:] = uk / norm(uk)
    end
    U
end

# draws a sample from the multivariate Gaussian N(x | μ, Σ)
# while imposing constraint u^T x = z
function randConstMvn(μ, Σ, u, z)
    n = length(μ)
    # construct orthogonal matrix U such that U[1,:] = u
    U = orth_matrix(u)
    # use y = U x, such that y ~ N(U μ, U Σ U^T), and sample y[2:end] | y[1] = z
    μy = U * μ
    Σy = U * Σ * U'
    y = Vector{Float64}(undef, n)
    y[1] = z
    # find Λdd^-1, where Σy^-1 = [Λaa Λbb; Λcc Λdd], using matrix inv. lemma
    invΛdd = Σy[2:end, 2:end] .- (Σy[2:end, 1] * Σy[1, 2:end]') / Σy[1, 1]
    μd = μy[2:end] .+ Σy[2:end, 1] / Σy[1, 1] * (z .- μy[1])
    y[2:end] = rand(MvNormal(μd, invΛdd))
    # map back into x
    return U \ y
end

#NOTE TO SELF: z0 here needs to be re-scaled by σ2 / σμ2 already!

# simulates diffusion and returns (response, t, x)
# dx ~ N( c dt, I dt )     (c is a vector)
# and diffusion model bounds {-θ(t)-z0, θ(t)-z0} on u^T x, with both u and z0
# given as parameters. (θ(t)-z0)/||u|| is specified through b.
# The returned response is either 1 (upper) or 0 (lower bound)
function _simdiffusion(c, u, z0, b, δt, maxt)
    # sample (z, t) from 1D diffusion model, bounds at { θ-z0, -θ-z0 }
    μz = dot(u, c)
    σz = norm(u)
    t, resp = 0.0, true
    try
        t, resp = rand(sampler(ConstDrift(μz/σz, δt), b))
        t = min(t, maxt) # in case t > tmax
    catch y
        println("c = $(c), u = $(u)")
        println("μz = $(μz), σz = $(σz), θ = $(θ), z0 = $(z0)")
        println("μz/σz = $(μz/σz), (θ-z0)/σz = $((θ-z0)/σz)")
        flush(STDOUT)
        throw(y)
    end
    # get distance travelled z(t) = ±θ(t) - z0 at decision time t
    n = floor(Int, t / δt) + 1
    z = σz * (resp ?                       # as b[n] = (θ(t) - z0) / σz
              DiffModels.getubound(b, n) : # we get θ(t) - z0 by
              DiffModels.getlbound(b, n))  # b[n] * σz
    # find x(t) for given z(t)
    if length(c) == 1
        # single dimension, x is just scaled z = u x
        x = z ./ u
    else
        # draw x ~ N(μx(t), Σx(t)), subject to constraing u^T x(t) = z(t)
        μx = c * t
        Σx = Matrix(t * I, length(c), length(c))
        x = randConstMvn(μx, Σx, u, z)
    end
    # returns (response, t, x)
    return resp ? 1 : 0, t, x
end
# makes sure that ||u|| ≥ minwnorm * √(length(u)) and returns u, ||u||
function _lowerboundu(u)
    Dc = length(u)
    unorm2 = dot(u, u)
    if unorm2 < Dc * minwnorm^2
        return minwnorm * ones(Dc), minwnorm * √(Dc)
    else
        return u, √(unorm2)
    end
end
# simulate diffusion with bound vector θt in steps of δt, with θ(0) = θt[1]
function simdiffusionvec(c, u, z0::Real, θt::Vector, δt::Real, maxt::Real)
    z0 < θt[1]  || return 1, 0.0, zero(c)  # z0 above upper boundary
    z0 > -θt[1] || return 0, 0.0, zero(c)  # z0 below lower boundary
    u, σz = _lowerboundu(u)
    b = AsymBounds{VarBound, VarBound}(
        VarBound([(θi - z0)/σz for θi in θt], δt),
        VarBound([(θi + z0)/σz for θi in θt], δt))
    return _simdiffusion(c, u, z0, b, δt, maxt)
end
# simulate diffusion with constant bound at θ
function simdiffusionconst(c, u, z0::Real, θ::Real, maxt::Real)
    z0 < θ  || return 1, 0.0, zero(c)  # z0 above upper boundary
    z0 > -θ || return 0, 0.0, zero(c)  # z0 below lower boundary
    δt = 1e-5  # required for sampler, but not used for const. boundaries
    u, σz = _lowerboundu(u)
    b = ConstAsymBounds((θ-z0)/σz, -(θ+z0)/σz, δt)
    return _simdiffusion(c, u, z0, b, δt, maxt)
end
# simulate diffusion with linearly changing boundary θ + t * θslo
function simdiffusionlin(c, u, z0, θ::Real, θslo::Real, δt::Real, maxt::Real)
    z0 < θ  || return 1, 0.0, zero(c)  # z0 above upper boundary
    z0 > -θ || return 0, 0.0, zero(c)  # z0 below lower boundary
    u, σz = _lowerboundu(u)
    b = AsymBounds{LinearBound, LinearBound}(
        LinearBound((θ-z0)/σz, θslo/σz, δt, maxt),
        LinearBound((θ+z0)/σz, θslo/σz, δt, maxt))
    return _simdiffusion(c, u, z0, b, δt, maxt)
end

# simulates leaky competing accumulator model and returns (response, t, x)
# dxi = (ci - τinvi xi - winhi xj) dt + dWi, where i ∈ {1, 2}, j = 3-i,
# and bounds {θ1, θ2} on {x1, x2}. The returned response is either 0
# (x1 "wins") or 1 (x2 "wins").
function simlca(c, τinvi, winhi, θ1, θ2, δt, maxt)
    # pre-compute constants
    θ1n, θ2n, maxn = length(θ1), length(θ2), ceil(Int, maxt/δt)
    cδt = c .* δt
    τinviδt = 1 .- τinvi .* δt
    winhiδt = winhi .* δt
    sqrtδt = √(δt)

    # iterate in steps of size δt
    x = zeros(2)
    n = 0
    while n ≤ maxn
        x .= max.(0.0,
            τinviδt .* x .+ cδt .- winhiδt .* (@view x[end:-1:1]) .+
            sqrtδt .* randn(2))
        n += 1
        if x[1] ≥ θ1[min(θ1n, n)]
            # choose randomly if both races crossed the bound simultaneously
            return (x[2] ≥ θ2[min(θ2n, n)] && rand() < 0.5 ? 1 : 0), (n-1)*δt, x
        elseif x[2] ≥ θ2[min(θ2n, n)]
            return 1, (n-1)*δt, x
        end
    end
    return (x[2] > x[1] ? 1 : 0), (n-1)*δt, x
end


text:
function strengthened_comb_inequalities!(
    demand::Vector{Int},
    capacity::Int,
    edge_tail::Vector{Int},
    edge_head::Vector{Int},
    edge_x::Vector{Float64};
    max_n_cuts = 1000::Cint,
    Q_min = -1,
    K = -1
)
    n_customers, _demand, n_edges, _edge_tail, _edge_head, _edge_x =
        input_conversion(demand, edge_tail, edge_head, edge_x)

    new_cmp = init_CMP()

    _max_violation = Ref{Cdouble}()

    if Q_min == -1
        if K == -1
            K = n_customers
        end
        Q_min = sum(_demand) - (K - 1) * capacity
    end

    # @show n_customers, _demand, capacity, Q_min
    # @show n_edges
    # @show _edge_tail
    # @show _edge_head
    # @show _edge_x

    ccall(
        (:COMBSEP_SeparateCombs, LIB_CVRPSEP),
        Cvoid,
        (
            Cint, Ptr{Cint}, Cint, Cint,
            Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
            Cint, Ref{Cdouble}, CnstrMgrPointer
        ),
            Cint(n_customers), Cint.(_demand), Cint(capacity), Cint(Q_min),
            Cint(n_edges), Cint.(_edge_tail), Cint.(_edge_head), Cdouble.(_edge_x),
            Cint(max_n_cuts), _max_violation, new_cmp
    )

    new_cuts = convert_CMP(new_cmp)

    n_new_cuts = length(new_cuts)

    handles = Vector{Vector{Int}}(undef, n_new_cuts)
    teeth = Vector{Vector{Vector{Int}}}(undef, n_new_cuts)
    rhs = Vector{Float64}(undef, n_new_cuts)

    for i in 1:n_new_cuts
        cr = new_cuts[i]
        n_teeth = cr.Key

        # Handle
        handles[i] = cr.IntList .+ 1
        depot = n_customers + 2
        replace!(x -> x == depot ? 1 : x, handles[i])

        # Teeth
        teeth[i] = Vector{Int}[]
        for t in 1:n_teeth
            min_idx = cr.ExtList[t]
            if t == n_teeth
                max_idx = length(cr.ExtList)
            else
                max_idx = cr.ExtList[t + 1] - 1
            end
            tooth = cr.ExtList[min_idx:max_idx] .+ 1
            replace!(x -> x == depot ? 1 : x, tooth)
            push!(teeth[i], tooth)
        end

        # RHS
        rhs[i] = cr.RHS
    end

    return handles, teeth, rhs
end


text:
"""
csv_to_dataframe(csv_path::String,file_regexp::String)::DataFrame

Convert CSV files to a DataFrame.

# Arguments
- `csv_path::String`: path of csv files
- `file_regexp::String`: regular expression of csv files to convert DataFrame

# Example

    csv_to_dataframe("../data/csv/monthly/"; file_regexp=r"(.*)")

"""
function csv_to_dataframe(csv_path::String;
                        file_regexp::Regex=r"(.*)",delim::String=",",
                        header::Bool=false)::DataFrame
    # get all list of directory
    lst_dir = readdir(abspath(csv_path),join=true)
    # select file only
    lst_files = lst_dir[isfile.(lst_dir)]
    df = DataFrame()
    for f in lst_files
        m = match(file_regexp,f)
        append!(df,DataFrame(CSV.File(m.match;header=header,delim=delim)))
    end
    return df
end

max_stars_repo_name: mrchaos/fxrobot

max_issues_repo_name: mrchaos/fxrobot

max_forks_repo_name: mrchaos/fxrobot


text:
"plotVF!(fx,fy,x,y) plot vector field"
function plotVF!(p,fx::Function,fy::Function,xbase::Union{Array,UnitRange,StepRangeLen},ybase::Union{Array,UnitRange,StepRangeLen};map_speed = true,log_speed = false)

    dxbase = xbase[2] - xbase[1]
    dybase = ybase[2] - ybase[1]
    dxybase = sqrt(dxbase^2+dybase^2)
    cprod = collect.(collect(Base.product(xbase,ybase))) # not writing for-loop, use Cartesian product of x base and y base, cprod is a permutation set
    cprod = vec(cprod)
    x = map(x->x[1],cprod)
    y = map(x->x[2],cprod)
    dx=map(x->fx(x[1],x[2]),cprod)
    dy=map(x->fy(x[1],x[2]),cprod)
    speed = sqrt.(dx.^2 .+ dy.^2)#[sqrt(dx[i]^2+dy[i]^2) for i in 1:length(dx)]

    if log_speed
        speed = log.(speed)
        dx=log.(dx)
        dy = log.(dy)
    end

    #normalization
    dx ./=maximum(abs.(dx))
    dx .*=dxbase
    dy ./=maximum(abs.(dy))
    dy .*=dybase
    #show(dx);show(dy)
    #how(speed)

    if map_speed
        return quiver!(p,x,y,quiver = (dx,dy))#;markerstrokecolor = qcMap_def(speed))
    end
    return quiver!(p,x,y,quiver = (dx,dy))
end

#fx(x,y)=x;fy(x,y)=y
#plotVF!(plot(),fx,fy,collect(-3:1.0:3),collect(-3:1.0:3))

#quiver(1:10,1:10,quiver=(10:-1:1,10:-1:1),markerstrokecolor = qcMap_def(1:10))

"plotGradient!(p,x) - autonomous"
function plotGradient!(p,fx::Function,xbase;t = 0:1.0:5,map_speed = true,log_speed = false,normalize = true)
    if t==Missing
        t = 1:1.0:length(x)
    end
    dt = t[2]-t[1]
    dx = xbase[2] - xbase[1]

    if isinteger(t[1])
        t = Float64.(t)
    end
    if isinteger(xbase[1])
        xbase = Float64.(xbase)
    end

    dtbase = t[2] - t[1]
    dxbase = xbase[2] - xbase[1]
    dbase = sqrt(dt^2+dt^2)
    cprod = collect.(collect(Base.product(t,xbase))) # not writing for-loop, use Cartesian product of x base and y base, cprod is a permutation set
    cprod = vec(cprod)
    x = map(x->x[2],cprod)
    t = map(x->x[1],cprod)
    dt= repeat([dtbase],length(t))
    dx=map(x->fx(x[2]),cprod)
    speed = sqrt.(dt.^2 .+ dx.^2)#[sqrt(dt[i]^2+dx[i]^2) for i in 1:length(dt)]

    if log_speed
        speed = log.(speed)
        dt=log.(dt)
        dx = log.(dx)
    end

    #normalization
    if normalize
        dx ./=maximum(abs.(dx))
        dx .*=dxbase
        #shrinks t vec also so that gradient is preserved

        dt ./=maximum(abs.(dx))
        dt .*-dtbase
    end
    #show(dt);show(dx)
    #how(speed)

    if map_speed
        return quiver!(p,t,x,quiver = (dt,dx))#;markerstrokecolor = qcMap_def(speed))
    end
    return quiver!(p,t,x,quiver = (dt,dx))
end

"plotNullCline(fx,fy)"
function plotNullCline!(p,fx::Function,fy::Function,xbase,ybase;tol = 1e-3)
    p=contour!(p,xbase,ybase,[fx(xi,yi) for yi in ybase, xi in xbase],levels = [-tol,tol],color = :red,alpha = .3)#,key = "x NullCline est")
    p=contour!(p,xbase,ybase,[fy(xi,yi) for yi in ybase, xi in xbase],levels = [-tol,tol],color =:blue ,alpha = .3)#,key = "y NullCline est")
    return p
end

#contour(1:2,1:2,[.11 .02;.1 0],levels = [-.1,.1])

"plotValongCycle!(p,fx::Function,fy::Function,cycle with r) works for 2d system only"
function plotValongCycle!(p,fx::Function,fy::Function,Cycle::Function;theta = 0.0:.3:2*pi)
    xl = zeros(length(theta))
    yl = similar(xl)
    dxl =similar(yl)
    dyl =similar(dxl)
    for i in 1:length(theta)
        err = ""
        try
            x,y = Cycle(theta[i])
            xl[i] = x
            yl[i] = y
            dxl[i] = fx(x,y)
            dyl[i] = fy(x,y)
        catch err
            println("Singularity?")
        end
    end
    p = quiver!(p,xl,yl,quiver = (dxl,dyl))
    p = plot!(p,xl,yl;alpha = .3,key = false)
    return p
end

function circularCycle(theta;center = (0.0,0.0),r = 1.0)
    [center[1]+ r * cos(theta), center[2] + r * sin(theta) ]
end

function plotValongCycle!(p,fx::Function,fy::Function;r = 1.0,center = (0.0,0.0),theta = .3:.3:2*pi)
    cycle(theta) = circularCycle(theta;center = center,r=r)
    return plotValongCycle!(p,fx,fy,cycle;theta = .3:.3:2*pi)
end


text:
module FreeSurfacePotentialFlowTests

using Gridap
using Gridap.Geometry
using GridapODEs.TransientFETools
using GridapODEs.ODETools
using Test

# Parameters
L = 2*π
H = 1.0
n = 8
order = 2
g = 9.81
ξ = 0.1
λ = L/2
k = 2*π/L
h = L/n
ω = √(g*k*tanh(k*H))
t₀ = 0.0
tf = 2*π
Δt = h/(2*λ*ω)
θ = 0.5

# Exact solution
ϕₑ(x,t) = ω/k * ξ * (cosh(k*(x[2]))) / sinh(k*H) * sin(k*x[1] - ω*t)
ηₑ(x,t) = ξ * cos(k*x[1] - ω*t)
ϕₑ(t::Real) = x -> ϕₑ(x,t)
ηₑ(t::Real) = x -> ηₑ(x,t)

# Domain
domain = (0,L,0,H)
partition = (n,n)
model = CartesianDiscreteModel(domain,partition;isperiodic=(true,false))

# Boundaries
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"bottom",[1,2,5])
add_tag_from_tags!(labels,"free_surface",[3,4,6])

# Triangulation
Ω = Interior(model)
Γ = Boundary(model,tags="free_surface")
dΩ = Measure(Ω,2*order)
dΓ = Measure(Γ,2*order)

# FE spaces
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(Ω,reffe,conformity=:H1)
V_Γ = TestFESpace(Γ,reffe,conformity=:H1)
U = TransientTrialFESpace(V)
U_Γ = TransientTrialFESpace(V_Γ)
X = TransientMultiFieldFESpace([U,U_Γ])
Y = MultiFieldFESpace([V,V_Γ])

# Weak form
α = 2/Δt

# Optimal transient FE Operator:
m((ϕt,ηt),(w,v)) = ∫( 0.5*(α/g*(w*ϕt) + v*ϕt) - (w*ηt) )dΓ
a((ϕ,η),(w,v)) = ∫( ∇(ϕ)⋅∇(w) )dΩ + ∫( 0.5*(α*(w*η) + g*v*η) )dΓ
b((w,v)) = ∫( 0.0*w )dΓ
op_const = TransientConstantFEOperator(m,a,b,X,Y)

# TransientFEOperator exploiting automatic differentiation (testing purposes)
res(t,x,y) = m(∂t(x),y) + a(x,y) - b(y)
jac(t,x,dx,y) = a(dx,y)
jac_t(t,x,dxt,y) = m(dxt,y)
op_trans = TransientFEOperator(res,jac,jac_t,X,Y)
op_ad = TransientFEOperator(res,X,Y)

# TransientFEOperator exploiting time derivative of separate fields (TransientMultiFieldCellField)
res2(t,(ϕ,η),y) = m((∂t(ϕ),∂t(η)),y) + a((ϕ,η),y) - b(y)
op_multifield = TransientFEOperator(res2,jac,jac_t,X,Y)

# Solver
ls = LUSolver()
ode_solver = ThetaMethod(ls,Δt,θ)

# Initial solution
x₀ = interpolate_everywhere([ϕₑ(0.0),ηₑ(0.0)],X(0.0))

function test_operator(op)
  # Solution
  sol_t = solve(ode_solver,op,x₀,t₀,tf)

  # Post-process
  l2_Ω(w) = √(∑(∫(w*w)dΩ))
  l2_Γ(v) = √(∑(∫(v*v)dΓ))
  E_kin(w) = 0.5*∑( ∫(∇(w)⋅∇(w))dΩ )
  E_pot(v) = g*0.5*∑( ∫(v*v)dΓ )
  Eₑ = 0.5*g*ξ^2*L

  tol = 1.0e-2
  for ((ϕn,ηn),tn) in sol_t
    E = E_kin(ϕn) + E_pot(ηn)
    error_ϕ = l2_Ω(ϕn-ϕₑ(tn))
    error_η = l2_Γ(ηn-ηₑ(tn))
    @test abs(E/Eₑ-1.0) <= tol
    @test error_ϕ <= tol
    @test error_η <= tol
  end
end

test_operator(op_const)
test_operator(op_trans)
test_operator(op_multifield)
# test_operator(op_ad) # Not working yet

end


text:
__precompile__()

module NetalignMeasures

using MatrixNetworks
using Distances
using StatsBase
using Distributions
import Distributions: dim

export NetalignMeasure, NetalignScore, measure, score, dim

"""
    abstract type NetalignMeasure end
    abstract type NetalignScore end

Sub-types of `NetalignMeasure` are network alignment measures.
A type instance will contain information such that given the
type instance and an alignment, we are able to calculate
the alignment score.

Sub-types of `NetalignScore` contain network alignment score information.
The type instance will contain information of a particular network alignment measure and
an alignment. The `score` field in the type instance will contain the network alignment score.
`score` measures network similarity and will be between `0` and `1`.

### Two networks

Namely, let `S3Measure <: NetalignMeasure` be the S3 network alignment measure.
It will be instantiated as
    meas = S3Measure(G1,G2)
where `G1` and `G2` are the adjacency matrices of the two relevant
networks such that `n1 <= n2`, where `n1` and `n2` are the number of nodes
in `G1` and `G2`, respectively.

Then, given an alignment (described using a vector of `Int`s, `f`, such that when
node `i` in graph `G1` is mapped to node `j` in graph `G2`, `f[i] = j`),
we can calculate the alignment score using the `measure` and `score` functions.

Given a type `NAM <: NetalignMeasure`, the interface functions such a type
under `NetalignMeasure` include:

    measure(meas::NAM, f::Vector{Int}) -> NAS
Returns a struct of type `NAS <: NetalignScore`, which corresponds to type `NAM`,
containing information regarding the alignment score
of `f` with respect the alignment measure `NAM`.

Namely, given an instance of `NAM`, for example, `meas = NAM(G1,G2)`, then
    x = measure(meas,f)
`x.score` will contain the alignment score.
`x` might also contain other information such as the numerator
and denominator of the alignment score, or the number of conserved
and non-conserved edges, etc.

    score(meas::NAM, f::Vector{Int}) -> Float64
Returns the alignment score of `f` with respect the alignment measure `NAM`.

### Multiple networks

Not all `NetalignMeasure`s work for more than two networks. One that does it the
LCCS (largest common conserved sub-graph) measure.
This is an example with 3 networks but it is easily extended to `k` networks.
Namely, let `LCCSMeasure <: NetalignMeasure` be the LCCS network alignment measure.
It will be instantiated as
    meas = LCCSMeasure(G1,G2,G3)
where `G1,G2,G3` are the adjacency matrices of the four relevant
networks such that `n1 <= n2 <= n3`, where `n1,n2,n3` are the number of nodes
in `G1,G2,G3`, respectively.

Then, given an alignment we can calculate the alignment score using the `measure` and `score` functions.
An alignment between `k` networks, `G[1],...,G[k]`, is described using a tuple, `fs`,
of `k-1` vectors of `Int`s.
For a given vector of `Int`s, `fs[l]`, fs[l] is 1-1 mapping from graph `G[l]` to graph `G[k]`.
And similar to the two network case, each `fs[l]` is such that when
node `i` in graph `G[l]` is mapped to node `j` in graph `G[k]`, `fs[l][i] = j`).

Given a type `NAM <: NetalignMeasure`, the interface functions such a type
under `NetalignMeasure` include:

    measure(meas::NAM, fs::Tuple{Vector{Int}}) -> NAS
    measure(meas::NAM, fs::Vector{Int}...) -> NAS
Returns a struct of type `NAS <: NetalignScore`, which corresponds to type `NAM`,
containing information regarding the alignment score
of `fs` with respect the alignment measure `NAM`.

    score(meas::NAM, fs::Tuple{Vector{Int}}) -> Float64
    score(meas::NAM, f::Vector{Int}...) -> Float64
Returns the alignment score of `f` with respect the alignment measure `NAM`.

    score(x::NAS) -> Float64
Returns `x.score`.

    dim(meas::NAM) -> Int
Returns the number of nodes in the largest network among the relevant networks.
Note: this is always the last network since the networks
_must_ be ordered w.r.t. node size.

    dim : NAM x Int -> Int
Returns the number of nodes in the kth network.

### Other functions

    score(x::NAS) -> Float64
Returns `x.score`.

    dim(meas::NAM) -> Int
Returns the number of nodes in the largest network among the relevant networks.
Note: this is always the last network since the networks
_must_ be ordered w.r.t. node size.

    dim : NAM x Int -> Int
Returns the number of nodes in the kth network.

Note: The input adjacency matrices _must_ be symmetric.
"""
abstract type NetalignMeasure end
function dim end

abstract type NetalignScore end
score(x::NetalignScore) = x.score
score(m::NetalignMeasure,fs::Vector{<:Integer}...) = score(measure(m,fs...))

function measure end
measure(meas::NetalignMeasure, fs::Tuple{Vector{<:Integer}}) =
    error("measure function not implemented for $(typeof(meas))")
measure(meas::NetalignMeasure, fs::Vector{<:Integer}...) = measure(meas,fs)

include("mergesorted.jl")
include("networks.jl")
include("dynamicnetworks.jl")
include("measures.jl")
include("nodemeasures.jl")
include("aln.jl")

end


text:
"""
    load_timeseriesdata(data_path::String; T-#Segments,
years::Array{Int64,1}=# years to be selected for the time series, att::Array{String,1}=# attributes to be loaded)
- Loading all `*.csv` files in the folder or the file `data_path`
The `*.csv` files shall have the following structure and must have the same length:
|Timestamp |[column names...]|
|[iterator]|[values]         |
The first column should be called `Timestamp` if it contains a time iterator
The other columns can specify the single timeseries like specific geolocation.
Each column in `[file name].csv` file will be added to the ClustData.data called `"[file name]-[column name]"`
Loads all attributes if the `att`-Array is empty or only the ones specified in it
"""
function load_timeseries_data(data_path;
                              region::String="none",
                              T::Int64=24,
                              years::Array{Int64,1}=[2016],
                              att::Array{String,1}=Array{String,1}())
  dt = Dict{String,Array}()
  num=0
  K=0
  #Check if data_path is directory or file
  if isdir(data_path)
      for full_data_name in readdir(data_path)
          if split(full_data_name,".")[end]=="csv"
              data_name=split(full_data_name,".")[1]
              if isempty(att) || data_name in att
                  # Add the
                  K=add_timeseries_data!(dt, data_name, data_path; K=K, T=T, years=years)
              end
          end
      end
  elseif isfile(data_path)
      full_data_name=splitdir(data_path)[end]
      data_name=split(full_data_name,".")[1]
      K=add_timeseries_data!(dt, data_name, dirname(data_path); K=K, T=T, years=years)
  else
      throw(@error("The path $data_path is neither recognized as a directory nor as a file"))
  end
  # Store the data
  ts_input_data =  ClustData(FullInputData(region, years, num, dt),K,T)
  return ts_input_data
end #load_timeseries_data

"""
    add_timeseries_data!(dt::Dict{String,Array}, data::DataFrame; K::Int64=0, T::Int64=24, years::Array{Int64,1}=[2016])
selects first the years and second the data_points so that their number is a multiple of T and same with the other timeseries
"""
function add_timeseries_data!(dt::Dict{String,Array},
                            data_name::SubString,
                            data_path::String;
                            K::Int64=0,
                            T::Int64=24,
                            years::Array{Int64,1}=[2016])
    #Load the data
    data_df=CSV.read(joinpath(data_path,data_name*".csv");allowmissing=:none)
    # Add it to the dictionary
    return add_timeseries_data!(dt,data_name, data_df; K=K, T=T, years=years)
end

"""
    add_timeseries_data!(dt::Dict{String,Array}, data::DataFrame; K::Int64=0, T::Int64=24, years::Array{Int64,1}=[2016])
selects first the years and second the data_points so that their number is a multiple of T and same with the other timeseries
"""
function add_timeseries_data!(dt::Dict{String,Array},
                            data_name::SubString,
                            data::DataFrame;
                            K::Int64=0,
                            T::Int64=24,
                            years::Array{Int64,1}=[2016])
    # find the right years to select
    data_selected=data[in.(data[:year],[years]),:]
    for column in eachcol(data_selected, true)
        # check that this column isn't time or year
        if !(column[1] in [:Timestamp,:time,:Time,:Zeit,:year])
            K_calc=Int(floor(length(column[2])/T))
            if K_calc!=K && K!=0
                @error("The time_series $(column[1]) has K=$K_calc != K=$K of the previous")
            else
                K=K_calc
            end
            dt[data_name*"-"*string(column[1])]=Float64.(column[2][1:(Int(T*K))])
        end
    end
    return K
end

"""
        combine_timeseries_weather_data(ts::ClustData,ts_weather::ClustData)
-`ts` is the shorter timeseries with e.g. the demand
-`ts_weather` is the longer timeseries with the weather information
The `ts`-timeseries is repeated to match the number of periods of the longer `ts_weather`-timeseries.
If the number of periods of the `ts_weather` data isn't a multiple of the `ts`-timeseries, the necessary number of the `ts`-timeseries periods 1 to x are attached to the end of the new combined timeseries.
"""
function combine_timeseries_weather_data(ts::ClustData,
                                        ts_weather::ClustData)
    ts.T==ts_weather.T || throw(@error "The number of timesteps per period is not the same: `ts.T=$(ts.T)≢$(ts_weather.T)=ts_weather.T`")
    ts.K<=ts_weather.K || throw(@error "The number of timesteps in the `ts`-timeseries isn't shorter or equal to the ones in the `ts_weather`-timeseries.")
    ts_weather.K%ts.K==0 || @warn "The number of periods of the `ts_weather` data isn't a multiple of the other `ts`-timeseries: periods 1 to $(ts_weather.K%ts.K) are attached to the end of the new combined timeseries."
    ts_data=deepcopy(ts_weather.data)
    ts_mean=deepcopy(ts_weather.mean)
    ts_sdv=deepcopy(ts_weather.sdv)
    for (k,v) in ts.data
        ts_data[k]=repeat(v, 1, ceil(Int,ts_weather.K/ts.K))[:,1:ts_weather.K]
    end
    for (k,v) in ts.mean
        ts_mean[k]=v
    end
    for (k,v) in ts.sdv
        ts_sdv[k]=v
    end

    return ClustData(ts.region, ts_weather.years, ts_weather.K, ts_weather.T, ts_data, ts_weather.weights, ts_mean, ts_sdv, ts_weather.deltas, ts_weather.k_ids)
end