text:
puzzleinput_a = readlines("inputs/day19a.txt")

Symbl = Union{Int,Char}
Rule = Array{Array{Symbl,1},1}

import Base.parse

parse(::Type{Symbl}, str::AbstractString) = '"' ∈ str ? str[2] : parse(Int, str) + 1

function parse(::Type{Rule}, str::AbstractString)
    f(s) = parse.(Symbl, split(strip(s), " "))
    f.(split(str, "|"))
end

function parse_line!(rules, l)
    s = split(l, ":")
    rules[parse(Int, s[1]) + 1] = parse(Rule, s[2])
end

mutable struct State
    symlst::Vector{Symbl}
    ltrpos::Int
end

function topdown(rules, str)
    # Top down parser

    posslst = [State([1], 1)]
    # Set up possibility list

    while true
        if isempty(posslst)
            # If possibility list is empty the parsing is no possible

            return false
        end
        c = pop!(posslst)
        # Get the state c: first element of the possibility list

        if c.ltrpos == length(str) + 1 && !isempty(c.symlst)
            # If the algorithm reached the last char but the possibility
            # list is not empty the parsing is not possible

            return false
        end

        if isempty(c.symlst)
            if c.ltrpos == length(str) + 1
                # If the symbol list in state c is is empty and algorithm
                # reached the last char in string then parsing succeeds

                return true
            end
        else
            fstsym = popfirst!(c.symlst)
            if isa(fstsym, Char)
                if fstsym == str[c.ltrpos]
                    # The first symbol is a lexical symbol

                    newc = State(c.symlst, c.ltrpos + 1)
                    push!(posslst, newc)
                end
            else
                # The first symbol is a nom-terminal

                for r in rules[fstsym]
                    newc = State([r; c.symlst], c.ltrpos)
                    push!(posslst, newc)
                end
            end
        end
    end
end

# First Half

rules1 = fill(Rule(), 132)
for l in puzzleinput_a
    parse_line!(rules1, l)
end

puzzleinput_b = readlines("inputs/day19b.txt")
sum([topdown(rules1, s) for s in puzzleinput_b])

# Second Half

rules2 = fill(Rule(), 132)
for l in puzzleinput_a
    parse_line!(rules2, l)
end

rules2[9] = [[43], [43, 9]]
rules2[12] = [[43, 32], [43, 12, 32]]

sum([topdown(rules2, s) for s in puzzleinput_b])


text:
@testset "data/utils" begin
    X = rand(5, 2)
    wX, T = AGP.wrap_X(X)
    y = rand(5)
    ys = [rand(5) for _ in 1:3]
    indices = rand(1:5, 3)
    data = AGP.wrap_data(wX, y)
    @test AGP.view_x(data, indices) == view(wX, indices)
    l = GaussianLikelihood()
    @test AGP.view_y(l, data, indices) == view(y, indices)
    @test AGP.view_y(l, y, indices) == view(y, indices)
    y_softmax = falses(5, 3)
    data = AGP.wrap_data(wX, y_softmax)
    l = SoftMaxLikelihood(3)
    @test AGP.view_y(l, data, indices) == view(y_softmax, indices, :)

    l = GaussianLikelihood()
    data = AGP.wrap_modata(wX, ys)
    @test AGP.view_y(l, data, indices) == view.(ys, Ref(indices))
    @test AGP.view_y((l,), data, indices) == view.(ys, Ref(indices))
end

ext: jl


text:
load("core.jl")

srand(1)
means = [0.1, 0.1, 0.1, 0.1, 0.9]
n_arms = length(means)
means = means[randperm(n_arms)]
arms = map(mu -> BernoulliArm(mu), means)
println("Best arm is $(ind_max(means))")

total_results = zeros(0, 6)

for epsilon in [0.1:0.1:0.5]
  algo = EpsilonGreedy(epsilon, zeros(Int64, n_arms), zeros(n_arms))
  results = test_algorithm(algo, arms, 5000, 250)
  results = hcat(repmat([epsilon], size(results, 1), 1), results)
  total_results = vcat(total_results, results)
end

csvwrite("algorithms/epsilon_greedy/standard_results.csv", total_results)

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
export DisplayRow, render

struct DisplayRow
    size::Int
    background::Char
    content::Vector{Char}
end

function DisplayRow(len::Int; background=' ')
    isone(textwidth(background)) || (throw(
        "Invalid background character. " *
        "Text width of background character must be 1."
    ))

    return DisplayRow(len, background, fill(background, len))
end

function Base.show(io::IO, row::DisplayRow)
    Base.print(io, "DisplayRow(" *
        "size=$(length(row)), " *
        "background char=Char($(convert(UInt16, row.background)))" *
        ")")
end

Base.length(row::DisplayRow) = row.size

Base.lastindex(row::DisplayRow) = row.size

function get_element_index(row::DisplayRow, display_index::Int)
    (mapreduce(textwidth, +, row.content) < display_index) && (throw(BoundsError))

    accumulate_width = index = pre = post = 0
    for c in row.content
        width = textwidth(c)
        accumulate_width += width
        index += 1
        if accumulate_width >= display_index
            pre = display_index - (accumulate_width-width) - 1
            post = accumulate_width-display_index
            break
        end
    end

    return (i=index, pre=pre, post=post)
end

function Base.deleteat!(row::DisplayRow, display_index::Int)
    index = get_element_index(row, display_index)
    deleteat!(row.content, index.i)
    for i=1:(index.pre+index.post+1)
        insert!(row.content, index.i, row.background)
    end
end

function Base.setindex!(row::DisplayRow, c::Char, display_index::Int)
    start = display_index
    stop = display_index + textwidth(c) - 1

    i1, pre, _ = get_element_index(row, start)
    i2, _, post = get_element_index(row, stop)

    deleteat!(row.content, collect(i1:i2))

    for i=1:post insert!(row.content, i1, row.background) end
    insert!(row.content, i1, c)
    for i=1:pre insert!(row.content, i1, row.background) end
end

function Base.setindex!(row::DisplayRow, str::AbstractString, display_range::UnitRange{Int64})
    start, stop = display_range.start, display_range.stop
    str = padding(str, length(display_range), background=row.background)

    i1, pre, _ = get_element_index(row, start)
    i2, _, post = get_element_index(row, stop)

    deleteat!(row.content, collect(i1:i2))

    for i=1:post insert!(row.content, i1, row.background) end
    for c in reverse(str) insert!(row.content, i1, c) end
    for i=1:pre insert!(row.content, i1, row.background) end
end

function render(stream::IO, row::DisplayRow; pos=(-1, -1))
    T.buffered(stream, row, pos) do buffer, row, pos
        (pos != (-1, -1)) && T.cmove(buffer, pos[1], pos[2])
        T.join(buffer, row.content)
    end
end

render(row::DisplayRow; pos=(-1, -1)) = render(T.out_stream, row, pos=pos)


text:
mutable struct State{T}
    order::Integer

    # Continuous states
    xc::SVector{3,T}
    qc::UnitQuaternion{T}
    vc::SVector{3,T}
    ωc::SVector{3,T}

    # Knot points
    xk::Vector{SVector{3,T}}
    qk::Vector{UnitQuaternion{T}}
    Fk::Vector{SVector{3,T}}
    τk::Vector{SVector{3,T}}

    # Current solution estimate [before step;after step]
    xsol::Vector{SVector{3,T}}
    qsol::Vector{UnitQuaternion{T}}
    vsol::Vector{SVector{3,T}}
    ωsol::Vector{SVector{3,T}}

    # Current equations of motion evaluation
    d::SVector{6,T}

    function State{T}() where T
        xc = zeros(T, 3)
        qc = one(UnitQuaternion{T})
        vc = zeros(T, 3)
        ωc = zeros(T, 3)

        xk = [zeros(T, 3)]
        qk = [one(UnitQuaternion{T})]
        Fk = [zeros(T, 3)]
        τk = [zeros(T, 3)]

        xsol = [zeros(T, 3) for i=1:2]
        qsol = [one(UnitQuaternion{T}) for i=1:2]
        vsol = [zeros(T, 3) for i=1:2]
        ωsol = [zeros(T, 3) for i=1:2]

        d = zeros(T, 6)

        new{T}(0, xc, qc, vc, ωc, xk, qk, Fk, τk, xsol, qsol, vsol, ωsol, d)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, state::State{T}) where {T}
    summary(io, state)
    println(io,"")
    println(io,"xc:   "*string(state.xc))
    println(io,"qc:   "*string(state.qc))
    println(io,"ωc:   "*string(state.ωc))
    println(io,"xk:   "*string(state.xk))
    println(io,"qk:   "*string(state.qk))
    println(io,"Fk:   "*string(state.Fk))
    println(io,"τk:   "*string(state.τk))
    println(io,"xsol: "*string(state.xsol))
    println(io,"qsol: "*string(state.qsol))
    println(io,"vsol: "*string(state.vsol))
    println(io,"ωsol: "*string(state.ωsol))
    println(io,"d:    "*string(state.d))
end

function initknotpoints!(state::State, order)
    state.order = order

    state.xk = [state.xk[1] for i = 1:order]
    state.qk = [state.qk[1] for i = 1:order]
    state.Fk = [state.Fk[1] for i = 1:order]
    state.τk = [state.τk[1] for i = 1:order]

    return
end


text:
# formatting specification

# formatting specification language
#
#  spec  ::= [[fill]align][sign][#][0][width][,][.prec][type]
#  fill  ::= <any character>
#  align ::= '<' | '>'
#  sign  ::= '+' | '-' | ' '
#  width ::= <integer>
#  prec  ::= <integer>
#  type  ::= 'b' | 'c' | 'd' | 'e' | 'E' | 'f' | 'F' | 'g' | 'G' |
#            'n' | 'o' | 'x' | 'X' | 's'
#
# Please refer to http://docs.python.org/2/library/string.html#formatspec
# for more details
#

## FormatSpec type

const _numtypchars = Set(['b', 'd', 'e', 'E', 'f', 'F', 'g', 'G', 'n', 'o', 'x', 'X'])

_tycls(c::Char) =
    (c == 'd' || c == 'n' || c == 'b' || c == 'o' || c == 'x') ? 'i' :
    (c == 'e' || c == 'f' || c == 'g') ? 'f' :
    (c == 'c') ? 'c' :
    (c == 's') ? 's' :
    error("Invalid type char $(c)")

immutable FormatSpec
    cls::Char    # category: 'i' | 'f' | 'c' | 's'
    typ::Char
    fill::Char
    align::Char
    sign::Char
    width::Int
    prec::Int
    ipre::Bool   # whether to prefix 0b, 0o, or 0x
    zpad::Bool   # whether to do zero-padding
    tsep::Bool   # whether to use thousand-separator

    function FormatSpec(typ::Char;
               fill::Char=' ',
               align::Char='\0',
               sign::Char='-',
               width::Int=-1,
               prec::Int=-1,
               ipre::Bool=false,
               zpad::Bool=false,
               tsep::Bool=false)

        if align=='\0'
            align = (typ in _numtypchars) ? '>' : '<'
        end
        cls = _tycls(lowercase(typ))
        if cls == 'f' && prec < 0
            prec = 6
        end
        new(cls, typ, fill, align, sign, width, prec, ipre, zpad, tsep)
    end

    # copy constructor with overrides
    function FormatSpec(spec::FormatSpec;
               fill::Char=spec.fill,
               align::Char=spec.align,
               sign::Char=spec.sign,
               width::Int=spec.width,
               prec::Int=spec.prec,
               ipre::Bool=spec.ipre,
               zpad::Bool=spec.zpad,
               tsep::Bool=spec.tsep)
        new(spec.cls, spec.typ, fill, align, sign, width, prec, ipre, zpad, tsep)
    end
end

function show(io::IO, fs::FormatSpec)
    println(io, "$(typeof(fs))")
    println(io, "  cls   = $(fs.cls)")
    println(io, "  typ   = $(fs.typ)")
    println(io, "  fill  = $(fs.fill)")
    println(io, "  align = $(fs.align)")
    println(io, "  sign  = $(fs.sign)")
    println(io, "  width = $(fs.width)")
    println(io, "  prec  = $(fs.prec)")
    println(io, "  ipre  = $(fs.ipre)")
    println(io, "  zpad  = $(fs.zpad)")
    println(io, "  tsep  = $(fs.tsep)")
end

## parse FormatSpec from a string

const _spec_regex = r"^(.?[<>])?([ +-])?(#)?(\d+)?(,)?(.\d+)?([bcdeEfFgGnosxX])?$"

function FormatSpec(s::AbstractString)
    # default spec
    _fill = ' '
    _align = '\0'
    _sign = '-'
    _width = -1
    _prec = -1
    _ipre = false
    _zpad = false
    _tsep = false
    _typ = 's'

    if !isempty(s)
        m = match(_spec_regex, s)
        if m == nothing
            error("Invalid formatting spec: $(s)")
        end
        (a1, a2, a3, a4, a5, a6, a7) = m.captures

        # a1: [[fill]align]
        if a1 != nothing
            if length(a1) == 1
                _align = a1[1]
            else
                _fill = a1[1]
                _align = a1[2]
            end
        end

        # a2: [sign]
        if a2 != nothing
            _sign = a2[1]
        end

        # a3: [#]
        if a3 != nothing
            _ipre = true
        end

        # a4: [0][width]
        if a4 != nothing
            if a4[1] == '0'
                _zpad = true
                if length(a4) > 1
                    _width = parse(Int,a4[2:end])
                end
            else
                _width = parse(Int,a4)
            end
        end

        # a5: [,]
        if a5 != nothing
            _tsep = true
        end

        # a6 [.prec]
        if a6 != nothing
            _prec = parse(Int,a6[2:end])
        end

        # a7: [type]
        if a7 != nothing
            _typ = a7[1]
        end
    end

    return FormatSpec(_typ;
                      fill=_fill,
                      align=_align,
                      sign=_sign,
                      width=_width,
                      prec=_prec,
                      ipre=_ipre,
                      zpad=_zpad,
                      tsep=_tsep)
end

## formatted printing using a format spec

type _Dec end
type _Oct end
type _Hex end
type _HEX end
type _Bin end

_srepr(x) = repr(x)
_srepr(x::AbstractString) = x
_srepr(x::Char) = string(x)

if isdefined(:Enum)
    _srepr(x::Enum) = string(x)
end

function printfmt(io::IO, fs::FormatSpec, x)
    cls = fs.cls
    ty = fs.typ
    if cls == 'i'
        ix = Integer(x)
        ty == 'd' || ty == 'n' ? _pfmt_i(io, fs, ix, _Dec()) :
        ty == 'x' ? _pfmt_i(io, fs, ix, _Hex()) :
        ty == 'X' ? _pfmt_i(io, fs, ix, _HEX()) :
        ty == 'o' ? _pfmt_i(io, fs, ix, _Oct()) :
        _pfmt_i(io, fs, ix, _Bin())
    elseif cls == 'f'
        fx = float(x)
        if isfinite(fx)
            ty == 'f' || ty == 'F' ? _pfmt_f(io, fs, fx) :
            ty == 'e' || ty == 'E' ? _pfmt_e(io, fs, fx) :
            error("format for type g or G is not supported yet (use f or e instead).")
        else
            _pfmt_specialf(io, fs, fx)
        end
    elseif cls == 's'
        _pfmt_s(io, fs, _srepr(x))
    else # cls == 'c'
        _pfmt_s(io, fs, Char(x))
    end
end

printfmt(fs::FormatSpec, x) = printfmt(STDOUT, fs, x)

pyfmt(fs::FormatSpec, x) = (buf = IOBuffer(); printfmt(buf, fs, x); bytestring(buf))
pyfmt(spec::AbstractString, x) = pyfmt(FormatSpec(spec), x)


text:
function ProfileDifferentiable(f::F, x::PaddedMatrices.AbstractMutableFixedSizeVector{Pm1,T,R,R}, ::Val{P}) where {F,T,Pm1,P,R}

    result = DifferentiableObjects.GradientDiffResult(x)
    chunk = DifferentiableObjects.Chunk(Val{Pm1}())
    gconfig = ForwardDiff.GradientConfig(nothing, x, chunk, ForwardDiff.Tag(nothing, T))
    profile = ProfileWrapper(f,
        MutableFixedSizeVector{P,T}(undef),
        MutableFixedSizeVector{P,eltype(gconfig)}(undef),
        Ref(P), Ref(zero(T))
    )

    gradconfig = GradientConfiguration(profile, result, gconfig)

    OnceDifferentiable(similar(x), similar(x), gradconfig)#, Val{P}())
end

@generated function PrePaddedMatrix(::Val{M}, ::Val{N}, ::Type{T}) where {M,N,T}
    # R, L = calculate_L_from_size((M,N), T)
    R = PaddedMatrices.calc_padding(M, T)
    :(zero(MutableFixedSizeMatrix{$R,$N,$T})) # we want offset
end

function AsymptoticPosteriorFD(f, map_::MAPFD, initial_x::MutableFixedSizeVector{P,T}) where {T,P}

    ValP = Val{P}()
    ValPm1 = DifferentiableObjects.ValM1(Val{P}())

    backtrack = DifferentiableObjects.BackTracking2(Val(2)) # order 2 backtrack vs 3
    # state = DifferentiableObjects.uninitialized_state(nuisance)
    state = DifferentiableObjects.BFGSState(ValPm1)
    od = ProfileDifferentiable(f, DifferentiableObjects.ref_x_old(state), Val{P}())

    Lsmall = PrePaddedMatrix(ValPm1,ValPm1,T)

    AsymptoticPosteriorFD(od, backtrack, state, map_, Lsmall)
end
function AsymptoticPosteriorFD(od::D, method::M, state::S, map_::MAP, Lsmall::SM) where {P,T,Pm1,R,D<:AbstractDifferentiableObject{Pm1,T},M,S,MAP <: MAPFD{P,T,R},SV <: PaddedMatrices.AbstractMutableFixedSizeVector{Pm1,T},R2,SM <: MutableFixedSizeMatrix{R2,Pm1,T}}

    AsymptoticPosteriorFD{P,T,Pm1,R,D,M,S,MAP,R2,SM}(od, method, state, map_, Ref{T}(), Ref{T}(), Ref{T}(), Lsmall)#, @view(hessian(map_)[1:P-1,1:P-1]))
end

# Convenience function does a dynamic dispatch if we don't start with a MutableFixedSizeVector
# If we have a static arrays dependency, we may also support static dispatches with that?
function AsymptoticPosteriorFD(f, initial_x::AbstractArray{T}) where {T}
    AsymptoticPosteriorFD(f, MutableFixedSizeVector{length(initial_x),T}(initial_x))
end

function AsymptoticPosteriorFD(f, initial_x::MutableFixedSizeVector{P,T}) where {P,T}
    map_ = MAPFD(f, initial_x)
    fit!(map_, initial_x)
    AsymptoticPosteriorFD(f, map_, initial_x)
end
function AsymptoticPosteriorFD(::UndefInitializer, f, initial_x::MutableFixedSizeVector{P,T}) where {P,T}
    map_ = MAPFD(f, initial_x)
    # fit!(map_, initial_x)
    AsymptoticPosteriorFD(f, map_, initial_x)
end


text:
module testsuite

    using Test
    using JuliaDB
    using DataFrames
    using Query
    using CSVFiles
    using Parameters
    #
    # local imports; see ~/.julia/config/startup.jl
    #
    using GlobalDecls
    using Utils
    using ONSCodes
    using DataEditing
    using LAModelData
    using StayingPutSim
    import CareData
    using CareData

    code    = "E06000036"
    name    = "Bracknell Forest"
    regionc = "E12000008"
    regionn = "South East England"

    CCODE_SUSSEX="E10000030"

    @testset "Data Edit Tests" begin

        rc = region_code_from_code( code )
        @test rc == regionc

        rcc3 = region_code_from_name( name )
        @test rcc3 == regionc

        cc = code_from_name( name )
        @test cc == code

        rcc = region_code_from_name( name  )
        @test rcc == regionc

        rcc2 = region_code_from_name( regionn  )
        @test rcc2 == regionc

        rcn = region_name_from_name( name )
        @test rcn == regionn
    end

    @testset "ladata tests" begin
        x=do_exit_rate_query( 2017, "E09000002", 18 )
        exits = get_staying_put_rates( "E09000003", local_authority, true )
    end

    alldata = CareData.load_all()

    carer_dataset = CareData.make_carer_frame(0)
    yp_dataset = CareData.make_yp_frame(0)

    @testset "Population Creation" begin
        settings = CareData.default_data_settings()
        println( "settings=$settings")
        npy = 1000
        no = (2025-2019+1)*npy
        carer = CareData.Carer( -1, 40,0.0,0.1,0.2, 1,3 )
        yp = CareData.YP( 18, 1.0, 2.0, 3.0, OtherEd )
        @time for year in 2019:2025
            for cno in 1:npy
                cn = cno <= 200 ? "C1" : "C2"
                carer.id = cno
                CareData.add_carer_to_frame!( carer_dataset,  year, cn, carer )
                CareData.add_yp_to_frame!( yp_dataset, year, cn, carer.id, yp )
            end
        end
        println( size(carer_dataset))
        println( "no $no")
        # println( carer_dataset )
        @assert size( carer_dataset )[1] == no
        @time for year in 2019:2025
            carers = CareData.get_carers( carer_dataset, "C1", year )
            @test size( carers )[1] == 200
            for c in carers
                yp = CareData.get_yp( yp_dataset, year, c )
                if( c.id % 10 ) == 0
                    @test yp.age == 18
                end
            end
        end
        # try with JuliaDB and indexing
        tcarer_dataset = JuliaDB.table( carer_dataset, pkey=[:year,:id])
        typ_dataset = JuliaDB.table( yp_dataset, pkey=[:year,:carer])
        @time for year in 2019:2025
            carers = CareData.get_carers( tcarer_dataset, "C1", year )
            @test size( carers )[1] == 200
            for c in carers
                yp = CareData.get_yp2( typ_dataset, year, c )
                if( c.id % 200 ) == 0
                    @test yp.age == 18
                end
            end
        end

        for t = 1:3
            carer_data = CareData.make_carer_frame(0)
            yp_data = CareData.make_yp_frame(0)
            ladata = Dict()
            if t == 1
                stayingrates = [0.0,0.0,0.0]
            elseif t == 2
                stayingrates = [0.8,0.5,0.8]
            elseif t == 3
                stayingrates = [1.0,1.0,1.0]
            end
            pid = 0
            for year in 2019:2030
                for i in 1:1000
                    pid += 1
                    CareData.add_aged_data!(
                        yp_dataset = yp_data,
                        carer_dataset = carer_data,
                        pid=pid,
                        ccode="COUNCIL01",
                        ladata = Dict(),
                        stayingrates=stayingrates,
                        startyear=year
                    )
                end # add 1000
            end # years
            len = size( carer_data )[1]
            println( "case: $t  len ")
            if t == 3
                @assert size( carer_data )[1]==3*12000
                @assert size( yp_data )[1]==3*12000
            elseif t == 2
                @assert size( carer_data )[1]>=0
                @assert size( yp_data )[1]>=0
            elseif t == 1
                @assert size( carer_data )[1]==0
                @assert size( yp_data )[1]==0
            end
        end # t loop

    end # popn creation testsa

    @testset "Calculation tests" begin

        carer = CareData.Carer( -1, 40,0.0,0.1,0.2, 1,3 )
        yp = CareData.YP( 18, 1.0, 2.0, 3.0, OtherEd )

        cq = CareData.load_all()
        council = cq.ofdata[50,:]
        println(typeof( council ))
        println( council.rcode )
        outcomes = CarerOutcomes( 100.0, 50.0, 25.0, 150.0 )
        baseoutcomes = CarerOutcomes( 100.0, 50.0, 25.0, 150.0 )
        params = get_default_params()
        settings = default_data_settings()
        year = 2019
        StayingPutSim.override_outcomes!(
            outcomes,
            council.ccode,
            year,#         :: Integer,
            carer, #       :: Carer,
            yp, #           :: YP,.yp_contrib_type = no_contribution
            council, # :: DataFrameRow,
            params ) #       :: Params )
        @test outcomes ≈ baseoutcomes
        baseoutcomes = CarerOutcomes( 100.0, 50.0, 25.0, 150.0 )
        params = get_default_params()
        settings = default_data_settings()
        year = 2019
        StayingPutSim.override_outcomes!(
            outcomes,
            council.ccode,
            year,#         :: Integer,
            carer, #       :: Carer,
            yp, #           :: YP,.yp_contrib_type = no_contribution
            council, # :: DataFrameRow,
            params ) #       :: Params )
        params.yp_contrib_type = no_contribution
        println( "outcomees before no_yp_contrib=$outcomes")
        StayingPutSim.override_outcomes!(
            outcomes,
            council.ccode,
            year,#         :: Integer,
            carer, #       :: Carer,
            yp, #           :: YP,.yp_contrib_type = no_contribution
            council, # :: DataFrameRow,
            params ) #       :: Params )

        println( "outcomes after no_yp_contrib=$outcomes")
        @test outcomes.contributions_from_yp == 0
        @test outcomes.income_recieved == baseoutcomes.income_recieved
        @test outcomes.payments_from_la == baseoutcomes.income_recieved

        params.yp_contrib_type = flat_rate
        params.yp_contrib = 10.0
        StayingPutSim.override_outcomes!(
            outcomes,
            council.ccode,
            year,#         :: Integer,
            carer, #       :: Carer,
            yp, #           :: YP,.yp_contrib_type = no_contribution
            council, # :: DataFrameRow,
            params ) #       :: Params )

        @test outcomes.contributions_from_yp == 10
        @test outcomes.income_recieved == baseoutcomes.income_recieved
        @test outcomes.payments_from_la == baseoutcomes.income_recieved-10.0

        params = get_default_params()
        params.payment = min_payment
        rcode = pay_class_from_region_code( council.rcode )
        StayingPutSim.override_outcomes!(
            outcomes,
            council.ccode,
            year,#         :: Integer,
            carer, #       :: Carer,
            yp, #           :: YP,.yp_contrib_type = no_contribution
            council, # :: DataFrameRow,
            params ) #       :: Params )

        println( "outcomees after no_yp_contrib=$outcomes")
        @test outcomes.income_recieved==194.0

        # payments_from_la :: Real
        # contributions_from_yp :: Real
        # housing_cont :: Real
        # income_recieved :: Real

    end

end


text:
function γ(x::AbstractFloat) # faster if only one wrap can occur
    if x < -π
        x+typeof(x)(2π)
    elseif x > π
        x-typeof(x)(2π)
    else
        x
    end
end

getdimoffsets(A) = (1, cumprod(size(A)[1:end-1])...)


text:
export InertialFrame, SynodicFrame, convert_to_frame

#------------------#
# REFERENCE FRAMES #
#------------------#

struct InertialFrame <: Abstract_ReferenceFrame end
struct SynodicFrame{IsNormalised} <: Abstract_ReferenceFrame
    SynodicFrame(isnormalised::Bool=true) = new{isnormalised}()
end

default_reference_frame(model::Abstract_DynamicalModel) = error("$(model) has no default reference frame assigned!")
default_synodic_reference_frame(::Abstract_DynamicalModel) = SynodicFrame()

#----------------------------------#
# CONVERT BETWEEN REFERENCE FRAMES #
#----------------------------------#
convert_to_frame(state::State, frame::Abstract_ReferenceFrame) = error("Cannot convert $(nameof(typeof(state))){$(state.model),$(state.frame)} to $(frame).")
convert_to_frame(state::State{<:Abstract_DynamicalModel,T}, ::T) where {T<:Abstract_ReferenceFrame} = state
convert_to_frame(traj::Trajectory{<:Abstract_DynamicalModel,T}, ::T) where {T<:Abstract_ReferenceFrame} = traj
function convert_to_frame(traj::Trajectory, frame::Abstract_ReferenceFrame)
    prob0 = traj.prob
    times = traj.t

    u1 = deepcopy(traj.u)
    convert_u!(u1, times, traj, frame)

    # TODO: Convert this into a type so that we can make an interp_summary(::interpType) dispatch (to print out converted sols)
    function interp_and_convert(times,idxs,deriv::Type{Val{0}},p,continuity::Symbol=:left)
        u = deepcopy(traj.interp(times, idxs, deriv, p, continuity))
        convert_u!(u, times, traj, frame)
        return u
    end

    new_sol = DiffEqBase.build_solution(prob0, traj.alg, times, u1; interp=interp_and_convert, retcode=traj.retcode)

    return Trajectory(traj.model, frame, new_sol)
end
function convert_u!(u::AbstractArray, times::AbstractArray, traj, frame)
    for (i, t) in enumerate(times)
        convert_u!(u[i], t, traj, frame)
    end
end
function convert_u!(u, t, traj, frame)
    prob1 = remake(traj.sol.prob; u0=u, tspan=(t, t))
    new_state = convert_to_frame(State(traj.model, traj.frame, prob1), frame)
    u .= new_state.prob.u0
end

const CRASHED_RETCODE = :Crashed
function collision(system::Abstract_DynamicalModel, body, dist=bodvrd(String(body), "RADII")[1]; radii=1., interp_points=10)
    diam = radii * dist
    ContinuousCallback((integrator) -> terminate!(integrator, CRASHED_RETCODE); interp_points) do u, t, integrator
        check_distance(u, t, system, body, diam)
    end
end
check_distance(u, t, system::Abstract_DynamicalModel, body) = error("check_distance not defined for $(nameof(typeof(system)))")
crashed(sol::Trajectory) = sol.sol.retcode == CRASHED_RETCODE
crashed(sol::ODESolution) = sol.retcode == CRASHED_RETCODE


text:
# ==================== struct of model parameters ===============
export get_gamma_profile,
       Mvz,
       Rvz,
       Mvx,
       Rvx,
       Mpz,
       Rpz,
       Mpx,
       Rpx,
       TdParams,
       show,
       model_padding,
       model_smooth,
       density2buoyancy,
       vp2lambda,
       finite_difference_coefficients,
# ==================== Snapshot =================
       Snapshot,
       Wavefield,
       sample_spt2wfd,
       sample_adjoint2pre!,
       snapshot_header,
       append_one_snapshot,
       read_one_snapshot,
# ==================== Wavefield =================
       wavefield_header,
       append_one_wavefield,
       write_wavefield,
       read_one_wavefield,
# ==================== pressure =================
       pressure_header,
       append_one_pressure,
       read_one_pressure,
       copy_snapshot!,
       copy_wavefield!,
# ==================== Boundary =================
       WavefieldBound,
       boundary_header,
       write_boundary,
       append_one_boundary,
       read_boundary,
       read_one_boundary!,
       l2norm_snapshot,
       minus_snapshot,
       add_snapshot,
       reverse_order,
# # ==================== Source =================
       Source,
       ricker,
       convert2miniphase,
       tapered_sinc,
       scaled_sinc,
       get_multi_sources,
       add_source!,
       subtract_source!,
       time_range_source,
       write_source,
       read_source,
       source_isequal,
# # ==================== Recordings ===============
       Recordings,
       write_recordings,
       read_recordings,
       sample_spt2rec!,
       inject_rec2spt!,
       recordings_isequal,
# ============= forward time stepping =============
       A_mul_b!,
       At_mul_b!,
# ============= forward time stepping =============
       vertical_partial_derivative!,
       horizontal_partial_derivative!,
       one_step_forward!,
       multi_step_forward,
       pressure_reconstruct_forward,
# ============= adjoint time stepping =============
       one_step_adjoint!,
       multi_step_adjoint,
       one_step_backward!,
       pressure_reconstruct_backward

include("pml.jl")
include("fdstencil.jl")
include("tdparams.jl")
include("snapshot.jl")
include("source.jl")
include("recordings.jl")
include("spmatvec.jl")
include("forward.jl")
include("adjoint.jl")

# # ==================== struct of model parameters ===============
# export ModelParams,
#        show,
#        model_padding,
#        model_smooth,
#        finite_difference_coefficients,
# # ==================== PML ======================
#        PMLCoefficients,
# # ==================== FD stencil ===============
#        ObsorbFDStencil,
#        RigidFDStencil,
# # ==================== Snapshot =================
#        Snapshot,
#        Wavefield,
#        sample_spt2wfd,
#        snapshot_header,
#        append_one_snapshot,
#        read_one_snapshot,
# # ==================== Wavefield =================
#        wavefield_header,
#        append_one_wavefield,
#        write_wavefield,
#        read_one_wavefield,
# # ==================== pressure =================
#        pressure_header,
#        append_one_pressure,
#        read_one_pressure,
#        copy_snapshot!,
#        copy_wavefield!,
# # ==================== Boundary =================
#        WavefieldBound,
#        boundary_header,
#        write_boundary,
#        append_one_boundary,
#        read_boundary,
#        read_one_boundary,
#        l2norm_snapshot,
#        minus_snapshot,
#        add_snapshot,
# # # ==================== Source =================
#        Source,
#        ricker,
#        convert2miniphase,
#        tapered_sinc,
#        scaled_sinc,
#        get_multi_sources,
#        add_source!,
#        add_multi_sources!,
#        subtract_source!,
#        time_range_multisources,
# # # ==================== Recordings ===============
#        Recordings,
#        write_recordings,
#        read_recordings,
#        sample_spt2rec!,
#        inject_rec2spt!,
# # ============= forward time stepping =============
#        A_mul_b!,
#        At_mul_b!,
# # ============= forward time stepping =============
#        one_step_forward!,
#        one_step_backward!,
#        multi_step_forward,
#        multi_step_forward!,
#        get_boundary_wavefield,
#        pressure_reconstruct_forward,
# # ============= adjoint time stepping =============
#        one_step_adjoint!,
#        one_step_backward,
#        multi_step_adjoint,
#        pressure_reconstruct_backward
#
#
# include("modelparams.jl")
# include("pml.jl")
# include("fdstencil.jl")
# include("snapshot.jl")
# include("source.jl")
# include("recordings.jl")
# include("spmatvec.jl")
# include("forward.jl")
# include("adjoint.jl")


text:
"""
    Biquad(filterType::Symbol, f, dbGain, Q)

Return peak, low shelf or high shelf biquad depending on respective `filterType`
 `:PK`, `:LS` or `:HS`. Parameters are center frequency `f` in Hz,
 Gain `dbGain` in dB and the quality factor `Q`.

The sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function Biquad(filterType::Symbol, f, dbGain, Q)
    if filterType == :PK
        return PK(f, dbGain, Q)
    elseif filterType == :LS
        return LS(f, dbGain, Q)
    elseif filterType == :HS
        return HS(f, dbGain, Q)
    else
        error("three-argument filterType not recognized")
    end
end

"""
    Biquad(filterType::Symbol, f, Q)

Return lowpass, highpass, bandpass, notch or allpass biquad depending on
respective `filterType` `:LP`, `:HP`, `:BP`, `:NO` or `:AP`.
Parameters are frequency `f` in Hz and the quality factor `Q`.

The sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function Biquad(filterType::Symbol, f, Q)
    if filterType == :LP
        return LP(f, Q)
    elseif filterType == :HP
        return HP(f, Q)
    elseif filterType == :BP
        return BP(f, Q)
    elseif filterType == :NO
        return NO(f, Q)
    elseif filterType == :AP
        return AP(f,Q)
    else
        error("two-argument filterType not recognized")
    end
end

"""
    Biquad(filterType::Symbol, dbGain)

When `filterType` is `:G` for Gain return the field level ratio given by the
parameter `dbGain` in dB. Throw error for other `filterType`.
"""
function Biquad(filterType::Symbol, dbGain)
    if filterType == :G
        return 10^(dbGain/20)
    else
        error("single-argument filterType expects :G for Gain in dB")
    end
end

"""
    Biquad(x::Tuple{Symbol, Vararg})

Return Biquad of a filter specified in a tuple format, where the first argument
is a symbol and defines the filter type:
`([:LP, :HP, :BP, :NO, :AP, :PK, :LS, :HS, :G], f, dbGain, Q)`

For non-gain filter types [:LP, :HP, :BP, :NO, :AP] the argument `dbGain` is to
be omitted. For the simple gain-filter `:G` the arguments `f` and `Q` shall be
omitted.

`f` is in Hz and `dbGain` in dB. The quality factor `Q` is dimensionless.

The filter types refert to a lowpass, highpass, bandpass, notch, allpass, peaking,
low shelf, high shelf and gain filter respectively.

The sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function Biquad(x::Tuple{Symbol, Vararg})
    Biquad(x...)
end

"""
    SecondOrderSections(x::Tuple{Symbol, Vararg})

Return a filter specified in a tuple format as SecondOrderSections, where the first argument
is a symbol and defines the filter type:
`([:LP, :HP, :BP, :NO, :AP, :PK, :LS, :HS, :G], f, dbGain, Q)`

For non-gain filter types [:LP, :HP, :BP, :NO, :AP] the argument `dbGain` is to
be omitted. For the simple gain-filter `:G` the arguments `f` and `Q` shall be
omitted.

`f` is in Hz and `dbGain` in dB. The quality factor `Q` is dimensionless.

The filter types refert to a lowpass, highpass, bandpass, notch, allpass, peaking,
low shelf, high shelf and gain filter respectively.

The sampling rate defaults to 48 kHz or can be
set with [`EqualizerFilters.sampling_rate`](@ref).
"""
function SecondOrderSections(x::Tuple{Symbol, Vararg})
    SecondOrderSections(Biquad(x...))
end

"""
    SecondOrderSections(x::Array{T} where {T <: Tuple{Symbol, Vararg}})

Return SecondOrderSections that represent the filters defined in the array
in a tuple format. For details of the tuple format see
[`Biquad(x::Tuple{Symbol, Vararg})`](@ref).
"""
function SecondOrderSections(x::Array{T} where {T <: Tuple{Symbol, Vararg}})
    biquads = [Biquad(x[i]...) for i in 1:length(x)]
    *(biquads...)
end


text:
@reexport module _Math

import Nemo: add!, sub!, mul!

export
    checked_add, checked_sub,
    setbits,
    @activate_mathematica,
    @activate_matlab

using Reexport: @reexport

checked_add(a::Number, b::Number) = checked_add(promote(a, b)...)
checked_add(a::T, b::T) where {T<:Number} = Base.checked_add(a, b)
checked_add(a::BigInt, b::BigInt) = a + b
checked_add(a::BigFloat, b::BigFloat) = a + b

checked_sub(a::Number, b::Number) = checked_sub(promote(a, b)...)
checked_sub(a::T, b::T) where {T<:Number} = Base.checked_sub(a, b)
checked_sub(a::BigInt, b::BigInt) = a - b
checked_sub(a::BigFloat, b::BigFloat) = a - b

setbits(n::Integer, idx) = begin
    m = copy(n)
    for i in idx
        m |= one(n) << i
    end
    m
end

macro activate_mathematica(); Expr(:using, :Mathematica); end
macro activate_matlab(); Expr(:using, :MATLAB); end

include("gmp.jl")
include("mpfr.jl")
include("Series/series.jl")
include("GroupTheory/grouptheory.jl")
include("NumberTheory/numbertheory.jl")
include("Combinatorics/combinatorics.jl")
include("Algebra/algebra.jl")
include("EllipticCurves/ellipticcurves.jl")

end


text:
using PalmerPenguins
using CSV
using DataFrames
using Tables

using Dates
using Test

# Always accept download in CI tests
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

@testset "PalmerPenguins.jl" begin
    # column names
    colnames_simplified = [
        "species",
        "island",
        "bill_length_mm",
        "bill_depth_mm",
        "flipper_length_mm",
        "body_mass_g",
        "sex",
    ]
    coleltypes_simplified = [
        String,
        String,
        Union{Missing,Float64},
        Union{Missing,Float64},
        Union{Missing,Int},
        Union{Missing,Int},
        Union{Missing,String},
    ]
    firstrow_simplified = ["Adelie", "Torgersen", 39.1, 18.7, 181, 3750, "male"]

    colnames_raw = [
        "studyName",
        "Sample Number",
        "Species",
        "Region",
        "Island",
        "Stage",
        "Individual ID",
        "Clutch Completion",
        "Date Egg",
        "Culmen Length (mm)",
        "Culmen Depth (mm)",
        "Flipper Length (mm)",
        "Body Mass (g)",
        "Sex",
        "Delta 15 N (o/oo)",
        "Delta 13 C (o/oo)",
        "Comments",
    ]
    coleltypes_raw = [
        String,
        Int,
        String,
        String,
        String,
        String,
        Bool,
        Date,
        Union{Missing,Float64},
        Union{Missing,Float64},
        Union{Missing,Int},
        Union{Missing,Int},
        Union{Missing,String},
        Float64,
        Float64,
        String,
    ]
    firstrow_raw = [
        "PAL0708",
        1,
        "Adelie Penguin (Pygoscelis adeliae)",
        "Anvers",
        "Torgersen",
        "Adult, 1 Egg Stage",
        "N1A1",
        true,
        Date(2007, 11, 11),
        39.1,
        18.7,
        181,
        3750,
        "MALE",
        missing,
        missing,
        "Not enough blood for isotopes.",
    ]

    @testset "load" begin
        @testset "simplified" begin
            table = PalmerPenguins.load()
            table2 = PalmerPenguins.load(; raw=false)
            for i in 1:length(table)
                @test all(table[i] .=== table2[i])
            end

            # Check some properties
            @test table isa CSV.File
            @test length(table) == 344
            @test Tables.columnnames(Tables.columns(table)) == Symbol.(colnames_simplified)
            for (name, T) in zip(colnames_simplified, coleltypes_simplified)
                eltype(Tables.getcolumn(Tables.columns(table), Symbol(name))) === T
            end

            # Check first row
            firstrow = first(table)
            for i in 1:length(firstrow)
                @test isequal(firstrow[i], firstrow_simplified[i])
            end
        end

        @testset "raw" begin
            table = PalmerPenguins.load(; raw=true)

            # Check some properties
            @test table isa CSV.File
            @test length(table) == 344
            @test Tables.columnnames(Tables.columns(table)) == Symbol.(colnames_raw)
            for (name, T) in zip(colnames_raw, coleltypes_raw)
                eltype(Tables.getcolumn(Tables.columns(table), Symbol(name))) === T
            end

            # Check first row
            firstrow = first(table)
            for i in 1:length(firstrow)
                @test isequal(firstrow[i], firstrow_raw[i])
            end
        end
    end

    @testset "DataFrames" begin
        @testset "simplified" begin
            df = DataFrame(PalmerPenguins.load())
            df2 = DataFrame(PalmerPenguins.load(; raw=false))
            for i in 1:size(df, 2)
                @test all(df[!, i] .=== df2[!, i])
            end

            # Check some properties
            @test df isa DataFrame
            @test nrow(df) == 344
            @test names(df) == colnames_simplified
            for (name, T) in zip(colnames_simplified, coleltypes_simplified)
                eltype(df[!, name]) === T
            end

            # Check first row
            firstrow = df[1, :]
            for i in 1:length(firstrow)
                @test isequal(firstrow[i], firstrow_simplified[i])
            end
        end

        @testset "raw" begin
            df = DataFrame(PalmerPenguins.load(; raw=true))

            # Check some properties
            @test df isa DataFrame
            @test nrow(df) == 344
            @test names(df) == colnames_raw
            for (name, T) in zip(colnames_raw, coleltypes_raw)
                eltype(df[!, name]) === T
            end

            # Check first row
            firstrow = df[1, :]
            for i in 1:length(firstrow)
                @test isequal(firstrow[i], firstrow_raw[i])
            end
        end
    end
end


text:
using OffsetArrays
using Test
ambs = detect_ambiguities(Base, Core)  # in case these have ambiguities of their own
using PaddedViews
using PaddedViews: filltype
@testset "ambiguities" begin
    @test isempty(setdiff(detect_ambiguities(PaddedViews, Base, Core), ambs))
end

if VERSION >= v"1.2"
    # array summary changes after v1.2
    using Documenter
    doctest(PaddedViews, manual = false)
end

@testset "PaddedView" begin
    for n = 0:5
        a = @inferred(PaddedView(0, ones(Int,ntuple(d->1,n)), ntuple(x->x+1,n)))
        @test axes(a) == ntuple(x->1:x+1,n)
        @test axes(a, ndims(a)+1) == Base.OneTo(1)
        @test @inferred(a[1]) == 1
        n > 0 && @test @inferred(a[2]) == 0
        @test @inferred(a[ntuple(x->1,n)...]) == 1
        n > 0 && @test @inferred(a[2, ntuple(x->1,n-1)...]) == 0
    end
    a0 = reshape([3])
    a = @inferred(PaddedView(-1, a0, ()))
    @test axes(a) == ()
    @test ndims(a) == 0
    @test a[] == 3

    a = reshape(1:9, 3, 3)
    A = @inferred(PaddedView(0.0, a, (Base.OneTo(4), Base.OneTo(5))))
    @test eltype(A) == Int
    @test ndims(A) == 2
    @test size(A) === (4,5)
    @test @inferred(axes(A)) === (Base.OneTo(4), Base.OneTo(5))
    @test @inferred(axes(A, 3)) === Base.OneTo(1)
    @test A == [1 4 7 0 0;
                2 5 8 0 0;
                3 6 9 0 0;
                0 0 0 0 0]
    A32 = @inferred(PaddedView{Float32}(0.0, a, (Base.OneTo(4), Base.OneTo(5))))
    @test eltype(A32) == Float32
    @test A32 == A
    @test A32[1, 1] === 1.0f0

    A = @inferred(PaddedView(0.0, a, (0:4, -1:5)))
    @test eltype(A) == Int
    @test ndims(A) == 2
    @test size(A) == (5, 7)
    @test @inferred(axes(A)) === (0:4, -1:5)
    @test @inferred(axes(A, 3)) === 1:1
    @test A == OffsetArray([0 0 0 0 0 0 0;
                            0 0 1 4 7 0 0;
                            0 0 2 5 8 0 0;
                            0 0 3 6 9 0 0;
                            0 0 0 0 0 0 0], 0:4, -1:5)
    A32 = @inferred(PaddedView{Float32}(0.0, a, (0:4, -1:5)))
    @test eltype(A32) == Float32
    @test A32 == A
    @test A32[1, 1] === 1.0f0

    A = @inferred(PaddedView(0.0, a, (Base.OneTo(5), Base.OneTo(5)), (2:4, 2:4)))
    @test A == [0 0 0 0 0;
                0 1 4 7 0;
                0 2 5 8 0;
                0 3 6 9 0;
                0 0 0 0 0]
    A32 = @inferred(PaddedView{Float32}(0.0, a, (Base.OneTo(5), Base.OneTo(5)), (2:4, 2:4)))
    @test eltype(A32) == Float32
    @test A32 == A
    @test A32[2, 2] === 1.0f0

    @test A == @inferred(PaddedView(0.0, a, (5, 5), (2, 2)))
    A32 = @inferred(PaddedView{Float32}(0.0, a, (5, 5), (2, 2)))
    @test eltype(A32) == Float32
    @test A32 == A
    @test A32[2, 2] === 1.0f0

    B = @inferred(PaddedView(0.0, OffsetArray(a, (2:4, 2:4)), (Base.OneTo(5), Base.OneTo(5))))
    @test B == A

    @test PaddedView(0, ones(4), (3,), (1,)) == PaddedView(0, ones(4), (3,))

    let a = reshape(1:6, 2, 3) # [1 3 5; 2 4 6]
        @test @inferred(PaddedView(-1, a, (2, 2))) == [1 3
                                                       2 4]

        @test @inferred(PaddedView(-1, a, (1, 4)))  == [1 3 5 -1 ]

        @test @inferred(PaddedView(-1, a, (2, 4), (1, 1)))  == [1 3 5 -1
                                                                2 4 6 -1]

        @test @inferred(PaddedView(-1, a, (2, 4), (1, 2)))  == [-1 1 3 5
                                                                -1 2 4 6]

        @test @inferred(PaddedView(-1, a, (2, 3), (0, 2)))  == [-1  2  4
                                                                -1 -1 -1]

        @test @inferred(PaddedView(-1, a, (1, 4), (1, -1))) == [5 -1 -1 -1]

        # Create an OffsetArray with axes (0:1, 2:4)
        oa = OffsetArray(a, -1, 1)
        # Make a PaddedView of its elements (1:3, 1:2) (i.e pv[1, 1] == oa[1, 1] == "a[2, 0]" == -1)
        pv = @inferred PaddedView(-1, oa, (3, 2))
        @test axes(pv) === (Base.OneTo(3), Base.OneTo(2))
        @test pv[1:3, 1:2] == [-1 2; -1 -1; -1 -1]
        # Put oa[1, 1] in pv[2, 0] and make a PaddedView of the resulting (1:3, 1:2) elements
            # i.e. same as putting oa[0,2] in pv[1, 1]
        pv = @inferred PaddedView(-1, oa, (3, 2), (2, 0))
        @test pv[1:3, 1:2] == [1 3; 2 4; -1 -1]
    end
end

@testset "paddedviews" begin
    @testset "0-d array" begin
        a = reshape([1])
        pa, pb = paddedviews(-1, a, a)
        @test pa === pb
        @test pa == a
    end

    a1 = reshape([1,2], 2, 1)
    a2 = [1.0,2.0]'
    a1p, a2p = @inferred(paddedviews(0, a1, a2))
    @test a1p == [1 0; 2 0]
    @test a1p[CartesianIndices(a1)] == a1
    @test a2p == [1.0 2.0; 0.0 0.0]
    @test a2p[CartesianIndices(a2)] == a2
    @test eltype(a1p) == Int
    @test eltype(a2p) == Float64
    @test axes(a1p) === axes(a2p) === (Base.OneTo(2), Base.OneTo(2))

    a3 = OffsetArray([1.0,2.0]', (0,-1))
    a1p, a3p = @inferred(paddedviews(0, a1, a3))
    @test a1p == OffsetArray([0 1; 0 2], 1:2, 0:1)
    @test a1p[CartesianIndices(a1)] == a1
    @test a3p == OffsetArray([1.0 2.0; 0.0 0.0], 1:2, 0:1)
    @test a3p[CartesianIndices(a3)] == a3
    @test eltype(a1p) == Int
    @test eltype(a3p) == Float64
    @test axes(a1p) === axes(a3p) === (1:2, 0:1)

    @test @inferred(paddedviews(3)) == ()
    @test_throws ErrorException PaddedViews.outerinds()
    # But a zero-dimensional input should not trigger that error
    a = reshape([5])
    @test paddedviews(-1, a) == (a,)
end

@testset "sym_paddedviews" begin
    @testset "0-d array" begin
        a = reshape([1])
        pa, pb = sym_paddedviews(-1, a, a)
        @test pa === pb
        @test pa == a
    end

    # even case
    a1 = reshape([1,2], 2, 1)
    a2 = [1.0,2.0]'
    a1p, a2p = @inferred(sym_paddedviews(0, a1, a2))
    @test a1p == [1 0; 2 0]
    @test a1p[CartesianIndices(a1)] == a1
    @test a2p == [1.0 2.0; 0.0 0.0]
    @test a2p[CartesianIndices(a2)] == a2
    @test eltype(a1p) == Int
    @test eltype(a2p) == Float64
    @test axes(a1p) === axes(a2p) === (1:2, 1:2)

    a1 = reshape([1,2,3], 3, 1)
    a2 = [1.0,2.0,3.0]'
    a1p, a2p = @inferred(sym_paddedviews(0, a1, a2))

    @test a1p == OffsetArray([0 1 0; 0 2 0; 0 3 0], (1:3, 0:2))
    @test a1p[CartesianIndices(a1)] == a1
    @test a2p == OffsetArray([0.0 0.0 0.0; 1.0 2.0 3.0; 0.0 0.0 0.0], (0:2, 1:3))
    @test a2p[CartesianIndices(a2)] == a2
    @test eltype(a1p) == Int
    @test eltype(a2p) == Float64
    @test axes(a1p) === (1:3, 0:2)
    @test axes(a2p) === (0:2, 1:3)

    a3 = OffsetArray([1.0,2.0,3.0]', (0,-1))
    a1p, a3p = @inferred(sym_paddedviews(0, a1, a3))

    @test a1p == OffsetArray([0 1 0; 0 2 0; 0 3 0], 1:3, 0:2)
    @test a1p[CartesianIndices(a1)] == a1
    @test a3p == OffsetArray([0.0 0.0 0.0; 1.0 2.0 3.0; 0.0 0.0 0.0], 0:2, 0:2)
    @test a3p[CartesianIndices(a3)] == a3
    @test eltype(a1p) == Int
    @test eltype(a3p) == Float64
    @test axes(a1p) == (1:3, 0:2)
    @test axes(a3p) == (0:2, 0:2)

    @test @inferred(sym_paddedviews(3)) == ()
end

@testset "showarg" begin
    a = collect(reshape(1:9, 3, 3))
    pv = PaddedView(-1, a, (0:4, 1:3))
    io = IOBuffer()
    show(IOContext(io, :displaysize=>(1000,1000)), MIME("text/plain"), pv)
    str = String(take!(io))
    @test endswith(str, "PaddedView(-1, ::$(Array{Int,2}), (0:4, 1:3)) with eltype $Int with indices 0:4×1:3:\n -1  -1  -1\n  1   4   7\n  2   5   8\n  3   6   9\n -1  -1  -1")
end

@testset "similar" begin
    for (T, v) in ((Int, 0),
                   (Int, 0.),
                   (Float64, 0.),
                   (Missing, missing),
                   (Nothing, nothing))
        A = reshape(1:9, 3, 3)
        Ap = @inferred(PaddedView(v, A, (0:4, 0:4)))
        B = similar(Ap)
        @test eltype(B) == eltype(Ap)
        @test size(B) == size(Ap)
        @test axes(B) == axes(Ap)

        B = similar(Ap, (4, 4))
        @test eltype(B) == eltype(Ap)
        @test size(B) == (4, 4)
        @test axes(B) == (Base.OneTo(4), Base.OneTo(4))

        B = similar(Ap, (4, 4))
        @test eltype(B) == eltype(Ap)
        @test size(B) == (4, 4)
        @test axes(B) == (Base.OneTo(4), Base.OneTo(4))

        B = similar(Ap, Int, (4, 4))
        @test eltype(B) == Int
        @test size(B) == (4, 4)
        @test axes(B) == (Base.OneTo(4), Base.OneTo(4))
    end
end

@testset "nothing/missing" begin
    for (FT, T) in ((Missing, Float32),
                    (Nothing, Float32),
                    (Float32, Nothing),
                    (Float32, Missing),
                    (Nothing, Missing),
                    (Missing, Nothing))
        @test @inferred(filltype(FT, T)) === Union{FT, T}
        @test @inferred(filltype(T, Union{FT, T})) === Union{FT, T}
        @test @inferred(filltype(FT, Union{FT, T})) === Union{FT, T}
    end

    for (T, v) in ((Missing, missing),
                 (Nothing, nothing))
        A = reshape(1:9, 3, 3)
        Ap = @inferred(PaddedView(v, A, (0:4, 0:4)))
        @test axes(Ap) == (OffsetArrays.IdOffsetRange(0:4), OffsetArrays.IdOffsetRange(0:4))
        @test eltype(Ap) === Union{T, eltype(A)}
        @test Ap[0, 0] === v
        @test Ap[1, 1] === 1
        @test Ap[axes(A)...] == A

        Ap = @inferred(PaddedView(v, A, (5, 5)))
        @test axes(Ap) === (Base.OneTo(5), Base.OneTo(5))
        @test eltype(Ap) === Union{T, eltype(A)}
        @test Ap[4, 4] === v
        @test Ap[1, 1] === 1
        @test Ap[axes(A)...] == A
    end
end

@testset "getindex/setindex!" begin
    @testset "1d" begin
        A = collect(1:3)
        Ap = PaddedView(missing, A, (0:4, ))
        @test axes(Ap) == (0:4,)
        @test ismissing(Ap[0])
        @test ismissing(first(Ap))
        @test ismissing(last(Ap))
        @test Ap[1] == 1
        @test_throws ArgumentError Ap[0] = 1
        Ap[1] = 10
        @test Ap[1] == 10

        # test setindex! result
        @test setindex!(Ap, 20, 1) === Ap
        @test 10 === (Ap[1] = 10)

        A = collect(1:3)
        Ap = PaddedView(missing, A, (4:6, ))
        @test axes(Ap) == (4:6, )
        @test all(ismissing, Ap)
        @test_throws BoundsError Ap[1] = 1

        A = collect(1:3)
        Ap = PaddedView(missing, OffsetArray(A, -1), (0:4, ))
        @test axes(Ap) == (0:4, )
        @test first(Ap) == Ap[0] === 1
        @test ismissing(Ap[4])
        Ap[0] = 10
        @test Ap[0] == 10
        @test_throws ArgumentError Ap[4] = 1
    end

    @testset "2d" begin
        A = reshape(collect(1:6), 2, 3)
        Ap = PaddedView(missing, A, (0:4, 0:4))
        @test axes(Ap) == (0:4, 0:4)
        @test ismissing(Ap[0, 0])
        @test ismissing(Ap[2])
        @test Ap[7] == 1
        Ap[1, 1] = 10
        @test Ap[7] == Ap[1, 1] == 10

        A = reshape(collect(1:6), 2, 3)
        Ap = PaddedView(missing, A, (3:4, 4:6))
        @test axes(Ap) == (3:4, 4:6)
        @test all(ismissing, collect(Ap))
        @test_throws BoundsError Ap[1, 1] = 1

        A = reshape(collect(1:6), 2, 3)
        Ap = PaddedView(missing, OffsetArray(A, -1, -1), (0:4, 0:4))
        @test axes(Ap) == (0:4, 0:4)
        @test Ap[0, 0] === 1
        @test ismissing(Ap[4, 4])
        @test Ap[6] == 3
        Ap[0, 0] = 10
        @test Ap[0, 0] == 10
        @test_throws ArgumentError Ap[4, 4] = 1
    end
end


text:
# functions to compare unrooted networks
# Claudia & Cecile November 2015

# Traverses a tree postorder and modifies matrix M
# with edges as rows and species as columns (see tree2Matrix)
# S should be sorted --fixit: are you sure?
function traverseTree2Matrix!(node::Node, edge::Edge, ie::Vector{Int}, M::Matrix{Int}, S::Union{Vector{String},Vector{Int}})
    child = getOtherNode(edge,node) # must not be a leaf
    indedge = ie[1]
    M[indedge,1] = edge.number
    ie[1] += 1 # mutable array: to modify edge index 'ie' outside function scope
    for e in child.edge #postorder traversal
        if(!isEqual(e,edge)) # assumes a tree here
            grandchild = getOtherNode(e,child)
            if grandchild.leaf
                indsp = findfirst(isequal(grandchild.name), S)
                indsp != nothing || error("leaf $(grandchild.name) not in species list $(S)")
                M[indedge,indsp+1] = 1 #indsp+1 bc first column is edge numbers
            else
                inde = ie[1];
                traverseTree2Matrix!(child,e,ie,M,S)
                M[indedge,2:size(M,2)] .|= M[inde,2:size(M,2)]
            end
        end
    end
end

# takes a tree and a list of species as input,
# and produces a matrix M with edges as rows and species as columns:
# Mij=1 if species j is descendant of edge i, 0 ow.
# allows for missing taxa:
# Mij=0 if species not present in tree. This is handled in calculateObsCFAll with sameTaxa function
function tree2Matrix(T::HybridNetwork, S::Union{Vector{String},Vector{Int}}; rooted=true::Bool)
    length(T.hybrid)==0 || error("tree2Matrix only works on trees. Network has $(T.numHybrids) hybrid nodes.")
    # sort!(S) # why sort 'taxa', again and again for each tree? Benefits?
    ne = length(T.edge)-T.numTaxa # number of internal branch lengths
    if (T.node[T.root].leaf)      # the root is a leaf: the 1 edge stemming from the root is an external edge
        ne += 1                   # and will need to get a row in the matrix, to be deleted later.
    end
    M = zeros(Int,ne,length(S)+1)
    # M[:,1] = sort!([e.number for e in T.edge])
    ie = [1] # index of next edge to be documented: row index in M
    for e in T.node[T.root].edge
        child = getOtherNode(e,T.node[T.root])
        if (!child.leaf)
            traverseTree2Matrix!(T.node[T.root],e,ie,M,S)
        end
    end
    if (!rooted && length(T.node[T.root].edge)<3)
        # remove first row of M: 1st edge at the root, duplicated edge
        # if the tree is to be considered as unrooted, or just leaf.
        M = M[2:size(M,1),:] # makes a copy, too bad.
    end
    return M
end

"""
`hardwiredClusters(net::HybridNetwork, S::Union{Vector{String},Vector{Int}})`

Returns a matrix describing all the hardwired clusters in a network.
Warnings: Clusters are rooted, so the root must be correct.
          Allows for missing taxa, with entries all 0.

Each row corresponds to one internal edge, that is, external edges are excluded.
If the root is a leaf node, the external edge to that leaf is included (first row).
Both parent hybrid edges to a given hybrid node only contribute a single row (they share the same hardwired cluster).

- first column: edge number
- next columns: 0/1. 1=descendant of edge, 0=not a descendant, or missing taxon.
- last column:  10/11 values. 10=tree edge, 11=hybrid edge
"""
function hardwiredClusters(net::HybridNetwork, S::Union{Vector{String},Vector{Int}})
    ne = length(net.edge)-net.numTaxa # number of internal branch lengths
    ne -= length(net.hybrid)          # to remove duplicate rows for the 2 parent edges of each hybrid
    if (net.node[net.root].leaf)      # root is leaf: the 1 edge stemming from the root is an external edge
        ne += 1                       #               but needs to be included still (rooted clusters).
    end
    M = zeros(Int,ne,length(S)+2)
    ie = [1] # index of next edge to be documented: row index in M
    for e in net.node[net.root].edge
        hardwiredClusters!(net.node[net.root],e,ie,M,S)
    end
    return M
end

function hardwiredClusters!(node::Node, edge::Edge, ie::Vector{Int}, M::Matrix{Int},
                            S::Union{Vector{String},Vector{Int}})
    child = getOtherNode(edge,node)

    !child.leaf || return 0 # do nothing if child is a leaf.

    if (edge.hybrid) # edge is a hybrid. Need to find its partner.
        (edge.isChild1 ? edge.node[1] == child : edge.node[2] == child) || error(
        "inconsistency during network traversal: node $(child.number) supposed to be child of hybrid edge $(edge.number), inconsistent with isChild1.")
        partner = nothing
        partnerVisited = true
        indpartner = 0
        for e in child.edge
            if (e.hybrid && e != edge && (e.isChild1 ? e.node[1] == child : e.node[2] == child))
                partner = e
                indpartner = findfirst(isequal(partner.number), M[:,1])
                if isnothing(indpartner)
                    partnerVisited = false # will need to continue traversal
                end
                break # partner hybrid edge was found
            end
        end
        partner !== nothing || error("partner hybrid edge not found for edge $(edge.number), child $(child.number)")
        !partnerVisited || return indpartner
    end

    indedge = ie[1]
    M[indedge,1] = edge.number
    M[indedge,end] = edge.hybrid ? 11 : 10
    ie[1] += 1 # mutable array

    for e in child.edge # postorder traversal
        if (e != edge && (!edge.hybrid || e!=partner)) # do not go back to (either) parent edge.
            grandchild = getOtherNode(e,child)
            if grandchild.leaf
                indsp = findfirst(isequal(grandchild.name), S)
                indsp != nothing || error("leaf $(grandchild.name) not in species list $(S)")
                M[indedge,indsp+1] = 1 #indsp+1 because first column is edge numbers
            else
                inde = hardwiredClusters!(child,e,ie,M,S)
                M[indedge,2:end-1] .|= M[inde,2:end-1]
            end
        end
    end
    return indedge
end

"""
    hardwiredCluster(edge::Edge,taxa::Union{Vector{String},Vector{Int}})
    hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}})
    hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}},
                      visited::Vector{Int})

Calculate the hardwired cluster of `node`, coded a vector of booleans:
true for taxa that are descendent of nodes, false for other taxa (including missing taxa).

The node should belong in a rooted network for which `isChild1` is up-to-date.
Run `directEdges!` beforehand. This is very important, otherwise one might enter an infinite loop,
and the function does not test for this.

visited: vector of node numbers, of all visited nodes.

# Examples:
```jldoctest
julia> net5 = "(A,((B,#H1),(((C,(E)#H2),(#H2,F)),(D)#H1)));" |> readTopology |> directEdges! ;

julia> taxa = net5 |> tipLabels # ABC EF D
6-element Array{String,1}:
 "A"
 "B"
 "C"
 "E"
 "F"
 "D"

julia> hardwiredCluster(net5.edge[12], taxa) # descendants of 12th edge = CEF
6-element Array{Bool,1}:
 0
 0
 1
 1
 1
 0
```
"""
function hardwiredCluster(edge::Edge,taxa::Union{Vector{String},Vector{Int}})
    v = zeros(Bool,length(taxa))
    hardwiredCluster!(v,edge,taxa)
    return v
end

hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}}) =
    hardwiredCluster!(v,edge,taxa,Int[])

function hardwiredCluster!(v::Vector{Bool},edge::Edge,taxa::Union{Vector{String},Vector{Int}},
                           visited::Vector{Int})
    n = getChild(edge)
    if n.leaf
        j = findall(isequal(n.name), taxa)
        length(j)==1 || error("taxon $(n.name) was not found in taxon list, or more than once")
        v[j[1]]=true
        return nothing
    end
    if n.number in visited
        return nothing  # n was already visited: exit. avoid infinite loop is isChild1 was bad.
    end
    push!(visited, n.number)
    for ce in n.edge
        if n == getParent(ce)
            hardwiredCluster!(v,ce,taxa,visited)
        end
    end
    return nothing
end

"""
    descendants(edge::Edge, internal::Bool=false)

Return the node numbers of the descendants of a given edge:
all descendant nodes if `internal` is true (internal nodes and tips),
or descendant tips only otherwise (defaults).

`edge` should belong in a rooted network for which `isChild1` is up-to-date.
Run `directEdges!` beforehand. This is very important, otherwise one might enter an infinite loop,
and the function does not test for this.

## Examples
```jldoctest
julia> net5 = "(A,((B,#H1),(((C,(E)#H2),(#H2,F)),(D)#H1)));" |> readTopology |> directEdges! ;

julia> PhyloNetworks.descendants(net5.edge[12], true) # descendants of 12th edge: all of them
7-element Array{Int64,1}:
 -6
 -7
  4
  6
  5
 -9
  7

julia> PhyloNetworks.descendants(net5.edge[12]) # descendant leaves only
3-element Array{Int64,1}:
 4
 5
 7
```
"""
function descendants(edge::Edge, internal::Bool=false)
    visited = Int[]
    des = Int[]
    descendants!(des, visited, edge, internal)
    return des
end

function descendants!(des::Vector{Int}, visited::Vector{Int}, edge::Edge, internal::Bool=false)
    n = getChild(edge)
    if n.hybrid # only need to check previous visits for hybrid nodes
        n.number in visited && return nothing
        push!(visited, n.number)
    end
    if internal || n.leaf
        push!(des, n.number)
    end
    for ce in n.edge
        if n == getParent(ce)
            descendants!(des, visited, ce, internal)
        end
    end
    return nothing
end

"""
    isdescendant(des:Node, anc::Node)

Return true if `des` is a strict descendant of `anc`, using `isChild1` fields
to determine the direction of edges. See [`isdescendant_undirected`](@ref)
for a version that does not use `isChild1`.
"""
function isdescendant(des::Node, anc::Node)
    visited = Int[]
    for e in anc.edge
        anc !== getChild(e) || continue # skip parents of anc
        if isdescendant!(visited, des, e)
            return true
        end
    end
    return false
end
function isdescendant!(visited::Vector{Int}, des::Node, e::Edge)
    n = getChild(e)
    if n == des
        return true
    end
    if n.hybrid # only need to check previous visits for hybrid nodes
        if n.number in visited # n & its descendants were already visited: exit
            return false
        end
        push!(visited, n.number)
    end
    for ce in n.edge
        if n == getParent(ce)
            if isdescendant!(visited, des, ce) return true; end
        end
    end
    return false
end

"""
    isdescendant_undirected(des:Node, ancestor::Node, parentedge)

Return `true` if `des` is a strict descendant of `ancestor` when starting
from edge `parentedge` and going towards `ancestor` onward, regardless
of the field `isChild1` of tree edges; `false` otherwise.

This is useful to know how descendant relationships would change as a result
of reverting the direction of a tree edge, without actually
modifying the direction (`isChild1`) of any edge.

`parentedge` should be connected to `ancestor` (not checked).
The direction of hybrid edges is respected (via `isChild1`), that is,
the traversal does not go from the child to the parent of a hybrid edge.
"""
function isdescendant_undirected(des::Node, anc::Node, parentedge::Edge)
    visited = Int[]
    isdescendant_undirected!(visited, des, anc, parentedge)
end
function isdescendant_undirected!(visited::Vector{Int}, des::Node, anc::Node, parentedge::Edge)
    for e in anc.edge
        e !== parentedge || continue # do not go back up where we came from
        !e.hybrid || getParent(e) === anc || continue # do not go back up a parent hybrid edge of anc
        n = getOtherNode(e, anc)
        if n === des
            return true
        end
        if n.hybrid
            !(n.number in visited) || continue # skip to next edge is n already visited
            push!(visited, n.number)
        end
        if isdescendant_undirected!(visited, des, n, e)
            return true
        end
    end
    return false
end

"""
    ladderpartition(tree::HybridNetwork)

For each node in `tree`, calculate the clade below each child edge of the node,
and each clade moving up the "ladder" from the node to the root. The output is a
tuple of 2 vectors (node) of vector (clade) of vectors (taxon in clade):
`below,above`. More specifically, for node number `n`,
`below[n]` is generally of vector of 2 clades: one for the left child and one
for the right child of the node (unless the node is of degree 2 or is a polytomy).
`above[n]` contains the grade of clades above node number `n`.

WARNING: assumes that
1. node numbers and edge numbers can be used as indices, that is,
   be all distinct, positive, covering exactly 1:#nodes and 1:#edges.
2. edges are corrected directed (`isChild1` is up-to-date) and
   nodes have been pre-ordered already (field `nodes_changed` up-to-date).

# examples
```jldoctest
julia> tree = readTopology("(O,A,((B1,B2),(E,(C,D))));");

julia> PhyloNetworks.resetNodeNumbers!(tree; checkPreorder=true, type=:postorder)

julia> printNodes(tree)
node leaf  hybrid hasHybEdge name inCycle edges'numbers
1    true  false  false      O    -1      1
2    true  false  false      A    -1      2
3    true  false  false      B1   -1      3
4    true  false  false      B2   -1      4
8    false false  false           -1      3    4    5
5    true  false  false      E    -1      6
6    true  false  false      C    -1      7
7    true  false  false      D    -1      8
9    false false  false           -1      7    8    9
10   false false  false           -1      6    9    10
11   false false  false           -1      5    10   11
12   false false  false           -1      1    2    11

julia> below, above = PhyloNetworks.ladderpartition(tree);

julia> below
12-element Array{Array{Array{Int64,1},1},1}:
 [[1]]
 [[2]]
 [[3]]
 [[4]]
 [[5]]
 [[6]]
 [[7]]
 [[3], [4]]
 [[6], [7]]
 [[5], [6, 7]]
 [[3, 4], [5, 6, 7]]
 [[1], [2], [3, 4, 5, 6, 7]]

julia> for n in 8:12
         println("clades below node ", n, ": ", join(below[n], " "))
       end
clades below node 8: [3] [4]
clades below node 9: [6] [7]
clades below node 10: [5] [6, 7]
clades below node 11: [3, 4] [5, 6, 7]
clades below node 12: [1] [2] [3, 4, 5, 6, 7]

julia> above[8:12] # clades sister to and above nodes 8 through 12:
5-element Array{Array{Array{Int64,1},1},1}:
 [[5, 6, 7], [1], [2]]
 [[5], [3, 4], [1], [2]]
 [[3, 4], [1], [2]]
 [[1], [2]]
 []
```
"""
function ladderpartition(net::HybridNetwork)
    nnodes = length(net.node)
    nleaf  = length(net.leaf)
    for n in net.node
        n.number > 0 && n.number <= nnodes || error("node numbers must be in 1 - #nodes")
    end
    sort!([n.number for n in net.leaf]) == collect(1:nleaf) || error("leaves must come first")
    below = Vector{Vector{Vector{Int}}}(undef, nnodes)
    above = Vector{Vector{Vector{Int}}}(undef, nnodes)
    for nni in nnodes:-1:1 # post-order (not pre-order) traversal
        nn = net.nodes_changed[nni]
        above[nn.number] = Vector{Vector{Int}}(undef,0) # initialize
        if nn.leaf
            below[nn.number] = [[nn.number]]
            continue
        end
        below[nn.number]  = Vector{Vector{Int}}(undef,0) # initialize
        !nn.hybrid || error("ladder partitions not implemented for non-tree networks")
        children = [getChild(e) for e in nn.edge]
        filter!(n -> n!=nn, children)
        for cc in children
            allbelowc = union(below[cc.number]...)
            push!(below[nn.number], allbelowc)
            for cc2 in children
                cc2 !=cc || continue
                push!(above[cc2.number], allbelowc)
            end
        end
    end
    # so far: above[n] contains the clades *sister* to node number n, only.
    #         add those above = any above n's parent, using pre-order this time.
    for nni in 2:nnodes # avoid 1: it's the root, no parent, nothing to update
        nn = net.nodes_changed[nni]
        pn = getMajorParent(nn).number # parent number
        for clade in above[pn]  push!(above[nn.number], clade); end
    end
    return below,above
end

"""
    deleteHybridThreshold!(net::HybridNetwork, threshold::Float64,
                           nofuse=false, unroot=false, multgammas=false,
                           keeporiginalroot=false)

Deletes from a network all hybrid edges with heritability below a threshold gamma.
Returns the network.

- if threshold<0.5: delete minor hybrid edges with γ < threshold
  (or with a missing γ, for any threshold > -1.0)
- if threshold=0.5: delete all minor hybrid edges (i.e normally with γ < 0.5, if γ non-missing)
- `nofuse`: if true, do not fuse edges and keep original nodes.
- `unroot`: if false, the root will not be deleted if it becomes of degree 2.
- `multgammas`: if true, the modified edges have γ values equal to the
  proportion of genes that the extracted subnetwork represents. For an edge `e`
  in the modified network, the inheritance γ for `e` is the product of γs
  of all edges in the original network that have been merged into `e`.
-`keeporiginalroot`: if true, the root will be retained even if of degree 1.

Warnings:

- by default, `nofuse` is false, partner hybrid edges are fused with their child edge
  and have their γ changed to 1.0.
  If `nofuse` is true: the γ's of partner hybrid edges are unchanged.
- assumes correct `isMajor` fields, and correct `isChild1` fields to update `containRoot`.
"""
function deleteHybridThreshold!(net::HybridNetwork, gamma::Float64,
                                nofuse=false::Bool, unroot=false::Bool,
                                multgammas=false::Bool, keeporiginalroot=false::Bool)
    gamma <= 0.5 || error("deleteHybridThreshold! called with gamma = $(gamma)>0.5")
    for i = net.numHybrids:-1:1
    # starting from last because net.hybrid changes as hybrids are removed. Empty range if 0 hybrids.
        e = getMinorParentEdge(net.hybrid[i])
        # remove minor edge e if γ < threshold OR threshold=0.5
        # warning: no check if γ and isMajor are in conflict
        if e.gamma < gamma || gamma == 0.5 # note: γ=-1 if missing, so < gamma threshold
            # deleteHybrid!(net.hybrid[i],net,true,false) # requires non-missing edge lengths
            # deleteHybridizationUpdate! requires level-1 network with corresponding attributes
            deletehybridedge!(net, e, nofuse, unroot, multgammas, true, keeporiginalroot)
        end
    end
    return net
end

"""
    displayedNetworks!(net::HybridNetwork, node::Node, keepNode=false,
                       unroot=false, multgammas=false, keeporiginalroot=false)

Extracts the two networks that simplify a given network at a given hybrid node:
deleting either one or the other parent hybrid edge.
If `nofuse` is true, the original edges (and nodes) are kept in both networks,
provided that they have one or more descendant leaves.
If `unroot` is true, the root will be deleted if it becomes of degree 2.
If `keeporiginalroot` is true, the root is retained even if it is of degree 1.

- the original network is modified: the minor edge removed.
- returns one HybridNetwork object: the network with the major edge removed
"""
function displayedNetworks!(net::HybridNetwork, node::Node,
                            nofuse=false::Bool, unroot=false::Bool,
                            multgammas=false::Bool, keeporiginalroot=false::Bool)
    node.hybrid || error("will not extract networks from tree node $(node.number)")
    ind = findfirst(x -> x===node, net.node)
    ind !== nothing || error("node $(node.number) was not found in net")
    netmin = deepcopy(net)
    emin = getMinorParentEdge(node)
    deletehybridedge!(net   , emin, nofuse, unroot, multgammas, true, keeporiginalroot)
    emaj = getMajorParentEdge(netmin.node[ind]) # hybrid node & edge in netmin
    deletehybridedge!(netmin, emaj, nofuse, unroot, multgammas, true, keeporiginalroot)
    return netmin
end

"""
    displayedTrees(net::HybridNetwork, gamma::Float64; nofuse=false::Bool,
                   unroot=false::Bool, multgammas=false::Bool,
                   keeporiginalroot=false::Bool)

Extracts all trees displayed in a network, following hybrid edges
with heritability >= γ threshold (or >0.5 if threshold=0.5)
and ignoring any hybrid edge with heritability lower than γ.
Returns an array of trees, as HybridNetwork objects.

`nofuse`: if true, do not fuse edges (keep degree-2 nodes) during hybrid edge removal.
`unroot`: if false, the root will not be deleted if it becomes of degree 2 unless
  keeporiginalroot is true.
`multgammas`: if true, the edges in the displayed trees have γ values
  equal to the proportion of genes that the edge represents, even though all
  these edges are tree edges. The product of all the γ values across all edges
  is the proportion of genes that the tree represents. More specifically,
  edge `e` in a given displayed tree has γ equal to the product of γs
  of all edges in the original network that have been merged into `e`.
`keeporiginalroot`: if true, keep root even if of degree 1.

Warnings:

- if `nofuse` is true: the retained partner hybrid edges have
  their γ values unchanged, but their `isMajor` is changed to true
- assume correct `isMajor` attributes.
"""
function displayedTrees(net0::HybridNetwork, gamma::Float64;
                        nofuse=false::Bool, unroot=false::Bool,
                        multgammas=false::Bool, keeporiginalroot=false::Bool)
    trees = HybridNetwork[]
    net = deepcopy(net0)
    deleteHybridThreshold!(net,gamma,nofuse,unroot, multgammas, keeporiginalroot)
    displayedTrees!(trees,net,nofuse,unroot, multgammas, keeporiginalroot)
    return trees # should have length 2^net.numHybrids
end

"""
    inheritanceWeight(tree::HybridNetwork)

Return the *log* inheritance weight of a network or tree
(as provided by [`displayedTrees`](@ref) with `nofuse` = true for instance).
For a tree displayed in a network, its inheritance weight is the log of the product
of γ's of all edges retained in the tree. To avoid underflow, the log is calculated:
i.e. sum of log(γ) across retained edges.

If any edge has a negative γ, it is assumed to mean that its γ is missing,
and the function returns `missing`.

# Example

```julia-repl
julia> net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);");

julia> trees = displayedTrees(net,0.0; nofuse=true);

julia> PhyloNetworks.inheritanceWeight.(trees)
2-element Array{Float64,1}:
 -0.105361
 -2.30259
```
"""
function inheritanceWeight(tree::HybridNetwork)
    ltw = 0.0
    for e in tree.edge
        e.gamma != 1.0 || continue
        if e.gamma < 0.0 return missing; end # any negative γ assumed -1.0: missing
        ltw += log(e.gamma)
    end
    return ltw
end

"""
    majorTree(net::HybridNetwork; nofuse=false::Bool, unroot=false::Bool,
              keeporiginalroot=false::Bool)

Extract the major tree displayed in a network, keeping the major edge
and dropping the minor edge at each hybrid node.

`nofuse`: if true, edges and degree-2 nodes are retained during edge removal.
Otherwise, at each reticulation the child edge (below the hybrid node)
is retained: the major hybrid edge is fused with it.

`unroot`: is true, the root will be deleted if it becomes of degree 2.

`keeporiginalroot`: the network's root is kept even if it becomes of degree 1.

Warnings:

- if `nofuse` is true: the hybrid edges that are retained (without fusing)
  have their γ values unchanged, but their `isMajor` is changed to true
- assume correct `isMajor` attributes.
"""
majorTree(net::HybridNetwork; nofuse=false::Bool, unroot=false::Bool,
          keeporiginalroot=false::Bool) =
    displayedTrees(net,0.5; nofuse=nofuse, unroot=unroot,
                   keeporiginalroot=keeporiginalroot)[1]

# expands current list of trees, with trees displayed in a given network
function displayedTrees!(trees::Array{HybridNetwork,1}, net::HybridNetwork,
                        nofuse=false::Bool, unroot=false::Bool,
                        multgammas=false::Bool, keeporiginalroot=false::Bool)
    if isTree(net)
        # warning: no update of edges' containRoot (true) or edges' and nodes' inCycle (-1)
        push!(trees, net)
    else
        netmin = displayedNetworks!(net, net.hybrid[1], nofuse, unroot, multgammas, keeporiginalroot)
        displayedTrees!(trees, net, nofuse, unroot, multgammas, keeporiginalroot)
        displayedTrees!(trees, netmin, nofuse, unroot, multgammas, keeporiginalroot)
    end
end

"""
    minorTreeAt(net::HybridNetwork, hybindex::Integer, nofuse=false, unroot=false::Bool)

Extract the tree displayed in the network, following the major hybrid edge
at each hybrid node, except at the ith hybrid node (i=`hybindex`),
where the minor hybrid edge is kept instead of the major hybrid edge.
If `nofuse` is true, edges are not fused (degree-2 nodes are kept).
If `unroot` is true, the root will be deleted if it becomes of degree 2.

Warning: assume correct `isMajor` fields.
"""
function minorTreeAt(net::HybridNetwork, hybindex::Integer,
            nofuse=false::Bool, unroot=false::Bool)
    hybindex <= length(net.hybrid) || error("network has fewer hybrid nodes than index $(hybindex).")
    tree = deepcopy(net)
    hybedge = getMajorParentEdge(tree.hybrid[hybindex])
    deletehybridedge!(tree, hybedge, nofuse, unroot) # delete major hybrid edge at reticulation of interest
    return majorTree(tree; nofuse=nofuse, unroot=unroot) # all remaining minor edges removed: now it's a tree.
end

"""
    displayedNetworkAt!(net::HybridNetwork, node::Node, nofuse=false,
                        unroot=false, multgammas=false)

Delete all the minor hybrid edges, except at input node. The network is left
with a single hybridization, and otherwise displays the same major tree as before.
If `nofuse` is true, edges are not fused (degree-2 nodes are kept).

Warning: assume correct `isMajor` fields.
"""
function displayedNetworkAt!(net::HybridNetwork, node::Node,
                             nofuse=false::Bool, unroot=false::Bool,
                             multgammas=false::Bool)
    node.hybrid || error("will not extract network from tree node $(node.number)")
    for i = net.numHybrids:-1:1
    # starting from last because net.hybrid changes as hybrids are removed. Empty range if 0 hybrids.
        net.hybrid[i] != node || continue
        emin = getMinorParentEdge(net.hybrid[i])
        deletehybridedge!(net, emin, nofuse, unroot, multgammas)
    end
end

"""
    hardwiredClusterDistance(net1::HybridNetwork, net2::HybridNetwork, rooted::Bool)

Takes 2 networks and returns their hardwired cluster distance, that is,
the number of hardwired clusters found in one network and not in the other.
Note that this is not a distance per se on the full space of hybrid networks:
there are pairs of different networks for which this measure is 0.
But it is a distance on some network subspaces.

If the 2 networks are trees, this is the Robinson-Foulds distance.
If rooted=false, the trees are considered unrooted.

If rooted is false and one of the phylogenies is not a tree (1+ reticulations),
then all degree-2 nodes are removed before comparing the hardwired clusters,
and the minimum distance is returned over all possible ways to root the
networks at internal nodes.
"""
function hardwiredClusterDistance(net1::HybridNetwork, net2::HybridNetwork, rooted::Bool)
    bothtrees = (net1.numHybrids == 0 && net2.numHybrids == 0)
    rooted || bothtrees ||
        return hardwiredClusterDistance_unrooted(net1, net2) # tries all roots, but removes degree-2 nodes
    taxa = sort!(String[net1.leaf[i].name for i in 1:net1.numTaxa])
    length(setdiff(taxa, String[net2.leaf[i].name for i in 1:net2.numTaxa])) == 0 ||
        error("net1 and net2 do not share the same taxon set. Please prune networks first.")
    nTax = length(taxa)
    if bothtrees # even if rooted, different treatment at the root if root=leaf
        M1 = tree2Matrix(net1, taxa, rooted=rooted)
        M2 = tree2Matrix(net2, taxa, rooted=rooted)
    else
        M1 = hardwiredClusters(net1, taxa) # last row: 10/11 if tree/hybrid edge.
        M2 = hardwiredClusters(net2, taxa)
        #println("M1="); print(M1); println("\nM2="); print(M2); println("\n");
    end
    dis = 0
    n2ci = collect(1:size(M2, 1)) # cluster indices
    for i1 in 1:size(M1,1)
        found = false
        m1 = 1 .- M1[i1,2:end] # going to the end: i.e. we want to match a tree edge with a tree edge
                                # and hybrid edge with hybrid edge
        for j in length(n2ci):-1:1 # check only unmatched cluster indices, in reverse
            i2 = n2ci[j]
            if (M1[i1,2:end] == M2[i2,2:end] ||
                  ( !rooted && m1 == M2[i2,2:end])     )
                found = true
                deleteat!(n2ci, j) # a cluster can be repeated
                break
            end
        end
        if !found
            dis += 1
        end
    end # (size(M1)[1] - dis) edges have been found in net2, dis edges have not.
    # so size(M2)[1] - (size(M1)[1] - dis) edges in net2 are not in net1.
    dis + dis + size(M2)[1] - size(M1)[1]
end

"""
    hardwiredClusterDistance_unrooted(net1::HybridNetwork, net2::HybridNetwork)

Miminum hardwired cluster dissimilarity between the two networks, considered as
unrooted (or semi-directed). This dissimilarity is defined as the minimum
rooted distance, over all root positions that are compatible with the direction
of hybrid edges.
Called by [`hardwiredClusterDistance`](@ref).

To avoid repeating identical clusters, all degree-2 nodes
are deleted before starting the comparison.
Since rooting the network at a leaf creates a root node of degree 2 and
an extra cluster, leaves are excluded from possible rooting positions.
"""
function hardwiredClusterDistance_unrooted(net1::HybridNetwork, net2::HybridNetwork)
    return hardwiredClusterDistance_unrooted!(deepcopy(net1), deepcopy(net2))
end
function hardwiredClusterDistance_unrooted!(net1::HybridNetwork, net2::HybridNetwork)
    #= fixit: inefficient function, because r1 * r2 "M" matrices of
      hardwiredClusters() are calculated, where ri = # root positions in neti.
      Rewrite to calculate only r1 + r2 M's.
    =#
    removedegree2nodes!(net1) # because re-rooting would remove them in an
    removedegree2nodes!(net2) # unpredictable order
    # find all permissible positions for the root
    net1roots = [n.number for n in net1.node if !n.leaf]
    #= disallow the root at a leaf: adding a degree-2 node adds a cluster
       that could be artificially matched to a cluster from a degree-3 node
       sister to a hybrid edge, when a the leaf edge is the donor. =#
    for i in length(net1roots):-1:1 # reverse order, to delete some of them
        try
            rootatnode!(net1, net1roots[i]; verbose=false)
            # tricky: rootatnode adds a degree-2 node if i is a leaf,
            #         and delete former root node if it's of degree 2.
        catch e
            isa(e, RootMismatch) || rethrow(e)
            deleteat!(net1roots, i)
        end
    end
    net2roots = [n.number for n in net2.node if !n.leaf]
    for i in length(net2roots):-1:1
        try
            rootatnode!(net2, net2roots[i]; verbose=false)
        catch e
            isa(e, RootMismatch) || rethrow(e)
            deleteat!(net2roots, i)
        end
    end
    bestdissimilarity = typemax(Int)
    bestns = missing
    for n1 in net1roots
        rootatnode!(net1, n1; verbose=false)
        for n2 in net2roots
            rootatnode!(net2, n2; verbose=false)
            diss = hardwiredClusterDistance(net1, net2, true) # rooted = true now
            if diss < bestdissimilarity
                bestns = (n1, n2)
                bestdissimilarity = diss
            end
        end
    end
    # @info "best root nodes: $bestns"
    # warning: original roots (and edge directions) NOT restored
    return bestdissimilarity
end


text:
using Documenter
using AiEmpires

push!(LOAD_PATH,"../src/")
makedocs(sitename="AiEmpires.jl Documentation",
         pages = [
            "Index" => "index.md",
            "AiEmpires" => "api.md",
         ],
         format = Documenter.HTML(prettyurls = false)
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/AiEmpires/AiEmpires.jl.git",
    devbranch = "master"
)


text:
function process(::Type{StatPlotsRecipe}, a::Analysis)
    args = Tuple(Iterators.filter(i -> i !== nothing, (a.x, a.y, a.z)))
    t = a.data.table
    grp_cols = a.data.splitby
    cols = select(t, All(args, grp_cols)) |> dropna
    a.plot(columns(cols, args)...; group = columns(cols, grp_cols),  a.plot_kwargs...)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiFlowcontrolV1alpha1FlowSchema <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus } # spec name: status

    function IoK8sApiFlowcontrolV1alpha1FlowSchema(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchema, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1FlowSchema

const _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchema = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiFlowcontrolV1alpha1FlowSchema = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec", Symbol("status")=>"IoK8sApiFlowcontrolV1alpha1FlowSchemaStatus")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1FlowSchema))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1FlowSchema[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchema[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1FlowSchema)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchema }, name::Symbol, val)
end


text:
@testset "Renewable data misspecification" begin
    # See https://discourse.julialang.org/t/how-to-use-test-warn/15557/5 about testing for warning throwing
    warn_message = "The data doesn't include devices of type HydroDispatch, consider changing the device models"
    model = DeviceModel(PSY.HydroDispatch, PSI.HydroDispatchRunOfRiver)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Hydro, model)
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys14)
    @test_logs (:warn, warn_message) construct_device!(op_problem, :Hydro, model)
end

@testset "Hydro DCPLossLess FixedOutput" begin
    model = DeviceModel(PSY.HydroFix, PSI.HydroFixed)

    # Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hy; use_parameters = true)
    construct_device!(op_problem, :Hydro, model)
    moi_tests(op_problem, true, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hy)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hy; use_parameters = true, use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, true, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast - No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hy; use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 0, 0, 0, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

end

@testset "Hydro DCPLossLess HydroDispatch with HydroDispatchRunOfRiver formulations" begin
    model = DeviceModel(PSY.HydroDispatch, PSI.HydroDispatchRunOfRiver)

    # Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true)
    construct_device!(op_problem, :Hydro, model)
    moi_tests(op_problem, true, 24, 0, 24, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 24, 0, 24, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true, use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, true, 1, 0, 1, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast - No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 1, 1, 0, 0, 0, false)
    psi_checkobjfun_test(op_problem, GAEVF)

end

# @testset "Hydro DCPLossLess HydroDispatch with HydroDispatchSeasonalFlow Formulations" begin
#     model = DeviceModel(PSY.HydroDispatch, PSI.HydroDispatchSeasonalFlow)

#     # Parameters Testing
#     op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true)
#     construct_device!(op_problem, :Hydro, model)
#     moi_tests(op_problem, true, 24, 0, 24, 0, 0, false)
#     psi_checkobjfun_test(op_problem, GAEVF)

#     # No Parameters Testing
#     op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd)
#     construct_device!(op_problem, :Hydro, model);
#     moi_tests(op_problem, false, 24, 0, 24, 0, 0, false)
#     psi_checkobjfun_test(op_problem, GAEVF)

#     # No Forecast Testing
#     op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true, use_forecast_data = false)
#     construct_device!(op_problem, :Hydro, model);
#     moi_tests(op_problem, true, 1, 0, 1, 0, 0, false)
#     psi_checkobjfun_test(op_problem, GAEVF)

#     # No Forecast - No Parameters Testing
#     op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_forecast_data = false)
#     construct_device!(op_problem, :Hydro, model);
#     moi_tests(op_problem, false, 1, 1, 0, 0, 0, false)
#     psi_checkobjfun_test(op_problem, GAEVF)

# end

@testset "Hydro DCPLossLess HydroDispatch with HydroCommitmentRunOfRiver Formulations" begin
    model = DeviceModel(PSY.HydroDispatch, PSI.HydroCommitmentRunOfRiver)

    # Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true)
    construct_device!(op_problem, :Hydro, model)
    moi_tests(op_problem, true, 96, 0, 72, 0, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 96, 0, 48, 0, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true, use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, true, 4, 0, 3, 0, 1, true)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast - No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 4, 0, 2, 1, 1, true)
    psi_checkobjfun_test(op_problem, GAEVF)

end

@testset "Hydro DCPLossLess HydroDispatch with HydroCommitmentSeasonalFlow Formulations" begin
    model = DeviceModel(PSY.HydroDispatch, PSI.HydroCommitmentSeasonalFlow)

    # Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true)
    construct_device!(op_problem, :Hydro, model)
    moi_tests(op_problem, true, 96, 0, 72, 0, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 96, 0, 48, 0, 24, true)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_parameters = true, use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, true, 4, 0, 3, 0, 1, true)
    psi_checkobjfun_test(op_problem, GAEVF)

    # No Forecast - No Parameters Testing
    op_problem = OperationsProblem(TestOpProblem, DCPPowerModel, c_sys5_hyd; use_forecast_data = false)
    construct_device!(op_problem, :Hydro, model);
    moi_tests(op_problem, false, 4, 0, 2, 1, 1, true)
    psi_checkobjfun_test(op_problem, GAEVF)

end
#=
@testset " Hydro Tests" begin
    PSI.activepower_variables(ps_model, generators_hg, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, DCPPowerModel, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, ACPPowerModel, 1:24)
    PSI.reactivepower_variables(ps_model, generators_hg, 1:24)
    PSI.reactivepower_constraints(ps_model, generators_hg, PSI.HydroDispatchRunOfRiver, ACPPowerModel, 1:24)
end

@testset " Hydro Tests" begin
    ps_model = PSI._canonical_init(length(sys5b.buses), nothing, PM.AbstractPowerModel, sys5b.time_periods)
    PSI.activepower_variables(ps_model, generators_hg, 1:24)
    PSI.commitment_variables(ps_model, generators_hg, 1:24);
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, DCPPowerModel, 1:24)
    PSI.activepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, ACPPowerModel, 1:24)
    PSI.reactivepower_variables(ps_model, generators_hg, 1:24)
    PSI.reactivepower_constraints(ps_model, generators_hg, PSI.HydroCommitmentRunOfRiver, ACPPowerModel, 1:24)
end
=#


text:
@testset "h2fk5" begin
    rh  =  1.767794352
    dh  = -0.2917512594
    drh = -2.76413026e-6
    ddh = -5.92994449e-6
    pxh =  0.379210
    rvh = -7.6
    r5, d5, dr5, dd5, px5, rv5 = h2fk5(rh, dh, drh, ddh, pxh, rvh)
    @test r5 ≈ 1.767794455700065506 atol=1e-13
    @test d5 ≈ -0.2917513626469638890 atol=1e-13
    @test dr5 ≈ -0.27597945024511204e-5 atol=1e-18
    @test dd5 ≈ -0.59308014093262838e-5 atol=1e-18
    @test px5 ≈ 0.37921 atol=1e-13
    @test rv5 ≈ -7.6000001309071126 atol=1e-10
end

@testset "hd2ae" begin
    h = 1.1
    d = 1.2
    p = 0.3

    a, e = ERFA.hd2ae(h, d, p)

    @test a ≈ 5.916889243730066194 atol=1e-13
    @test e ≈ 0.4472186304990486228 atol=1e-14
end

@testset "hd2pa" begin
    h = 1.1
    d = 1.2
    p = 0.3

    q = ERFA.hd2pa(h, d, p)

    @test q ≈ 1.906227428001995580 atol=1e-13
end

@testset "hfk5z" begin
    rh =  1.767794352
    dh = -0.2917512594
    r5, d5, dr5, dd5 = hfk5z(rh, dh, 2400000.5, 54479.0)
    @test r5 ≈ 1.767794490535581026 atol=1e-13
    @test d5 ≈ -0.2917513695320114258 atol=1e-14
    @test dr5 ≈ 0.4335890983539243029e-8 atol=1e-22
    @test dd5 ≈ -0.8569648841237745902e-9 atol=1e-23
end