text:
const RExpr = Union{Expr,Symbol,Number}

replace_post(ex, s, s′) = MacroTools.postwalk(x -> x == s ? s′ : x, ex)
replace_post(ex, dict) = MacroTools.postwalk(x -> x in keys(dict) ? dict[x] : x, ex)

function free_symbols(ex::Expr)
    ls = Symbol[]
    MacroTools.postwalk(x -> x isa Symbol && Base.isidentifier(x) ? push!(ls, x) : x, ex)
    Base.unique(ls)
end

free_symbols(s::Symbol) = [s]
free_symbols(x) = []

multiplication_walk(f, expr) = MacroTools.prewalk(expr) do x
    if @capture(x, g_(a__))
        if g == :*
            f(a)
        elseif g == :^
            f(first(a))
        end
    end
    x
end

function find_multiplications(ex::RExpr)
    ex = convert(Expr, expand(convert(Basic, ex)))
    syms = Symbol[]
    multiplication_walk(ex) do vs
        if vs isa Expr
            append!(syms, free_symbols(vs))
        elseif vs isa Vector
            filter!(x->!isa(x, Number), vs)
            if length(vs) > 1
                for v in vs
                    append!(syms, free_symbols(v))
                end
            end
        elseif vs isa Symbol
            push!(syms, vs)
        else
            error("Something else: $vs")
        end
    end
    Base.unique(syms)
end

function split_assign(xs::Vector{Expr})
    ls = Symbol[]
    rs = RExpr[]
    for x in xs
        @capture(x, l_ = r_)
        push!(ls, unblock(l))
        push!(rs, unblock(r))
    end
    ls, rs
end

# function iscounter(l::Symbol, r::RExpr)
#     @match r begin
#         $l + c_ => c isa Number
#         c_ + $l => c isa Number
#         $l - c_ => c isa Number
#         c_ - $l => c isa Number
#         _      => false
#     end
# end

function _parallel(lhss::Vector{Symbol}, rhss::Vector{RExpr})
    del = Int[]
    dict = collect(zip(lhss, rhss))
    for (i,v) in enumerate(lhss)
        rhss[i] = replace_post(rhss[i], Dict(dict[1:i-1]))
    end
    _lhss, _rhss = Symbol[], RExpr[]
    for (l, r) in reverse(collect(zip(lhss, rhss)))
        if l ∉ _lhss
            pushfirst!(_lhss, l)
            pushfirst!(_rhss, r)
        end
    end
    _lhss, _rhss
end

# function lrs_sequential(exprs::Vector{Expr}, lc::Symbol = gensym_unhashed(:n))
#     lhss, rhss = _parallel(split_assign(exprs)...)
#     @debug "Splitted and parallel assignments" rhss lhss

#     lines = Expr[]
#     mvars = Symbol[]
#     loop_counters = Symbol[]
#     for (l, r) in zip(lhss, rhss)
#         _mvars = find_multiplications(r)
#         filter!(x->x!=l, _mvars)
#         append!(mvars, _mvars)
#         if iscounter(l, r)
#             push!(loop_counters, l)
#         end
#         _l = Recurrences.symbol_walk(x->:($x($(lc)+1)), l)
#         _r = Recurrences.symbol_walk(x->:($x($(lc))), r)
#         push!(lines, :($(_l) = $(_r)))
#     end
#     Recurrences.lrs(lines)
# end

function __transform_body(expr::Expr)
    expr = MacroTools.postwalk(x->x isa Expr && x.head == :elseif ? Expr(:if, x.args...) : x, expr)
    branches = _transform(MacroTools.striplines(expr))

    vars = Base.unique(Iterators.flatten(map(free_symbols, branches)))
    for b in branches
        lhss, _ = split_assign(Vector{Expr}(b.args))
        left = setdiff(vars, lhss)
        for v in left
            push!(b.args, :($v = $v))
        end
    end
    branches
end

merge(x::Expr, ys::Vector{Expr}) = [merge(x, y) for y in ys]
merge(xs::Vector{Expr}, y::Expr) = [merge(x, y) for x in xs]
merge(xs::Vector{Expr}, ys::Vector{Expr}) = [merge(x, y) for x in xs for y in ys]
merge(x) = x

function merge(x::Expr, y::Expr)
    @assert x.head == :block && y.head == :block
    Expr(:block, x.args..., y.args...)
end

function _transform(expr)
    res = @match expr begin
        if _ x_ end         => [_transform(x); Expr(:block)]
        if _ x_ else y_ end => [_transform(x); _transform(y)]
        (x_ = y_)           => [block(expr)]
    end
    if res === nothing
        if expr isa Expr && expr.head == :block
            res = merge(reduce(merge, _transform.(expr.args)))
        else
            @warn("Ignoring unsupported program construct: $expr")
            res = Expr(:block)
        end
    end
    return res
end

function __extract_init(xs::Vector{Expr})
    m = Dict{Symbol, RExpr}()
    for x in xs
        if @capture(x, l_tuple = r_tuple)
            for (ll, rr) in zip(l.args, r.args)
                m[ll] = rr
            end
        elseif @capture(x, l_ = r_)
            m[l] = r
        end
    end
    return m
end

function transform(expr::Expr)
    if @capture(expr, (init__; while _ body_ end))
        __extract_init(Vector{Expr}(init)), __transform_body(body)
    end
end


text:
Random.seed!(209)

@testset "0D grid" begin
  @test Agents.grid0D() == Agents.Graph(1)
  obj1 = grid()
  obj2 = Agents.grid0D()
  @test Agents.nv(obj1) == Agents.nv(obj2)
  @test Agents.ne(obj1) == Agents.ne(obj2)
  # @test gridsize(obj1) == 1
end

@testset "1D grid" begin
  obj1 = Agents.grid1D(5)
  obj2 = Agents.grid1D(5, periodic=true)
  @test Agents.nv(obj1) == 5
  @test Agents.nv(obj2) == 5
  @test Agents.ne(obj1) == 4
  @test Agents.ne(obj2) == 5
  last_edge1 = collect(Agents.edges(obj1))[end]
  last_edge2 = collect(Agents.edges(obj2))[2]
  @test last_edge1.src == 4
  @test last_edge1.dst == 5
  @test last_edge2.src == 1
  @test last_edge2.dst == 5
end

@testset "2D grid" begin
  obj1 = Agents.grid2D(3, 4)
  obj2 = Agents.grid2D(3, 4, periodic=true)
  @test Agents.nv(obj1) == 12
  @test Agents.nv(obj2) == 12
  @test Agents.ne(obj1) == 17
  @test Agents.ne(obj2) == 24
end

@testset "2D Moore" begin
  obj1 = Agents.grid2D_Moore(3, 3)
  @test Agents.nv(obj1) == 9
  @test Agents.ne(obj1) == 20
  obj1 = Agents.grid2D_Moore(3, 4)
  @test Agents.nv(obj1) == 12
  @test Agents.ne(obj1) == 29
  obj2 = Agents.grid2D_Moore(3, 2, periodic=true)
  @test Agents.nv(obj2) == 6
  @test Agents.ne(obj2) == 15
  obj2 = Agents.grid2D_Moore(3, 3, periodic=true)
  @test Agents.nv(obj2) == 9
  @test Agents.ne(obj2) == 36
end

@testset "3D grid" begin
  g1 = Agents.grid3D(2,3,2)
  g2 = Agents.grid3D(2,3,3)
  g3 = Agents.grid3D(2,3,2, periodic=true)
  g4 = Agents.grid3D(2,3,3, periodic=true)
  @test Agents.ne(g1) == 20
  @test Agents.ne(g2) == 33
  @test Agents.ne(g3) == 24
  @test Agents.ne(g4) == 39
end

@testset "grid coord/vertex conversions" begin
  @test coord_to_vertex((2,2,1), (3,4,1)) == 5
  @test coord_to_vertex((2,2), (3,4)) == 5
  @test coord_to_vertex(2,2,1, (3,4,1)) == 5
  @test coord_to_vertex(2,2, (3,4)) == 5
  @test vertex_to_coord(5, (3,4,1)) == (2,2,1)
  @test vertex_to_coord(5, (3,4)) == (2,2)

  @test coord_to_vertex((2,2,1), (2,3,3)) == 4
  @test coord_to_vertex((2,2,2), (2,3,3)) == 10
  @test coord_to_vertex((2,2,3), (2,3,3)) == 16
  @test coord_to_vertex((1,3,1), (2,3,3)) == 5
  @test coord_to_vertex((1,3,2), (2,3,3)) == 11
  @test coord_to_vertex((1,3,3), (2,3,3)) == 17

  @test vertex_to_coord(5, (2,3,3)) == (1,3,1)
  @test vertex_to_coord(7, (2,3,3)) == (1,1,2)
  @test vertex_to_coord(13, (2,3,3)) == (1,1,3)
  @test vertex_to_coord(13, (2,3,3)) == (1,1,3)
  @test vertex_to_coord(15, (2,3,3)) == (1,2,3)
  @test vertex_to_coord(18, (2,3,3)) == (2,3,3)
  @test vertex_to_coord(18, (2,3,3)) == (2,3,3)
end

@testset "grid" begin
  @test Agents.grid(1,1,1, false, true) == Agents.grid0D()
  @test Agents.grid(1,1,1, true, true) == Agents.grid0D()
  @test Agents.grid(6,1,1, true, true) == Agents.grid1D(6, periodic=true)
  @test Agents.grid(6,1,1, false, true) == Agents.grid1D(6, periodic=false)
  @test Agents.grid(6,5,1, false, true) == Agents.grid2D_Moore(6, 5, periodic=false)
  @test Agents.grid(6,5, false, true) == Agents.grid2D_Moore(6, 5, periodic=false)
  @test Agents.grid(6,5,1, true, true) == Agents.grid2D_Moore(6, 5, periodic=true)
  @test Agents.grid(6,5, true, true) == Agents.grid2D_Moore(6, 5, periodic=true)
  @test Agents.grid(6,5,1, true, false) == Agents.grid2D(6, 5, periodic=true)
  @test Agents.grid(6,5, true, false) == Agents.grid2D(6, 5, periodic=true)
  @test Agents.grid(6,5,1, false, false) == Agents.grid2D(6, 5, periodic=false)
  @test Agents.grid(6,5, false, false) == Agents.grid2D(6, 5, periodic=false)
  @test Agents.grid((3,2,1), false, true) == Agents.grid(3,2,1, false, true)
  @test Agents.grid((3,2), false, true) == Agents.grid(3,2, false, true)
end

@testset "gridsize" begin
  @test Agents.gridsize((3,4,6)) == 3*4*6
  @test Agents.gridsize((3,4)) == 12
  @test Agents.gridsize((3,4)) == Agents.gridsize((3,4,1))
end

@testset "all the rest" begin

  model = model_initiation(f=0.1, d=0.8, p=0.1, griddims=(20, 20), seed=2)

  agent = model.agents[1]
  move_agent!(agent, (3,4), model)  # node number 63
  @test agent.pos == (3,4)
  @test agent.id in model.space.agent_positions[63]

  agent = model.agents[2]
  move_agent!(agent, 83, model)  # pos (3,5)
  @test agent.pos == (3,5)
  @test agent.id in model.space.agent_positions[83]

  new_pos = move_agent!(agent, model)
  @test agent.id in model.space.agent_positions[coord_to_vertex(new_pos, model)]

  add_agent!(agent, (2,9), model)
  @test agent.pos == (2,9)
  @test agent.id in model.space.agent_positions[coord_to_vertex((2,9), model)]
  @test agent.id in model.space.agent_positions[coord_to_vertex(new_pos, model)]

  @test agent.id in get_node_contents(agent, model)

  ii = model.agents[end].id
  @test id_to_agent(ii, model) == model.agents[end]
end

ext: jl


text:
# Autogenerated wrapper script for libwatcher_jll for x86_64-linux-gnu-cxx11
export libwatcher

using LibUV_jll
JLLWrappers.@generate_wrapper_header("libwatcher")
JLLWrappers.@declare_library_product(libwatcher, "libwatcher.so")
function __init__()
    JLLWrappers.@generate_init_header(LibUV_jll)
    JLLWrappers.@init_library_product(
        libwatcher,
        "lib/libwatcher.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "CorrectedMutualInformation" begin

    Gaoetal2011 = joinpath(pwd(), "data", "Gaoetal2011.fasta")

    function gao11_buslje09(measure)
        filename = string("data_Gaoetal2011_soft_Busljeetal2009_measure_", measure, ".txt")
        joinpath(pwd(), "data", filename)
    end

    ## Column numbers for the output of Buslje et. al. 2009
    SCORE = 9
    ZSCORE = 12

    MIToS_SCORE = 2
    MIToS_ZSCORE = 1

    @testset "APC!" begin

        for MI in ([ 0. 2. 4.
                     2. 0. 6.
                     4. 6. 0. ],
                    PairwiseListMatrix([2., 4., 6.]),
                    PairwiseListMatrix([0., 2., 4., 0., 6., 0.], true),
                    NamedArray(PairwiseListMatrix([2., 4., 6.])),
                    NamedArray(PairwiseListMatrix([0., 2., 4., 0., 6., 0.],true)))

            if isa(MI, Matrix{Float64})
                @test Information._mean_column(MI) == [3., 4., 5.]
                @test Information._mean_total(MI) == 4.
            end

            MIp = APC!(MI)

            @test filter(x -> !isnan(x), vec(Matrix(MIp))) ≈ filter(x -> !isnan(x), vec(
                [   NaN -1.0  0.25
                    -1.0  NaN  1.00
                    0.25 1.00   NaN ]))
        end
    end

    @testset "Simple example I" begin

        aln = Residue[ 'A' 'A'
                       'A' 'R' ]

        @testset "MI" begin

            # Fill the Pij matrix of the example
            Pij = zeros(Float64, 20, 20);
            N = 2 # There are 2 sequences
            Pij[1,1] = (1/N) # A A
            Pij[1,2] = (1/N) # A R
            @test sum(Pij) ≈ 1.0 # Is the matrix correct?
            # Fill the marginals
            Pi = dropdims(Base.sum(Pij, dims=2), dims=2)
            @test Pi[1] == 1.0 # Is always A
            Pj = dropdims(Base.sum(Pij, dims=1), dims=1)
            @test Pj[1] == 0.5 # A
            @test Pj[2] == 0.5 # R
            # Start to sum with 0.0
            total = 0.0
            for i in 1:20, j in 1:20
                if Pij[i,j] != 0.0 && Pi[i] != 0.0 && Pj[j] != 0.0
                    total += (Pij[i,j] * log(Pij[i,j]/(Pi[i]*Pj[j]))) # 0.5 * log(0.5/(0.5 * 1.0)) == 0.0
                end
            end
            # Compare with MIToS result
            mi = mapcolpairfreq!(mutual_information, aln,
                    Probabilities(ContingencyTable(Float64,Val{2},UngappedAlphabet())))
            @test mi[1,2] == total
        end

        @testset "MI: Using pseudocount (0.05)" begin

            Pij = zeros(Float64, 20, 20);
            N = (400 * 0.05) + 2 # 0.05 is the pseudocount and there are 2 sequences
            fill!(Pij, 0.05/N);
            Pij[1,1] =(1.05/N) # A A
            Pij[1,2] =(1.05/N) # A R
            @test sum(Pij) ≈ 1.0
            Pi = dropdims(Base.sum(Pij, dims=2), dims=2)
            Pj = dropdims(Base.sum(Pij, dims=1), dims=1)
            total = 0.0
            for i in 1:20, j in 1:20
                total += (Pij[i,j] * log(Pij[i,j]/(Pi[i]*Pj[j])))
            end
            # Compare with MIToS result
            mi = mapcolpairfreq!(mutual_information, aln,
                    Probabilities(ContingencyTable(Float64,Val{2},UngappedAlphabet())),
                    pseudocounts=AdditiveSmoothing(0.05))
            @test mi[1,2] ≈ total

            @testset "APC!" begin

                APC!(mi)
                @test isnan(mi[1,1])
                @test isapprox(mi[1,2], 0.0, rtol=1e-16)

                zerodiagonal = Float64[ 0.0 total
                                        total 0.0 ]

                # MI.. == MI.j == MIi. == total
                # APC = (total * total) / total == total
                # MI - APC == total - total == 0.0
                APC!(zerodiagonal)
                @test zerodiagonal[1,2] ≈ 0.0
            end
        end

        @testset "Z-score" begin

            # 2 Possibilities:  A A   A A
            #                   A R   R A
            # Is almost the same MI, Z score should be 0.0
            results = buslje09(aln, lambda=0.05, clustering=false, apc=false)

            @test results[MIToS_ZSCORE][1,2] ≈ 0.0

            @test aln == Residue[ 'A' 'A'
                                  'A' 'R' ]
        end
    end

    @testset "Simple example II" begin

        aln = Residue[ 'R' 'A'
                       'A' 'R' ]

        @testset "MI" begin

            Pij = zeros(Float64, 20, 20);
            N = 2 # There are 2 sequences
            Pij[2,1] = (1/N) # R A
            Pij[1,2] = (1/N) # A R
            @test sum(Pij) ≈ 1.0
            Pi = dropdims(Base.sum(Pij, dims=2), dims=2)
            @test Pi[2] == 0.5 # R
            @test Pi[1] == 0.5 # A
            Pj = dropdims(Base.sum(Pij, dims=1), dims=1)
            @test Pj[1] == 0.5 # A
            @test Pj[2] == 0.5 # R
            total = 0.0
            for i in 1:20, j in 1:20
                if Pij[i,j] != 0.0 && Pi[i] != 0.0 && Pj[j] != 0.0
                    total += (Pij[i,j] * log(Pij[i,j]/(Pi[i]*Pj[j])))
                end
            end
            @test total == 0.5 * log(2) + 0.5 * log(2) # 0.5/(0.5*0.5) == 2
            # Compare with MIToS result
            mi = mapcolpairfreq!(mutual_information, aln,
                    Probabilities(ContingencyTable(Float64,Val{2},UngappedAlphabet())))
            @test mi[1,2] == total
        end

        @testset "MI: Using pseudocount (0.05)" begin

            Pij = zeros(Float64, 20, 20);
            N = (400 * 0.05) + 2 # 0.05 is the pseudocount and there are 2 sequences
            fill!(Pij, 0.05/N);
            Pij[2,1] =(1.05/N) # R A
            Pij[1,2] =(1.05/N) # A R
            @test sum(Pij) ≈ 1.0
            Pi = dropdims(Base.sum(Pij, dims=2), dims=2)
            Pj = dropdims(Base.sum(Pij, dims=1), dims=1)
            total = 0.0
            for i in 1:20, j in 1:20
                total += (Pij[i,j] * log(Pij[i,j]/(Pi[i]*Pj[j])))
            end
            # Compare with MIToS result
            mi = mapcolpairfreq!(mutual_information, aln,
                    Probabilities(ContingencyTable(Float64,Val{2},UngappedAlphabet())),
                    pseudocounts=AdditiveSmoothing(0.05))
            @test mi[1,2] ≈ total

            @testset "APC!" begin

                APC!(mi)
                @test isnan(mi[1,1])
                @test mi[1,2] ≈ 0.0

                zerodiagonal = Float64[ 0.0 total
                                        total 0.0 ]

                # MI.. == MI.j == MIi. == total
                # APC = (total * total) / total == total
                # MI - APC == total - total == 0.0
                APC!(zerodiagonal)
                @test zerodiagonal[1,2] ≈ 0.0
            end
        end

        @testset "Z-score" begin

            # 4 Possibilities:  R A   R A   A R   A R
            #                   A R   R A   A R   R A
            mi = mapcolpairfreq!(mutual_information, aln,
                    Probabilities(ContingencyTable(Float64,Val{2},UngappedAlphabet())),
                    pseudocounts=AdditiveSmoothing(0.05))
            other_posib = Residue[ 'A' 'R'
                                   'A' 'R' ]
            other_mi = mapcolpairfreq!(mutual_information, other_posib,
                        Probabilities(ContingencyTable(Float64,Val{2},UngappedAlphabet())),
                        pseudocounts=AdditiveSmoothing(0.05))
            # The mean should be:
            r_mean = 0.5 * ( mi[1,2] + other_mi[1,2] )
            # And the std should be similar to:
            r_std = 0.5 * ( sqrt((mi[1,2]-r_mean)^2) + sqrt((other_mi[1,2]-r_mean)^2) )

            results = buslje09(aln, lambda=0.05, clustering=false, apc=false, samples=100)
            @test isapprox(results[MIToS_ZSCORE][1,2], (mi[1,2]-r_mean)/r_std, atol=0.55)

            results = buslje09(aln, lambda=0.05, clustering=false, apc=false, samples=1000)
            @test isapprox(results[MIToS_ZSCORE][1,2], (mi[1,2]-r_mean)/r_std, atol=0.15)

            results = buslje09(aln, lambda=0.05, clustering=false, apc=false, samples=10000)
            @test isapprox(results[MIToS_ZSCORE][1,2], (mi[1,2]-r_mean)/r_std, atol=0.055)

            @test aln == Residue[ 'R' 'A'
                                  'A' 'R' ]
        end
    end

    @testset "Results from Buslje et. al. 2009" begin

        @testset "Simple" begin
            data = readdlm(joinpath(pwd(), "data",
                           "data_simple_soft_Busljeetal2009_measure_MI.txt"), comments=true)
            results = buslje09(joinpath(pwd(), "data", "simple.fasta"), FASTA,
                               lambda=0.0, clustering=false, apc=false)

            @test isapprox(Float64(data[1, SCORE]), results[MIToS_SCORE][1,2], atol=1e-6)
            @test isapprox(Float64(data[1, ZSCORE]), results[MIToS_ZSCORE][1,2], atol=2.)
        end

        @testset "Gaoetal2011" begin
            data = readdlm(gao11_buslje09("MI"), comments=true)
            results = buslje09(Gaoetal2011, FASTA, lambda=0.0, clustering=false, apc=false)

            @test isapprox([ Float64(x) for x in data[:,SCORE] ],
                           matrix2list(results[MIToS_SCORE]), atol=1e-6)
            @test isapprox([ Float64(x) for x in data[:,ZSCORE] ],
                           matrix2list(results[MIToS_ZSCORE]), atol=2.)

            # println(cor(convert(Vector{Float64},data[:,ZSCORE]),matrix2list(results[MIToS_ZSCORE])))

            @testset "cMI" begin
                result = copy(results[1])
                result[diagind(result)] .= 0.0
                @test cumulative(results[1], -Inf) == sum(result, dims=1)
                @test all((cumulative(results[2],0.5) .-
                           [0.0,0.0,1.38629,1.38629,1.38629,0.0]') .< 0.00001)
                           #        1.38629 = 0.693147 + 0.693147
            end

            result_0_05 = buslje09(Gaoetal2011,FASTA,lambda=0.05,clustering=false,apc=false)
            @test isapprox(result_0_05[MIToS_SCORE][1,2], 0.33051006116310444, atol=1e-14)
        end
    end

    @testset "MI + clustering" begin

        data = readdlm(gao11_buslje09("MI_clustering"), comments=true)
        results = buslje09(Gaoetal2011, FASTA, lambda=0.0, clustering=true, apc=false)

        @test isapprox([ Float64(x) for x in data[:,SCORE] ],
                       matrix2list(results[MIToS_SCORE]), atol=1e-6)
        @test isapprox([ Float64(x) for x in data[:,ZSCORE] ],
                       matrix2list(results[MIToS_ZSCORE]), atol=2.)
    end

    @testset "MIp" begin

        data = readdlm(gao11_buslje09("MI_APC"), comments=true)
        results = buslje09(Gaoetal2011, FASTA, lambda=0.0, clustering=false, apc=true)

        @test isapprox([ Float64(x) for x in data[:,SCORE] ],
                       matrix2list(results[MIToS_SCORE]), atol=1e-6)
        @test isapprox([ Float64(x) for x in data[:,ZSCORE] ],
                       matrix2list(results[MIToS_ZSCORE]), atol=2.0)

        @test isapprox(results[MIToS_SCORE][5,6], 0.018484, atol=0.000001)
    end

    @testset "MIp + clustering" begin

        data = readdlm(gao11_buslje09("MI_APC_clustering"), comments=true)
        results = buslje09(Gaoetal2011, FASTA, lambda=0.0, clustering=true, apc=true)

        @test isapprox([ Float64(x) for x in data[:,SCORE] ],
                       matrix2list(results[MIToS_SCORE]), atol=1e-6)
        @test isapprox([ Float64(x) for x in data[:,ZSCORE] ],
                       matrix2list(results[MIToS_ZSCORE]), atol=2.)

        @test isapprox(results[MIToS_SCORE][5,6], 0.018484, atol=0.000001)
    end

    @testset "BLMI" begin

        @testset "Simple" begin
            file = joinpath(pwd(), "data", "simple.fasta")
            busl = buslje09(file, FASTA)
            blmi = BLMI(file, FASTA)

            @test PairwiseListMatrices.getlist(busl[1]) ≈
                  PairwiseListMatrices.getlist(blmi[1])
            @test PairwiseListMatrices.getlist(busl[2]) ≈
                  PairwiseListMatrices.getlist(blmi[2])
        end

        @testset "Gaoetal2011" begin
            msa  = read(Gaoetal2011, FASTA)
            busl = buslje09(Gaoetal2011, FASTA, lambda=0.0, samples=0)
            blmi = BLMI(msa, lambda=0.0, beta=0.0, samples=5)
            # BLMI should be equal to Buslje09 if beta is zero

            @test PairwiseListMatrices.getlist(busl[2]) ≈
                  PairwiseListMatrices.getlist(blmi[2]) # MIapc
            @test msa == read(Gaoetal2011, FASTA)
        end

        @testset "Gaoetal2011, lambda 0.05" begin

            busl = buslje09(Gaoetal2011, FASTA, lambda=0.5, samples=0)
            blmi = BLMI(Gaoetal2011, FASTA, lambda=0.5, beta=0.0, samples=5)
            # BLMI should be equal to Buslje09 if beta is zero

            @test PairwiseListMatrices.getlist(busl[2]) ≈
                  PairwiseListMatrices.getlist(blmi[2]) # MIapc
        end
    end
end


text:
@testset "4.2.2.2 (g sin)^p (a+b cos)^m (c+d cos)^n" begin
    (a, b, d, e, f, g, x, ) = @variables a b d e f g x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^p*(a+a*cos(e+f*x))^m*(c-c*cos(e+f*x))^n=#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^p*(a+a*cos(e+f*x))^m*(c+d*cos(e+f*x))^n=#

    #= ::Title:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^p*(a+b*cos(e+f*x))^m*(c+d*cos(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^m*(a+b*cos(e+f*x))^n / (c+d*cos(e+f*x))=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(g*sin(e+f*x))^(m/2)*(d*cos(e+f*x))^(n/2) / (a+b*cos(e+f*x))=#

    @test_int [sqrt(g*sin(e + f*x))*sqrt(d*cos(e + f*x))/(a + b*cos(e + f*x)), x, 16, -((sqrt(d)*sqrt(g)*atan(1 - (sqrt(2)*sqrt(d)*sqrt(g*sin(e + f*x)))/(sqrt(g)*sqrt(d*cos(e + f*x)))))/(sqrt(2)*b*f)) + (sqrt(d)*sqrt(g)*atan(1 + (sqrt(2)*sqrt(d)*sqrt(g*sin(e + f*x)))/(sqrt(g)*sqrt(d*cos(e + f*x)))))/(sqrt(2)*b*f) + (2*sqrt(2)*a*d*sqrt(g)*sqrt(cos(e + f*x))*Elliptic.Pi(-(sqrt(-a + b)/sqrt(a + b)), asin(sqrt(g*sin(e + f*x))/(sqrt(g)*sqrt(1 + cos(e + f*x)))), -1))/(b*sqrt(-a + b)*sqrt(a + b)*f*sqrt(d*cos(e + f*x))) - (2*sqrt(2)*a*d*sqrt(g)*sqrt(cos(e + f*x))*Elliptic.Pi(sqrt(-a + b)/sqrt(a + b), asin(sqrt(g*sin(e + f*x))/(sqrt(g)*sqrt(1 + cos(e + f*x)))), -1))/(b*sqrt(-a + b)*sqrt(a + b)*f*sqrt(d*cos(e + f*x))) + (sqrt(d)*sqrt(g)*log(sqrt(g) - (sqrt(2)*sqrt(d)*sqrt(g*sin(e + f*x)))/sqrt(d*cos(e + f*x)) + sqrt(g)*tan(e + f*x)))/(2*sqrt(2)*b*f) - (sqrt(d)*sqrt(g)*log(sqrt(g) + (sqrt(2)*sqrt(d)*sqrt(g*sin(e + f*x)))/sqrt(d*cos(e + f*x)) + sqrt(g)*tan(e + f*x)))/(2*sqrt(2)*b*f)]
    @test_int [sqrt(d*cos(e + f*x))/(sqrt(g*sin(e + f*x))*(a + b*cos(e + f*x))), x, 4, (2*sqrt(2)*sqrt(d)*Elliptic.Pi(-(a/(b - sqrt(-a^2 + b^2))), asin(sqrt(d*cos(e + f*x))/(sqrt(d)*sqrt(1 + sin(e + f*x)))), -1)*sqrt(sin(e + f*x)))/(sqrt(-a^2 + b^2)*f*sqrt(g*sin(e + f*x))) - (2*sqrt(2)*sqrt(d)*Elliptic.Pi(-(a/(b + sqrt(-a^2 + b^2))), asin(sqrt(d*cos(e + f*x))/(sqrt(d)*sqrt(1 + sin(e + f*x)))), -1)*sqrt(sin(e + f*x)))/(sqrt(-a^2 + b^2)*f*sqrt(g*sin(e + f*x)))]

    @test_int [sqrt(g*sin(e + f*x))/(sqrt(d*cos(e + f*x))*(a + b*cos(e + f*x))), x, 5, -((2*sqrt(2)*sqrt(g)*sqrt(cos(e + f*x))*Elliptic.Pi(-(sqrt(-a + b)/sqrt(a + b)), asin(sqrt(g*sin(e + f*x))/(sqrt(g)*sqrt(1 + cos(e + f*x)))), -1))/(sqrt(-a + b)*sqrt(a + b)*f*sqrt(d*cos(e + f*x)))) + (2*sqrt(2)*sqrt(g)*sqrt(cos(e + f*x))*Elliptic.Pi(sqrt(-a + b)/sqrt(a + b), asin(sqrt(g*sin(e + f*x))/(sqrt(g)*sqrt(1 + cos(e + f*x)))), -1))/(sqrt(-a + b)*sqrt(a + b)*f*sqrt(d*cos(e + f*x)))]
    @test_int [1/(sqrt(g*sin(e + f*x))*sqrt(d*cos(e + f*x))*(a + b*cos(e + f*x))), x, 7, -((2*sqrt(2)*b*Elliptic.Pi(-(a/(b - sqrt(-a^2 + b^2))), asin(sqrt(d*cos(e + f*x))/(sqrt(d)*sqrt(1 + sin(e + f*x)))), -1)*sqrt(sin(e + f*x)))/(a*sqrt(-a^2 + b^2)*sqrt(d)*f*sqrt(g*sin(e + f*x)))) + (2*sqrt(2)*b*Elliptic.Pi(-(a/(b + sqrt(-a^2 + b^2))), asin(sqrt(d*cos(e + f*x))/(sqrt(d)*sqrt(1 + sin(e + f*x)))), -1)*sqrt(sin(e + f*x)))/(a*sqrt(-a^2 + b^2)*sqrt(d)*f*sqrt(g*sin(e + f*x))) + (Elliptic.F(e - pi/4 + f*x, 2)*sqrt(sin(2*e + 2*f*x)))/(a*f*sqrt(d*cos(e + f*x))*sqrt(g*sin(e + f*x)))]
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.2.2 (g sin)^p (a+b cos)^m (c+d cos)^n.jl


text:
const ID = Int

const uidcounter = Counter(0)

# # __init__() = global uidcounter = Counter(0) #= also works =#
# __init__() = reset!(uidcounter)

"Unique id"
# uid() = (global uidcounter; @show increment!(uidcounter))
uid() = increment!(uidcounter)

@spec :nocheck (x = [uid() for i = 1:Inf]; unique(x) == x)

"Construct globally unique id for indices for ω"
macro id()
  uid()
end


text:
function update_vehicle(ev::EVinstance,dt::Float64)

    ev.currentWorkload-=ev.currentPower*dt;
    ev.currentDeadline-=dt;
    ev.currentReportedDeadline-=dt;

end

function compute_average(f,T::Vector{Float64},X::Vector{UInt16})
    return sum(f(X[1:end-1]).*diff(T))/T[end]
end

function compute_statistics!(sim::EVSim,t_start=0.0,t_end=Inf)

    idx = findall(t_start.<=sim.timetrace.T.<t_end)
    sim.stats.avgX = compute_average(x->x,sim.timetrace.T[idx],sim.timetrace.X[idx]);
    sim.stats.avgY = compute_average(x->x,sim.timetrace.T[idx],sim.timetrace.Y[idx]);

    rangeX = collect(minimum(sim.timetrace.X[idx]):maximum(sim.timetrace.X[idx]))
    pX = [compute_average(x->x.==l,sim.timetrace.T[idx],sim.timetrace.X[idx]) for l in rangeX]
    sim.stats.rangeX = rangeX;
    sim.stats.pX=pX;

    rangeY = collect(minimum(sim.timetrace.Y[idx]):maximum(sim.timetrace.Y[idx]))
    pY = [compute_average(x->x.==l,sim.timetrace.T[idx],sim.timetrace.Y[idx]) for l in rangeY]
    sim.stats.rangeY = rangeY;
    sim.stats.pY=pY;

    EVs = filter(ev-> ev.arrivalTime>= t_start && ev.departureTime<=t_end, sim.EVs)
    sim.stats.avgW = mean([ev.departureWorkload for ev in EVs]);

    sim.stats.pD = sum([ev.departureWorkload>0 for ev in EVs])/length(EVs);

    return nothing
end

function get_vehicle_trajectories(sim::EVSim,t_start=0.0,t_end=Inf)

    d=OrderedDict{Float64,Array{Array{Float64,1},1}}();

    snaps = filter(snap -> t_start<=snap.t <= t_end ,sim.snapshots)

    for snapshot in snaps

        charging = snapshot.charging;

        for j=1:length(charging)
            index = charging[j].arrivalTime;
            if haskey(d,index)
                push!(d[index],[snapshot.t-charging[j].arrivalTime;charging[j].currentWorkload;charging[j].currentDeadline])
            else
                d[index] = [];
            end
        end
    end

    return d
end

function fairness_index(sa,s)

    J=1.0;

    if length(sa)>0 && sum(sa)>0
        J = mean(sa./s).^2/(mean((sa./s).^2));
    end
    return J
end

function compute_fairness(sim::EVSim,t::Vector{Float64},h::Float64)

    J=zeros(length(t))

    for i=1:length(t)

        evs = filter(ev->(ev.completionTime<=t[i])&&(ev.completionTime>=t[i]-h),sim.EVs);
        sa = [ev.requestedEnergy-ev.departureWorkload for ev in evs];
        s = [ev.requestedEnergy for ev in evs];

        J[i] = fairness_index(sa,s);

    end
    return J;
end

macro addpolicy(name::String)
    f1 = Symbol("ev_",name);
    f2 = Symbol("ev_",name,"_trace");
    f3 = Symbol("ev_",name,"_uncertain");
    policy = Symbol(name,"_policy");
    eval( quote

        function $f1(lambda,mu,gamma,Tfinal,C=Inf;snapshots=[Inf])
            ev_sim(lambda,mu,gamma,Tfinal,C,$policy,snapshots)
        end

        function $f2(arribos,demandas,salidas,potencias,C=Inf;snapshots=[Inf],salidaReportada=nothing)
            ev_sim_trace(arribos,demandas,salidas,potencias,$policy,C,snapshots,salidaReportada=salidaReportada)
        end

        function $f2(df::DataFrame,C=Inf;snapshots=[Inf])
            ev_sim_trace(df,$policy,C,snapshots)
        end

        function $f3(lambda,mu,gamma,Tfinal,C=Inf,uncertainity_paramter=0.0;snapshots=[Inf])
            ev_sim_uncertain(lambda,mu,gamma,Tfinal,C,$policy,uncertainity_paramter,snapshots)
        end

        export $f1, $f2, $f3

    end)
end

function get_policy_name(policy::Function)

    name = String(Symbol(policy));
    name = split(name,"_")[1];
    uppercase(name);

end

##Genera una traza de arribos Poisson como la que se usa en ev_sim para pasarle
#al ev_sim_trace. De este modo se puede fijar la traza de vehiculos. Devuelve un
#dataframe de arribos.
function generate_Poisson_stream(lambda,mu,gamma,Tfinal)

    t=0.0
    arr = Exponential(1/lambda)
    work = Exponential(1/mu)
    lax = Exponential(1/gamma)

    #initial approximate memory allocation
    n=round(Integer,1.2*lambda*Tfinal)
    arribos = zeros(n)
    demandas = zeros(n)
    salidas = zeros(n)
    #todas las potencias van a ser 1
    potencias = ones(n)

    i=0

    while t<Tfinal

        dt=rand(arr)
        t=t+dt
        i=i+1

        arribos[i] = t
        demandas[i] = rand(work)
        salidas[i] = t+demandas[i] + rand(lax)

    end

    df=DataFrame(:arribos=>arribos[1:i], :demandas=>demandas[1:i], :salidas=>salidas[1:i], :potencias=>potencias[1:i])
    return df
end

function Base.show(ev::EVinstance)

    println("An EV instance with:")
    println("Arrival time: $(ev.arrivalTime)")
    println("Departure time: $(ev.departureTime)")
    println("Self-reported departure time: $(ev.reportedDepartureTime)")
    println("Requested energy: $(ev.requestedEnergy)")
    println("Current remaining work: $(ev.currentWorkload)")
    println("Current remaining deadline: $(ev.currentDeadline)")
    println("Current remaining deadline as reported: $(ev.currentReportedDeadline)")
    println("Current charging rate: $(ev.chargingPower)")
    println("Remaining energy on departure (if departed): $(ev.departureWorkload)")
    println("Comppletion time (if completed): $(ev.completionTime)")

end

function Base.show(sim::EVSim)

    println("An EV simulation with:")
    for key in keys(sim.parameters)
        println("\t$key: \t\t $(sim.parameters[key])")
    end

end

function savesim(sim::EVSim, file::String)
    io=open(file,"w");
    serialize(io,sim);
    close(io);
end

function loadsim(file::String)
    io=open(file,"r");
    sim::EVSim = deserialize(io);
    close(io);
    return sim;
end


text:
"""
    Mountain Car

    https://en.wikipedia.org/wiki/Mountain_car_problem
"""

struct MountainCar{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    gravity
    mass
    steepness
    friction

    xl
    xu
    ul
    uu
end

function fd(model::MountainCar{Discrete, FixedTime}, x⁺, x, u, w, h, t)
    x⁺ - [x[1] + x[2] * h;
          x[2] + (model.gravity * model.mass * cos(model.steepness * x[1])
            + u[1] / model.mass - model.friction * x[2]) * h]
end

function fd(model::MountainCar{Discrete, FixedTime}, x, u, w, h, t)
    [x[1] + x[2] * h;
     x[2] + (model.gravity * model.mass * cos(model.steepness * x[1])
        + u[1] / model.mass - model.friction * x[2]) * h]
end

n, m, d = 2, 1, 0
gravity = 9.81
mass = 0.2
steepness = 3.0
friction = 0.3
xl = [-1.2; -1.5]
xu = [0.5; 1.5]
ul = [-2.0]
uu = [2.0]
model = MountainCar{Discrete, FixedTime}(n, m, d,
    gravity, mass, steepness, friction, xl, xu, ul, uu)

function visualize!(vis, model::Car, x;
       Δt = 0.1, traj = false)

	default_background!(vis)

	obj_path = joinpath(pwd(),"models/cybertruck/cybertruck.obj")
	mtl_path = joinpath(pwd(),"models/cybertruck/cybertruck.mtl")

	m = ModifiedMeshFileObject(obj_path,mtl_path,scale=0.05)
	setobject!(vis["car"],m)

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

    for t = 1:length(x)
        MeshCat.atframe(anim,t) do
			settransform!(vis["car"],
				compose(Translation(x[t][1], x[t][2], 0.0),
					LinearMap(RotZ(x[t][3] + pi) * RotX(pi / 2.0))))
        end
    end

	if traj
		pts_nom = collect(eachcol(hcat(x...)))
		material_nom = LineBasicMaterial(color = colorant"orange", linewidth = 3.0)
		setobject!(vis["traj"], Object(PointCloud(pts_nom), material_nom, "Line"))
	end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, 0.0, 0.0),LinearMap(RotY(-pi/2.5))))

    MeshCat.setanimation!(vis, anim)
end


text:
### factor_Zp code

## 14.3
## f in Zq[x]
## return g1, g2, ..., gs: g) is product of all monic irreducible polys in Fq[x] of degree i that divide f
function distinct_degree_factorization(f, q, x=gen(parent(f)))
    degree(f) > 0 || error("f is a constant poly")

    h, f0, i = x, f, 0
    gs = eltype(f)[]
    while true
        i += 1
        h = _powermod(h, q, f0)
        g = gcd(h - x, f)
        f = divexact(f, g)
        push!(gs, g)
        if degree(f) < 2(i+1)
            push!(gs, f)
            break
        end
    end
    gs
end

#function _equal_degree_splitting(f::AbstractAlgebra.Generic.Poly{AbstractAlgebra.gfelem{T}}, q, x, d) where {T}
function _equal_degree_splitting(f, q, x, d) #where {T}
    T = eltype(f)
    n = degree(f)
    n <= 1 && return (f, false)

    # random poly of degree < n
    q′ = Int(q)
    a = sum(convert(T, rand(0:(q′-1))) * x^i for i in 0:(n-1))
    degree(a) <= 0 && return (f, false)

    g1 = gcd(a, f)
    !isone(g1) && return (g1, true)

    b = _powermod(a, (q^d-1) ÷ 2, f)
    g2 = gcd(b-1, f)

    !isone(g2) && g2 != f && return (g2, true)

    return (f, false)
end

# Algo 14.8
# f square free, monic in Fq[x]
# q odd prime power
# d divides n, all irreducible factors of f have degree d
# this calls _equal_degree_splitting which has probability 2^(1-r) of success
# returns a proper monic factor (g,true) or (f,false)
function equal_degree_splitting(f, q, x, d)
    K = 50
    while K > 0
        g, val = _equal_degree_splitting(f, q, x, d)
        val && return (g, true)
        K -= 1
    end
    return (f, false)
end

# f square free in F_q[x]
# d divides degree(f)
# all irreducible factors of f have degree d
# return all monic factors of f in F_q[x]
function equal_degree_factorization(f, q, x, d)
    degree(f) == 0 && return (f, )
    degree(f) == d && return (f, )
    g, val = equal_degree_splitting(f, q, x, d)

    if val
        return tuplejoin(equal_degree_factorization(g, q, x, d),
                         equal_degree_factorization(divexact(f, g), q, x, d))
    else
        return (f,)
    end
end

# Cantor Zazzenhaus factoring of a polynomial of Zp[x]
# Algorithm 14.13 from von Zur Gathen and Gerhard, Modern Computer Algebra v1, 1999
# f in Zq[x], q prime
function factor_Zp(f, q, x=variable(f))

    h = x
    v = f * inv(lead(f)) # monic
    i = 0
    U = Dict{typeof(f), Int}()

    while true
        i += 1
        h = _powermod(h, q, f)
        g = gcd(h - x, v)

        if !isone(g)
            gs = equal_degree_factorization(g, q, x, i)
            for (j,gj) in enumerate(gs)
                e = 0
                while true
                    qu,re = divrem(v, gj)
                    if iszero(re)
                        v = qu
                        e += 1
                    else
                        break
                    end
                end
                U[gj] = e
            end
        end
        isone(v) && break
    end
    U
end

# slightly faster?
# this saves a `divrem` per factor
function factor_Zp_squarefree(f, q, x=variable(f))
    h = x
    v = f * inv(lead(f)) # monic
    i = 0

    facs = typeof(f)[]
    while true
        i += 1
        h = _powermod(h, q, f)
        g = gcd(h - x, v)

        if !isone(g)
            gs = equal_degree_factorization(g, q, x, i)
            for (j,gj) in enumerate(gs)
                qu,re = divrem(v, gj)
                if iszero(re)
                    push!(facs, gj)
                    v = qu
                end
            end
        end
        isone(v) && break
    end
    facs
end

function factormod(f::AbstractAlgebra.Generic.Poly{T}, p::Integer) where {T <: Integer}
    x = string(var(parent(f)))
    fp = as_poly_Zp(f, p, x)
    c = lead(fp)
    fp = fp * inv(c)
    us = factor_Zp(fp, p)
    if !isone(c)
        us[c*one(fp)] = 1
    end
    us
end

function factormod(fs::Vector{T}, p::Integer) where {T <: Integer}
    f = as_poly_Zp(fs, p, "x")
    factor_Zp(fp, p)
end


text:
using DualNumbers
import DualNumbers: Dual4, real, epsilon1, epsilon2, epsilon3, epsilon4
import NaNMath
using Base.Test

x = DualNumbers.Dual4(2, 1, 2, 3, 4)
y = x^3

@test_approx_eq real(y) 2.0^3
@test_approx_eq epsilon1(y) 3.0*2^2
@test_approx_eq epsilon2(y) 2*3.0*2^2
@test_approx_eq epsilon3(y) 3*3.0*2^2
@test_approx_eq epsilon4(y) 4*3.0*2^2

y = x^3.0

@test_approx_eq real(y) 2.0^3
@test_approx_eq epsilon1(y) 3.0*2^2
@test_approx_eq epsilon2(y) 2*3.0*2^2
@test_approx_eq epsilon3(y) 3*3.0*2^2
@test_approx_eq epsilon4(y) 4*3.0*2^2

y = NaNMath.pow(x,3)

@test_approx_eq real(y) 2.0^3
@test_approx_eq epsilon1(y) 3.0*2^2
@test_approx_eq epsilon2(y) 2*3.0*2^2
@test_approx_eq epsilon3(y) 3*3.0*2^2
@test_approx_eq epsilon4(y) 4*3.0*2^2

y = sin(x)+exp(x)
@test_approx_eq real(y) sin(2)+exp(2)
@test_approx_eq epsilon1(y) cos(2)+exp(2)
@test_approx_eq epsilon2(y) 2*(cos(2)+exp(2))
@test_approx_eq epsilon3(y) 3*(cos(2)+exp(2))
@test_approx_eq epsilon4(y) 4*(cos(2)+exp(2))

y = x/2
@test_approx_eq real(y) 1.0
@test_approx_eq epsilon1(y) 0.5
@test_approx_eq epsilon2(y) 1.0
@test_approx_eq epsilon3(y) 3/2
@test_approx_eq epsilon4(y) 2.0

@test x - 2 == Dual4(0,1,2,3,4)

@test x > 1
y = abs(-x)
@test_approx_eq real(y) 2.0
@test_approx_eq epsilon1(y) 1.0
@test_approx_eq epsilon2(y) 2.0
@test_approx_eq epsilon3(y) 3.0
@test_approx_eq epsilon4(y) 4.0

@test isequal(1.0,Dual4(1.0))

y = 1/x
@test_approx_eq real(y) 1/2
@test_approx_eq epsilon1(y) -1/2^2
@test_approx_eq epsilon2(y) -2/2^2
@test_approx_eq epsilon3(y) -3/2^2
@test_approx_eq epsilon4(y) -4/2^2

function squareroot(x)
    it = x
    while abs(it*it - x) > 1e-13
        it = it - (it*it-x)/(2it)
    end
    return it
end

sq = squareroot(Dual4(10000.0,1.0,2.0,3.0,4.0))
@test_approx_eq epsilon1(sq) 0.005
@test_approx_eq epsilon2(sq) 2*0.005
@test_approx_eq epsilon3(sq) 3*0.005
@test_approx_eq epsilon4(sq) 4*0.005

@test_approx_eq epsilon1(exp(1)^Dual4(1.0,1.0,2.0,3.0,4.0)) exp(1)
@test_approx_eq epsilon2(exp(1)^Dual4(1.0,1.0,2.0,3.0,4.0)) 2exp(1)
@test_approx_eq epsilon3(exp(1)^Dual4(1.0,1.0,2.0,3.0,4.0)) 3exp(1)
@test_approx_eq epsilon4(exp(1)^Dual4(1.0,1.0,2.0,3.0,4.0)) 4exp(1)
@test_approx_eq epsilon2(NaNMath.sin(Dual4(1.0,1.0,2.0,3.0,4.0))) 2cos(1)

@test_approx_eq epsilon1(NaNMath.pow(exp(1),Dual4(1.0,1.0,2.0,3.0,4.0))) exp(1)
@test_approx_eq epsilon2(NaNMath.pow(exp(1),Dual4(1.0,1.0,2.0,3.0,4.0))) 2exp(1)
@test_approx_eq epsilon3(NaNMath.pow(exp(1),Dual4(1.0,1.0,2.0,3.0,4.0))) 3exp(1)
@test_approx_eq epsilon4(NaNMath.pow(exp(1),Dual4(1.0,1.0,2.0,3.0,4.0))) 4exp(1)


text:
@testset "Helmholtz" begin

# Only on one processor...
# Start by importing both `PetscWrap`, for the distributed matrices, and `SlepcWrap` for the eigenvalues.
using PetscWrap
using SlepcWrap

# Number of mesh points and mesh step
n = 21
Δx = 1. / (n - 1)

# Initialize SLEPc. Either without arguments, calling `SlepcInitialize()` or using "command-line" arguments.
# To do so, either provide the arguments as one string, for instance
# `SlepcInitialize("-eps_max_it 100 -eps_tol 1e-5")` or provide each argument in
# separate strings : `PetscInitialize(["-eps_max_it", "100", "-eps_tol", "1e-5")`.
# Here we ask for the five closest eigenvalues to ``0``, using a non-zero pivot for the LU factorization and a
# "shift-inverse" process.
SlepcInitialize("-eps_target 0 -eps_nev 5 -st_pc_factor_shift_type NONZERO -st_type sinvert")

# Create the problem matrices, set sizes and apply "command-line" options. Note that we should
# set the number of preallocated non-zeros to increase performance.
A = MatCreate()
B = MatCreate()
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)
MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, n, n)
MatSetFromOptions(A)
MatSetFromOptions(B)
MatSetUp(A)
MatSetUp(B)

# Get rows handled by the local processor
A_rstart, A_rend = MatGetOwnershipRange(A)
B_rstart, B_rend = MatGetOwnershipRange(B)

# Fill matrix A  with second order derivative central scheme
for i in A_rstart:A_rend-1
    if (i == 0)
        MatSetValues(A, [0], [0, 1], [-2., 1] / Δx^2, INSERT_VALUES)
    elseif (i == n-1)
        MatSetValues(A, [n-1], [n-2, n-1], [1., -2.] / Δx^2, INSERT_VALUES)
    else
        MatSetValues(A, [i], i-1:i+1, [1., -2., 1.] / Δx^2, INSERT_VALUES)
    end
end

# Fill matrix B with identity matrix
for i in B_rstart:B_rend-1
    MatSetValue(B, i, i, -1., INSERT_VALUES)
end

# Set boundary conditions : u(0) = 0 and u(1) = 0. Only the processor
# handling the corresponding rows are playing a role here.
(A_rstart == 0) && MatSetValues(A, [0], [0,1], [1., 0.], INSERT_VALUES)
(B_rstart == 0) && MatSetValue(B, 0, 0, 0., INSERT_VALUES)

(A_rend == n) && MatSetValues(A, [n-1], [n-2,n-1], [0., 1.], INSERT_VALUES)
(B_rend == n) && MatSetValue(B, n-1, n-1, 0., INSERT_VALUES)

# Assemble the matrices
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)
MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)
MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)

# Now we set up the eigenvalue solver
eps = EPSCreate()
EPSSetOperators(eps, A, B)
EPSSetFromOptions(eps)
EPSSetUp(eps)

# Then we solve
EPSSolve(eps)

# And finally we can inspect the solution. Let's first get the number of converged eigenvalues:
nconv = EPSGetConverged(eps)

# Then we can get/display these eigenvalues (more precisely their square root, i.e ``\simeq \omega``)
expected = (
    3.138363829113791,
    6.2573786016092345,
    9.337814554236218,
    12.360679774997893
)
for ieig in 0:3
    vpr, vpi = EPSGetEigenvalue(eps, ieig)
    @test isapprox(√(real(vpr)), expected[ieig + 1])

end

# We can also play with eigen vectors. First, create two Petsc vectors to allocate memory
vecr, veci = MatCreateVecs(A)

# Then loop over the eigen pairs and retrieve eigenvectors
for ieig in 0:nconv-1
    vpr, vpi, vecpr, vecpi = EPSGetEigenpair(eps, ieig, vecr, veci)

    ## At this point, you can call VecGetArray to obtain a Julia array (see PetscWrap examples).
    ## If you are on one processor, you can even plot the solution to check that you have a sinus
    ## solution. On multiple processors, this would require to "gather" the solution on one processor only.
end

# Finally, let's free the memory
MatDestroy(A)
MatDestroy(B)
EPSDestroy(eps)

# And call finalize when you're done
SlepcFinalize()

# Test if we reached this point
@test true

end


text:
@testset "prepareImage() input varieties" begin
    yolomod = YOLO.v3_tiny_COCO(w=416, h=416, silent=true)
    batch = emptybatch(yolomod)
    @testset "2D float32" begin
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 416, 416), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 500, 500), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 200, 500), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 500, 200), yolomod)
        @test true
    end
    @testset "3D float32 with 1 channel" begin
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 416, 416, 1), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 500, 500, 1), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 200, 500, 1), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 500, 200, 1), yolomod)
        @test true
    end
    @testset "3D float32 with 3 channels" begin
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 416, 416, 3), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 500, 500, 3), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 200, 500, 3), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Float32, 500, 200, 3), yolomod)
        @test true
    end
    @testset "2D Gray" begin
        batch[:,:,:,1], padding = prepareImage(rand(Gray, 416, 416), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Gray, 500, 500), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Gray, 200, 500), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(Gray, 500, 200), yolomod)
        @test true
    end
    @testset "2D RGB" begin
        batch[:,:,:,1], padding = prepareImage(rand(RGB, 416, 416), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(RGB, 500, 500), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(RGB, 200, 500), yolomod)
        batch[:,:,:,1], padding = prepareImage(rand(RGB, 500, 200), yolomod)
        @test true
    end
end


text:
"""
    pcd!(rbm, data)

Trains the RBM on data using Persistent Contrastive divergence.
"""
function pcd!(
    rbm::RBM,
    data::AbstractArray;
    batchsize::Int = 1,
    epochs::Int = 1,
    wts::Union{AbstractVector, Nothing} = nothing, # data weights
    steps::Int = 1, # MC steps to update fantasy chains
    optim = default_optimizer(_nobs(data), batchsize, epochs), # optimization algorithm
    stats = suffstats(rbm, data; wts), # sufficient statistics for visible layer

    # regularization
    l2_fields::Real = 0, # visible fields L2 regularization
    l1_weights::Real = 0, # weights L1 regularization
    l2_weights::Real = 0, # weights L2 regularization
    l2l1_weights::Real = 0, # weights L2/L1 regularization

    # gauge
    zerosum::Bool = true, # zerosum gauge for Potts layers
    rescale::Bool = true, # normalize hidden units to var(h) = 1
    center::Bool = true, # center gradients

    # momentum for hidden activity statistics tracking
    ρh::Real = 99//100,
    ϵh::Real = 1//100, # prevent vanishing var(h) estimate

    callback = empty_callback, # called for every batch
    mode::Symbol = :pcd, # :pcd, :cd, or :exact

    vm = fantasy_init(rbm; batchsize, mode), # fantasy chains
    shuffle::Bool = true
)
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    @assert isnothing(wts) || size(data)[end] == length(wts)
    @assert ϵh ≥ 0

    # we center layers with their average activities
    ave_v = batchmean(visible(rbm), data; wts)
    ave_h, var_h = total_meanvar_from_inputs(hidden(rbm), inputs_h_from_v(rbm, data); wts)
    @assert all(var_h .+ ϵh .> 0)

    # gauge constraints
    zerosum && zerosum!(rbm)
    rescale && rescale_hidden!(rbm, sqrt.(var_h .+ ϵh))

    # store average weight of each data point
    wts_mean = isnothing(wts) ? 1 : mean(wts)

    for epoch in 1:epochs
        batches = minibatches(data, wts; batchsize, shuffle)
        for (batch_idx, (vd, wd)) in enumerate(batches)
            # compute positive gradient from data
            ∂d = ∂free_energy(rbm, vd; wts = wd, stats)

            # compute negative gradient
            ∂m = ∂logpartition(rbm; vd, vm, wd, mode, steps)

            # likelihood gradient is the difference of positive and negative parts
            ∂ = subtract_gradients(∂d, ∂m)

            # correct weighted minibatch bias
            batch_weight = isnothing(wts) ? 1 : mean(wd) / wts_mean
            ∂ = gradmult(∂, batch_weight)

            # extract hidden unit statistics from gradient
            ave_h_batch = grad2ave(rbm.hidden, ∂d.hidden)
            var_h_batch = grad2var(rbm.hidden, ∂d.hidden)

            #= Exponential moving average of mean and variance of hidden unit activations.
            The batchweight can be interpreted as an "effective number of updates". =#
            ρh_eff = ρh ^ batch_weight # effective damp after 'batch_weight' updates
            ave_h .= ρh_eff * ave_h .+ (1 - ρh_eff) * ave_h_batch
            var_h .= ρh_eff * var_h .+ (1 - ρh_eff) * var_h_batch
            @assert all(var_h .+ ϵh .> 0)

            # weight decay
            ∂regularize!(∂, rbm; l2_fields, l1_weights, l2_weights, l2l1_weights, zerosum)

            if center
                # gradient of the centered RBM (Melchior et al 2016)
                ∂ = center_gradient(rbm, ∂, ave_v, ave_h)
            end

            # compute parameter update step from gradient, according to optimizer algorithm
            update!(∂, rbm, optim)

            if center
                # transform step in centered coordinates to uncentered coordinates
                ∂ = uncenter_step(rbm, ∂, ave_v, ave_h)
            end

            # update parameters with update step computed above
            update!(rbm, ∂)

            # reset gauge
            zerosum && zerosum!(rbm)
            rescale && rescale_hidden!(rbm, sqrt.(var_h .+ ϵh))

            callback(; rbm, optim, epoch, batch_idx, vm, vd, wd)
        end
    end
    return rbm
end

# init fantasy chains
fantasy_init(rbm::RBM; batchsize::Int, mode::Symbol = :pcd) = fantasy_init(rbm.visible; batchsize, mode)

function fantasy_init(l::AbstractLayer; batchsize::Int, mode::Symbol = :pcd)
    @assert mode ∈ (:pcd, :cd, :exact)
    if mode === :exact
        @warn "Running extensive sampling; this can take a lot of RAM and time"
        return extensive_sample(l)
    else
        return sample_from_inputs(l, falses(size(l)..., batchsize))
    end
end

empty_callback(@nospecialize(args...); @nospecialize(kw...)) = nothing

function ∂logpartition(
    rbm::RBM;
    vd::AbstractArray,
    vm::AbstractArray,
    wd::Union{AbstractArray, Nothing},
    mode::Symbol,
    steps::Int
)
    @assert mode ∈ (:pcd, :cd, :exact)
    if mode === :pcd
        # update persistent fantasy chains
        vm .= sample_v_from_v(rbm, vm; steps)
        return ∂free_energy(rbm, vm)
    elseif mode === :cd
        # chains re-initialized from data
        vm .= sample_v_from_v(rbm, vd; steps)
        return ∂free_energy(rbm, vm; wts = wd)
    elseif mode === :exact
        # use exact RBM probabilities (assuming `vm` is an extensive sample)
        p = softmax(-free_energy(rbm, vm))
        return ∂free_energy(rbm, vm; wts = p)
    end
end

function extensive_sample(layer::Binary; maxlen::Int = 12)
    @assert length(layer) ≤ maxlen
    N = length(layer)
    v = reduce(hcat, BitVector.(digits.(Bool, 0:(2^N - 1), base=2, pad=N)))
    return reshape(v, size(layer)..., :)
end

function extensive_sample(layer::Spin; maxlen::Int = 12)
    σ = extensive_sample(Binary(layer.θ); maxlen)
    return Int8(2) * σ .- Int8(1)
end

function extensive_sample(layer::Potts; maxlen::Int = 12)
    q = size(layer, 1)
    N = prod(tail(size(layer)))
    @assert N * log2(q) ≤ maxlen && q < typemax(Int8)
    potts = reduce(hcat, digits.(Int8, 0:(q^N - 1), base=q, pad=N))
    onehot = reshape(potts, 1, :) .== 0:(q - 1)
    return reshape(onehot, size(layer)..., :)
end


text:
@testset "MultiArmBanditsEnv" begin

    rng = StableRNG(123)
    env = MultiArmBanditsEnv(; rng = rng)
    rewards = []

    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

    N = 50_000
    for _ in 1:N
        while !is_terminated(env)
            env(rand(rng, legal_action_space(env)))
        end
        push!(rewards, reward(env))
        reset!(env)
    end

    @test isapprox(mean(rewards), mean(env.true_values); atol = 0.01)

end


text:
################################################################################
# Young diagrams, partitions of unity and characters of the symmetric group Sn #
################################################################################

typealias Partition Vector{Int64}
typealias YoungDiagram Array{Int64,2}
typealias SkewDiagram (Partition, Partition)

export Partition,
       YoungDiagram, #represents shape of Young diagram
       SkewDiagram,  #skew diagrams
       partitionsequence,
       isrimhook,    #Check if skew diagram is rim hook
       leglength,
       character     #Computes character of irrep of Sn

import Base.\

#################
# Skew diagrams #
#################

#This uses a very simple internal representation for skew diagrams
\(λ::Partition, μ::Partition) = MakeSkewDiagram(λ, μ)
function MakeSkewDiagram(λ::Partition, μ::Partition)
    m, n = length(λ), length(μ)
    if n>m error("Cannot construct skew diagram") end
    (λ, μ)
end

#Checks if skew diagram is a rim hook
isrimhook(λ::Partition, μ::Partition)=isrimhook(λ \ μ)
function isrimhook(ξ::SkewDiagram)
    λ, μ = ξ
    m, n = length(λ), length(μ)
    if n>m error("Cannot construct skew diagram") end
    #Construct matrix representation of diagram
    #XXX This is a horribly inefficient way of checking condition 1!
    l = maximum(λ)
    youngdiagram=zeros(Int64, m, l)
    for i=1:n
        youngdiagram[i, μ[i]+1:λ[i]]=1
    end
    for i=n+1:m
        youngdiagram[i, 1:λ[i]]=1
    end
    #Condition 1. Must be edgewise connected
    youngdiagramList=[]
    for i=1:m
        for j=1:l
            if youngdiagram[i, j]==1 youngdiagramList = [youngdiagramList; (i, j)] end
        end
    end
    for k=1:length(youngdiagramList)
        i, j = youngdiagramList[k]
        numNeighbors = 0
        for kp=1:length(youngdiagramList)
            ip,jp= youngdiagramList[kp]
            if abs(i-ip) + abs(j-jp) == 1 numNeighbors += 1 end
        end
        if numNeighbors == 0 return false end #Found a cell with no adjacent neighbors
    end
    #Condition 2. Must not contain 2x2 square of cells
    for i=1:m-1
        for j=1:l-1
            if youngdiagram[i, j]== 0 continue end
            if youngdiagram[i, j+1] == youngdiagram[i+1, j] == youngdiagram[i+1, j+1] == 1
                return false
            end
        end
    end
    return true
end

#Strictly speaking, defined for rim hook only, but here we define it for all skew diagrams
leglength(λ::Partition, μ::Partition)=leglength((λ \ μ))
function leglength(ξ::SkewDiagram)
    λ, μ = ξ
    m, n = length(λ), length(μ)
    #Construct matrix representation of diagram
    l = maximum(λ)
    youngdiagram=zeros(Int64, m, l)
    for i=1:n
        youngdiagram[i, μ[i]+1:λ[i]]=1
    end
    for i=n+1:m
        youngdiagram[i, 1:λ[i]]=1
    end
    for i=m:-1:1
        if any([x==1 for x in youngdiagram[i,:]]) return i-1 end
    end
    return -1 #If entire matrix is empty
end

#######################
# partition sequences #
#######################

#Computes essential part of the partition sequence of lambda
function partitionsequence(lambda::Partition)
    Λ▔ = Int64[]
    λ = [lambda; 0]
    m = length(lambda)
    for i=m:-1:1
        for k=1:(λ[i]-λ[i+1])
            Λ▔ = [Λ▔; 1]
        end
        Λ▔ = [Λ▔ ; 0]
    end
    Λ▔
end

#This takes two elements of a partition sequence, with a to the left of b
isrimhook(a::Int64, b::Int64) = (a==1) && (b==0)

#############################
# Character of irreps of Sn #
#############################

#Computes recursively using the Murnaghan-Nakayama rule.
function MN1inner(R::Vector{Int64}, T::Dict, μ::Partition, t::Integer)
    s=length(R)
    χ::Integer=1
    if t<=length(μ)
        χ, σ::Integer = 0, 1
        for j=1:μ[t]-1
            if R[j]==0 σ=-σ end
        end
        for i=1:s-μ[t]
            if R[i] != R[i+μ[t]-1] σ=-σ end
            if isrimhook(R[i], R[i+μ[t]])
                R[i], R[i+μ[t]] = R[i+μ[t]], R[i]
                rhohat = R[i:i+μ[t]]
                if !haskey(T, rhohat) #Cache result in lookup table
                    T[rhohat] = MN1inner(R, T, μ, t+1)
                end
                χ += σ * T[rhohat]
                R[i], R[i+μ[t]] = R[i+μ[t]], R[i]
            end
        end
    end
    χ
end

#Computes character $χ^λ(μ)$ of the partition μ in the λth irrep of the
#symmetric group $S_n$
#
#Implements the Murnaghan-Nakayama algorithm as described in:
#    Dan Bernstein,
#    "The computational complexity of rules for the character table of Sn",
#    Journal of Symbolic Computation, vol. 37 iss. 6 (2004), pp 727-748.
#    doi:10.1016/j.jsc.2003.11.001
function character(λ::Partition, μ::Partition)
    T = {()=>0} #Sparse array implemented as dict
    Λ▔ = partitionsequence(λ)
    MN1inner(Λ▔, T, μ, 1)
end


text:
using ConvexTests, COSMO

@info "Starting COSMO tests"

do_tests("COSMO", () -> COSMO.Optimizer(verbose=false, eps_abs=1e-6, eps_rel=1e-6); exclude = [r"mip", r"dual"], description = "Tests run with `eps_abs=1e-6` and `eps_rel=1e-6`.")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors, Materials
parameters = ChabocheParameterState(E = 200.0e3,
                                    nu = 0.3,
                                    R0 = 100.0,
                                    Kn = 100.0,
                                    nn = 10.0,
                                    C1 = 10000.0,
                                    D1 = 100.0,
                                    C2 = 50000.0,
                                    D2 = 1000.0,
                                    Q = 0.0,
                                    b = 0.1)
material = Chaboche(parameters = parameters)

dtime = 0.25
times = [material.drivers.time]
stresses = [copy(tovoigt(material.variables.stress))]
strains = [copy(tovoigt(material.drivers.strain; offdiagscale=2.0))]
cumeqs = [copy(material.variables.cumeq)]

stresses_expected = [[       50.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [      100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [      150.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [      150.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [      100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [     -100.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
dstress = 50.0

dtimes      = [dtime, dtime, dtime, 1e3, dtime, 1e3]
dstresses11  = [dstress, dstress, dstress, 0.0, -dstress, -4*dstress]

for i in 1:length(dtimes)
    dstress11 = dstresses11[i]
    dtime = dtimes[i]
    stress_driven_uniaxial_increment!(material, dstress11, dtime)
    update_material!(material)
    push!(times, material.drivers.time)
    push!(stresses, copy(tovoigt(material.variables.stress)))
    push!(strains, copy(tovoigt(material.drivers.strain; offdiagscale=2.0)))
    push!(cumeqs, copy(material.variables.cumeq))
    @test isapprox(tovoigt(material.variables.stress), stresses_expected[i]; atol=1e-4)
end

dstrain_creep = strains[5] - strains[4]
@test isapprox(dstrain_creep[2], -dstrain_creep[1]*0.5; atol=1e-4)
@test isapprox(dstrain_creep[3], -dstrain_creep[1]*0.5; atol=1e-4)

dcumeq = cumeqs[end] - cumeqs[end-1]
@test dcumeq > 0


text:
# This file was generated, do not modify it. # hide
sort!(cap_sum_plot, :capacity_mw_sum, rev=true)

figure(figsize=(8,6))

plt.bar(cap_sum_plot.country, cap_sum_plot.capacity_mw_sum, width=0.35)
plt.xticks(rotation=90)

savefig(joinpath(@OUTPUT, "D0-processing-g1.svg")) # hide


text:
using Test
using StaticArrays
using RootedTrees

using RootedTrees.Latexify: latexify

using Plots: Plots, plot
Plots.unicodeplots()

@testset "RootedTrees" begin

@testset "comparisons etc." begin
  trees = (rootedtree([1, 2, 3]),
           rootedtree([1, 2, 3]),
           rootedtree([1, 2, 2]),
           rootedtree([1, 2, 3, 3]),
           rootedtree(Int[]))
  trees_shifted = (rootedtree([1, 2, 3]),
                   rootedtree([2, 3, 4]),
                   rootedtree([1, 2, 2]),
                   rootedtree([1, 2, 3, 3]),
                   rootedtree(Int[]))

  for (t1,t2,t3,t4,t5) in (trees, trees_shifted)
    @test t1 == t1
    @test t1 == t2
    @test !(t1 == t3)
    @test !(t1 == t4)
    @test !(t1 == t5)
    @test !(t2 == t5)
    @test !(t3 == t5)
    @test !(t4 == t5)
    @test t5 == t5

    @test hash(t1) == hash(t1)
    @test hash(t1) == hash(t2)
    @test !(hash(t1) == hash(t3))
    @test !(hash(t1) == hash(t4))
    @test hash(t1) != hash(t5)
    @test hash(t2) != hash(t5)
    @test hash(t3) != hash(t5)
    @test hash(t4) != hash(t5)
    @test hash(t5) == hash(t5)

    @test !(t1 < t1)
    @test !(t1 < t2)
    @test !(t2 < t1)
    @test !(t1 > t2)
    @test !(t2 > t1)
    @test t3 < t2    && t2 > t3
    @test !(t2 < t3) && !(t3 > t2)
    @test t1 < t4    && t4 > t1
    @test !(t4 < t1) && !(t1 > t4)
    @test t1 <= t2   && t2 >= t1
    @test t2 <= t2   && t2 >= t2
    @test t5 < t1
    @test t1 > t5
    @test !(t5 < t5)
    @test !(t1 < t5)

    println(devnull, t1)
    println(devnull, t2)
    println(devnull, t3)
    println(devnull, t4)
  end

  # more tests of the canonical representation
  t = rootedtree([1, 2, 3, 2, 3, 3, 2])
  @test t.level_sequence == [1, 2, 3, 3, 2, 3, 2]
  @test !isempty(t)

  t = rootedtree([1, 2, 3, 2, 3, 4, 2, 3])
  @test t.level_sequence == [1, 2, 3, 4, 2, 3, 2, 3]
  @test !isempty(t)

  t = rootedtree([1, 2, 3, 2, 3, 3, 2, 3])
  @test t.level_sequence == [1, 2, 3, 3, 2, 3, 2, 3]
  @test !isempty(t)

  @test isempty(rootedtree(Int[]))
  @test isempty(empty(t))

  level_sequence = zeros(Int, RootedTrees.BUFFER_LENGTH + 1)
  level_sequence[1] -= 1
  @inferred rootedtree(level_sequence)
end

@testset "hashing" begin
  hashes = [hash(rootedtree(Int[]))]
  for o in 1:12
    for t in RootedTreeIterator(o)
      new_hash = @inferred hash(t)
      @test !(new_hash in hashes)
      push!(hashes, new_hash)
    end
  end
  t = rootedtree([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 2])
  new_hash = @inferred hash(t)
  @test !(new_hash in hashes)
end

# see Table 301(I) etc. in butcher2016numerical
@testset "functions on trees" begin
  t1 = rootedtree([1])
  @test order(t1) == 1
  @test σ(t1) == 1
  @test γ(t1) == 1
  @test α(t1) == 1
  @test β(t1) == α(t1)*γ(t1)
  @test butcher_representation(t1) == "τ"
  latex_string = "\\rootedtree[]"
  @test RootedTrees.latexify(t1) == latex_string
  @test latexify(t1) == latex_string
  @test isempty(RootedTrees.subtrees(t1))
  @test RootedTrees.latexify(empty(t1)) == "\\varnothing"

  @inferred order(t1)
  @inferred σ(t1)
  @inferred γ(t1)
  @inferred α(t1)
  @inferred β(t1)
  @inferred t1 ∘ t1

  t2 = rootedtree([1, 2])
  @test order(t2) == 2
  @test σ(t2) == 1
  @test γ(t2) == 2
  @test α(t2) == 1
  @test β(t2) == α(t2)*γ(t2)
  @test t2 == t1 ∘ t1
  @test butcher_representation(t2) == "[τ]"
  latex_string = "\\rootedtree[[]]"
  @test RootedTrees.latexify(t2) == latex_string
  @test latexify(t2) == latex_string
  @test RootedTrees.subtrees(t2) == [rootedtree([2])]

  t3 = rootedtree([1, 2, 2])
  @test order(t3) == 3
  @test σ(t3) == 2
  @test γ(t3) == 3
  @test α(t3) == 1
  @test β(t3) == α(t3)*γ(t3)
  @test t3 == t2 ∘ t1
  @test butcher_representation(t3) == "[τ²]"
  latex_string = "\\rootedtree[[][]]"
  @test RootedTrees.latexify(t3) == latex_string
  @test latexify(t3) == latex_string
  @test RootedTrees.subtrees(t3) == [rootedtree([2]), rootedtree([2])]

  t4 = rootedtree([1, 2, 3])
  @test order(t4) == 3
  @test σ(t4) == 1
  @test γ(t4) == 6
  @test α(t4) == 1
  @test β(t4) == α(t4)*γ(t4)
  @test t4 == t1 ∘ t2
  @test butcher_representation(t4) == "[[τ]]"
  latex_string = "\\rootedtree[[[]]]"
  @test RootedTrees.latexify(t4) == latex_string
  @test latexify(t4) == latex_string
  @test RootedTrees.subtrees(t4) == [rootedtree([2, 3])]

  t5 = rootedtree([1, 2, 2, 2])
  @test order(t5) == 4
  @test σ(t5) == 6
  @test γ(t5) == 4
  @test α(t5) == 1
  @test β(t5) == α(t5)*γ(t5)
  @test t5 == t3 ∘ t1
  @test butcher_representation(t5) == "[τ³]"
  @test RootedTrees.subtrees(t5) == [rootedtree([2]), rootedtree([2]), rootedtree([2])]

  t6 = rootedtree([1, 2, 2, 3])
  @inferred RootedTrees.subtrees(t6)
  @test order(t6) == 4
  @test σ(t6) == 1
  @test γ(t6) == 8
  @test α(t6) == 3
  @test β(t6) == α(t6)*γ(t6)
  @test t6 == t2 ∘ t2 == t4 ∘ t1
  @test butcher_representation(t6) == "[[τ]τ]"
  @test RootedTrees.subtrees(t6) == [rootedtree([2, 3]), rootedtree([2])]

  t7 = rootedtree([1, 2, 3, 3])
  @test order(t7) == 4
  @test σ(t7) == 2
  @test γ(t7) == 12
  @test β(t7) == α(t7)*γ(t7)
  @test t7 == t1 ∘ t3
  @test butcher_representation(t7) == "[[τ²]]"

  t8 = rootedtree([1, 2, 3, 4])
  @test order(t8) == 4
  @test σ(t8) == 1
  @test γ(t8) == 24
  @test α(t8) == 1
  @test t8 == t1 ∘ t4
  @test butcher_representation(t8) == "[[[τ]]]"

  t9 = rootedtree([1, 2, 2, 2, 2])
  @test order(t9) == 5
  @test σ(t9) == 24
  @test γ(t9) == 5
  @test α(t9) == 1
  @test β(t9) == α(t9)*γ(t9)
  @test t9 == t5 ∘ t1
  @test butcher_representation(t9) == "[τ⁴]"

  t10 = rootedtree([1, 2, 2, 2, 3])
  @test order(t10) == 5
  @test σ(t10) == 2
  @test γ(t10) == 10
  @test α(t10) == 6
  @test β(t10) == α(t10)*γ(t10)
  @test t10 == t3 ∘ t2 == t6 ∘ t1
  @test butcher_representation(t10) == "[[τ]τ²]"

  t11 = rootedtree([1, 2, 2, 3, 3])
  @test order(t11) == 5
  @test σ(t11) == 2
  @test γ(t11) == 15
  @test α(t11) == 4
  @test t11 == t2 ∘ t3 == t7 ∘ t1
  @test butcher_representation(t11) == "[[τ²]τ]"

  t12 = rootedtree([1, 2, 2, 3, 4])
  @test order(t12) == 5
  @test σ(t12) == 1
  @test γ(t12) == 30
  @test α(t12) == 4
  @test β(t12) == α(t12)*γ(t12)
  @test t12 == t2 ∘ t4 == t8 ∘ t1
  @test butcher_representation(t12) == "[[[τ]]τ]"

  t13 = rootedtree([1, 2, 3, 2, 3])
  @test order(t13) == 5
  @test σ(t13) == 2
  @test γ(t13) == 20
  @test α(t13) == 3
  @test β(t13) == α(t13)*γ(t13)
  @test t13 == t4 ∘ t2
  @test butcher_representation(t13) == "[[τ][τ]]"

  t14 = rootedtree([1, 2, 3, 3, 3])
  @test order(t14) == 5
  @test σ(t14) == 6
  @test γ(t14) == 20
  @test α(t14) == 1
  @test β(t14) == α(t14)*γ(t14)
  @test t14 == t1 ∘ t5
  @test butcher_representation(t14) == "[[τ³]]"

  t15 = rootedtree([1, 2, 3, 3, 4])
  @test order(t15) == 5
  @test σ(t15) == 1
  @test γ(t15) == 40
  @test α(t15) == 3
  @test β(t15) == α(t15)*γ(t15)
  @test t15 == t1 ∘ t6
  @test butcher_representation(t15) == "[[[τ]τ]]"

  t16 = rootedtree([1, 2, 3, 4, 4])
  @test order(t16) == 5
  @test σ(t16) == 2
  @test γ(t16) == 60
  @test α(t16) == 1
  @test β(t16) == α(t16)*γ(t16)
  @test t16 == t1 ∘ t7
  @test butcher_representation(t16) == "[[[τ²]]]"

  t17 = rootedtree([1, 2, 3, 4, 5])
  @test order(t17) == 5
  @test σ(t17) == 1
  @test γ(t17) == 120
  @test α(t17) == 1
  @test β(t17) == α(t17)*γ(t17)
  @test t17 == t1 ∘ t8
  @test butcher_representation(t17) == "[[[[τ]]]]"

  # test non-canonical representation
  level_sequence = [1, 2, 3, 2, 3, 4, 2, 3, 2, 3, 4, 5, 6, 2, 3, 4]
  @test σ(rootedtree(level_sequence)) == σ(RootedTrees.RootedTree(level_sequence))
end

# see butcher2008numerical, Table 302(I)
@testset "number of trees" begin
  number_of_rooted_trees = [1, 1, 2, 4, 9, 20, 48, 115, 286, 719]
  for order in 1:10
    num = 0
    for t in RootedTreeIterator(order)
      num += 1
    end
    @test num == number_of_rooted_trees[order] == count_trees(order)
  end
end

# Runge-Kutta method SSPRK33 of order 3
@testset "Runge-Kutta order conditions" begin
  A = [0 0 0; 1 0 0; 1/4 1/4 0]
  b = [1/6, 1/6, 2/3]
  c = A * fill(1, length(b))
  for order in 1:3
    for t in RootedTreeIterator(order)
      @test residual_order_condition(t, A, b, c) ≈ 0 atol=eps()
    end
  end
  let order=4
    res = 0.0
    for t in RootedTreeIterator(order)
      res += abs(residual_order_condition(t, A, b, c))
    end
    @test res > 10*eps()
  end

  A = @SArray [0 0 0; 1 0 0; 1/4 1/4 0]
  b = @SArray [1/6, 1/6, 2/3]
  c = A * SVector(1, 1, 1)
  for order in 1:3
    @test all(RootedTreeIterator(order)) do t
      abs(residual_order_condition(t, A, b, c)) < eps()
    end
  end

  let order=4
    res = 0.0
    for t in RootedTreeIterator(order)
      res += abs(residual_order_condition(t, A, b, c))
    end
    @test res > 10*eps()
  end
end

# See Section 2.3 and Table 2 of
# - Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
#   Algebraic Structures of B-series.
#   Foundations of Computational Mathematics
#   [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
@testset "partitions" begin
  let t = rootedtree([1, 2, 3, 4, 3])
    edge_set = [true, true, false, false]
    reference_forest = [rootedtree([1, 2, 3]),
                        rootedtree([4]),
                        rootedtree([3])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 2])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 4, 3])
    edge_set = [false, true, true, false]
    reference_forest = [rootedtree([3]),
                        rootedtree([2, 3, 4]),
                        rootedtree([1])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 3])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 4, 3])
    edge_set = [false, true, false, false]
    reference_forest = [rootedtree([4]),
                        rootedtree([3]),
                        rootedtree([2, 3]),
                        rootedtree([1])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 3, 3])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 2, 2, 2])
    edge_set = [false, false, true, true]
    reference_forest = [rootedtree([2]),
                        rootedtree([2]),
                        rootedtree([1, 2, 2])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 2])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 2, 2])
    edge_set = [false, false, false, true]
    reference_forest = [rootedtree([3]),
                        rootedtree([2]),
                        rootedtree([2]),
                        rootedtree([1, 2])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 3, 2])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 2, 2])
    edge_set = [true, true, true, true]
    reference_forest = [rootedtree([1, 2, 3, 2, 2])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 2, 3])
    edge_set = [true, true, false, false]
    reference_forest = [rootedtree([3]),
                        rootedtree([2]),
                        rootedtree([1, 2, 3])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 3])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 2, 3])
    edge_set = [false, true, false, false]
    reference_forest = [rootedtree([2, 3]),
                        rootedtree([3]),
                        rootedtree([2]),
                        rootedtree([1])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 2, 3])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  let t = rootedtree([1, 2, 3, 3, 3])
    edge_set = [false, true, true, false]
    reference_forest = [rootedtree([3]),
                        rootedtree([2, 3, 3]),
                        rootedtree([1])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2, 3])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end

  # additional tests not included in the examples of the paper
  let t = rootedtree([1, 2, 3, 2, 3])
    edge_set = [true, false, true, true]
    reference_forest = [rootedtree([1, 2, 3, 2]),
                        rootedtree([3])]
    @test sort!(partition_forest(t, edge_set)) == sort!(reference_forest)
    @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

    reference_skeleton = rootedtree([1, 2])
    @test reference_skeleton == partition_skeleton(t, edge_set)
  end
end

# See Table 3 of
# - Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
#   Algebraic Structures of B-series.
#   Foundations of Computational Mathematics
#   [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
@testset "all_partitions" begin
  t = rootedtree([1, 2, 3, 3])
  forests, skeletons = all_partitions(t)
  for forest in forests
    for tree in forest
      RootedTrees.normalize_root!(tree)
    end
    sort!(forest)
  end
  sort!(forests)
  for tree in skeletons
    RootedTrees.normalize_root!(tree)
  end
  sort!(skeletons)

  reference_forests = [
    [rootedtree([1, 2, 3, 3]),],
    [rootedtree([1]), rootedtree([1, 2, 2])],
    [rootedtree([1]), rootedtree([1, 2, 3])],
    [rootedtree([1]), rootedtree([1, 2, 3]),],
    [rootedtree([1]), rootedtree([1]), rootedtree([1, 2])],
    [rootedtree([1]), rootedtree([1]), rootedtree([1, 2])],
    [rootedtree([1]), rootedtree([1]), rootedtree([1, 2])],
    [rootedtree([1]), rootedtree([1]), rootedtree([1]), rootedtree([1])],
  ]
  reference_skeletons = [
    rootedtree([1]),
    rootedtree([1, 2]),
    rootedtree([1, 2]),
    rootedtree([1, 2]),
    rootedtree([1, 2, 2]),
    rootedtree([1, 2, 3]),
    rootedtree([1, 2, 3]),
    rootedtree([1, 2, 3, 3]),
  ]
  for forest in reference_forests
    sort!(forest)
  end
  sort!(reference_forests)
  sort!(reference_skeletons)

  @test forests == reference_forests
  @test skeletons == reference_skeletons

  @testset "PartitionIterator" begin
    partitions = collect(PartitionIterator(t))
    iterator_forests = map(first, partitions)
    iterator_skeletons = map(last, partitions)
    for forest in iterator_forests
      sort!(forest)
    end
    sort!(iterator_forests)
    sort!(iterator_skeletons)
    @test iterator_forests == forests
    @test iterator_skeletons == skeletons

    for order in 1:8
      for t in RootedTreeIterator(order)
        forests, skeletons = all_partitions(t)
        @test collect(zip(forests, skeletons)) == collect(PartitionIterator(t))
      end
    end

    level_sequence = zeros(Int, RootedTrees.BUFFER_LENGTH + 1)
    level_sequence[1] -= 1
    t = rootedtree(level_sequence)
    @inferred PartitionIterator(t)
    t = @inferred rootedtree!(view(level_sequence, :))
    @inferred PartitionIterator(t)
  end
end

# See Section 2.2 and Table 1 of
# - Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
#   Algebraic Structures of B-series.
#   Foundations of Computational Mathematics
#   [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
@testset "splittings" begin
  t = rootedtree([1, 2, 3, 2, 2])
  splittings = all_splittings(t)
  forests_and_subtrees = sort!(collect(zip(splittings.forests, splittings.subtrees)))

  reference_forests_and_subtrees = sort!([
    (empty([rootedtree([1])]),                               rootedtree([1, 2, 3, 2, 2])),
    ([rootedtree([1])],                                      rootedtree([1, 2, 3, 2])),
    ([rootedtree([1])],                                      rootedtree([1, 2, 3, 2])),
    ([rootedtree([1, 2])],                                   rootedtree([1, 2, 2])),
    ([rootedtree([1])],                                      rootedtree([1, 2, 2, 2])),
    ([rootedtree([1, 2]), rootedtree([1])],                  rootedtree([1, 2])),
    ([rootedtree([1, 2]), rootedtree([1])],                  rootedtree([1, 2])),
    ([rootedtree([1]),    rootedtree([1])],                  rootedtree([1, 2, 3])),
    ([rootedtree([1]),    rootedtree([1])],                  rootedtree([1, 2, 2])),
    ([rootedtree([1]),    rootedtree([1])],                  rootedtree([1, 2, 2])),
    ([rootedtree([1]),    rootedtree([1]), rootedtree([1])], rootedtree([1, 2])),
    ([rootedtree([1, 2]), rootedtree([1]), rootedtree([1])], rootedtree([1])),
    ([rootedtree([1, 2, 3, 2, 2])],                          rootedtree(Int[])),
  ])

  @test forests_and_subtrees == reference_forests_and_subtrees

  # tested with the Python package BSeries
  t = rootedtree([1, 2, 3, 4, 4, 2, 3, 3, 2, 3, 2, 3])
  splittings = all_splittings(t)
  @test length(splittings.forests) == length(splittings.subtrees) == 271

  # consistency of all_splittings and the SplittingIterator
  for order in 1:8
    for i in RootedTreeIterator(order)
      forests, subtrees = all_splittings(t)
      @test collect(zip(forests, subtrees)) == collect(SplittingIterator(t))
    end
  end
end

@testset "plots" begin
  for order in 1:4
    for t in RootedTreeIterator(order)
      plot(t)
    end
  end
end

end # @testset "RootedTrees"


text:
"""
Force (re-)evaluation of the objective value at `x`.

Returns `f(x)` and stores the value in `obj.F`
"""
function value!!(obj::AbstractObjective, x)
    obj.f_calls .+= 1
    copyto!(obj.x_f, x)
    obj.F = obj.f(x)
    value(obj)
end
"""
Evaluates the objective value at `x`.

Returns `f(x)`, but does *not* store the value in `obj.F`
"""
function value(obj::AbstractObjective, x)
    obj.f_calls .+= 1
    return obj.f(x)
end
"""
Evaluates the objective value at `x`.

Returns `f(x)` and stores the value in `obj.F`
"""
function value!(obj::AbstractObjective, x)
    if x != obj.x_f
        value!!(obj, x)
    end
    value(obj)
end

"""
Evaluates the gradient value at `x`

This does *not* update `obj.DF` or `obj.x_df`.
"""
function gradient(obj::AbstractObjective, x)
    newdf = copy(obj.DF)
    obj.df(newdf, x)
    obj.df_calls .+= 1
    return newdf
end
"""
Evaluates the gradient value at `x`.

Stores the value in `obj.DF`.
"""
function gradient!(obj::AbstractObjective, x)
    if x != obj.x_df
        gradient!!(obj, x)
    end
    gradient(obj)
end
"""
Force (re-)evaluation of the gradient value at `x`.

Stores the value in `obj.DF`.
"""
function gradient!!(obj::AbstractObjective, x)
    obj.df_calls .+= 1
    copyto!(obj.x_df, x)
    obj.df(obj.DF, x)
    gradient(obj)
end

function value_gradient!(obj::AbstractObjective, x)
    if x != obj.x_f && x != obj.x_df
        value_gradient!!(obj, x)
    elseif x != obj.x_f
        value!!(obj, x)
    elseif x != obj.x_df
        gradient!!(obj, x)
    end
    value(obj), gradient(obj)
end
function value_gradient!!(obj::AbstractObjective, x)
    obj.f_calls .+= 1
    obj.df_calls .+= 1
    copyto!(obj.x_f, x)
    copyto!(obj.x_df, x)
    obj.F = obj.fdf(gradient(obj), x)
    value(obj), gradient(obj)
end

function hessian!(obj::AbstractObjective, x)
    if x != obj.x_h
        hessian!!(obj, x)
    end
    hessian(obj)
end
function hessian!!(obj::AbstractObjective, x)
    obj.h_calls .+= 1
    copyto!(obj.x_h, x)
    obj.h(obj.H, x)
    hessian(obj)
end

# Getters are without ! and accept only an objective and index or just an objective
"Get the most recently evaluated objective value of `obj`."
value(obj::AbstractObjective) = obj.F
"Get the most recently evaluated gradient of `obj`."
gradient(obj::AbstractObjective) = obj.DF
"Get the most recently evaluated Jacobian of `obj`."
jacobian(obj::AbstractObjective) = obj.DF
"Get the `i`th element of the most recently evaluated gradient of `obj`."
gradient(obj::AbstractObjective, i::Integer) = obj.DF[i]
"Get the most recently evaluated Hessian of `obj`"
hessian(obj::AbstractObjective) = obj.H

value_jacobian!(obj, x) = value_jacobian!(obj, obj.F, obj.DF, x)
function value_jacobian!(obj, F, J, x)
    if x != obj.x_f && x != obj.x_df
        value_jacobian!!(obj, F, J, x)
    elseif x != obj.x_f
        value!!(obj, F, x)
    elseif x != obj.x_df
        jacobian!!(obj, J, x)
    end
    F, J
end
value_jacobian!!(obj, x) = value_jacobian!!(obj, obj.F, obj.DF, x)
function value_jacobian!!(obj, F, J, x)
    obj.fdf(F, J, x)
    copyto!(obj.x_f, x)
    copyto!(obj.x_df, x)
    obj.f_calls .+= 1
    obj.df_calls .+= 1
    obj.df_calls
    F, J
end

function jacobian!(obj, x)
    if x != obj.x_df
        jacobian!!(obj, x)
    end
    jacobian(obj)
end

jacobian!!(obj, x) = jacobian!!(obj, obj.DF, x)
function jacobian!!(obj, J, x)
    obj.df(J, x)
    copyto!(obj.x_df, x)
    obj.df_calls .+= 1
    obj.df_calls
    J
end
function jacobian(obj::AbstractObjective, x)
    tmp = copy(obj.DF)
    jacobian!!(obj, x)
    newdf = copy(obj.DF)
    copyto!(obj.DF, tmp)
    return newdf
end

value(obj::NonDifferentiable{TF, TX}, x) where {TF<:AbstractArray, TX} = value(obj, copy(obj.F), x)
value(obj::OnceDifferentiable{TF, TDF, TX}, x) where {TF<:AbstractArray, TDF, TX} = value(obj, copy(obj.F), x)
function value(obj::AbstractObjective, F, x)
    obj.f_calls .+= 1
    return obj.f(F, x)
end

value!!(obj::NonDifferentiable{TF, TX}, x) where {TF<:AbstractArray, TX} = value!!(obj, obj.F, x)
value!!(obj::OnceDifferentiable{TF, TDF, TX}, x) where {TF<:AbstractArray, TDF, TX} = value!!(obj, obj.F, x)
function value!!(obj::AbstractObjective, F, x)
    obj.f(F, x)
    copyto!(obj.x_f, x)
    obj.f_calls .+= 1
    obj.f_calls
    F
end

function _clear_f!(d::NLSolversBase.AbstractObjective)
    d.f_calls .= 0
    if typeof(d.F) <: AbstractArray
        d.F .= eltype(d.F)(NaN)
    else
        d.F = typeof(d.F)(NaN)
    end
    d.x_f .= eltype(d.x_f)(NaN)
    nothing
end

function _clear_df!(d::NLSolversBase.AbstractObjective)
    d.df_calls .= 0
    d.DF .= eltype(d.DF)(NaN)
    d.x_df .= eltype(d.x_df)(NaN)
    nothing
end

function _clear_h!(d::NLSolversBase.AbstractObjective)
    d.h_calls .= 0
    d.H .= eltype(d.H)(NaN)
    d.x_h .= eltype(d.x_h)(NaN)
    nothing
end

function _clear_hv!(d::NLSolversBase.AbstractObjective)
    d.hv_calls .= 0
    d.Hv .= eltype(d.Hv)(NaN)
    d.x_hv .= eltype(d.x_hv)(NaN)
    d.v_hv .= eltype(d.v_hv)(NaN)
    nothing
end

clear!(d::NonDifferentiable)  = _clear_f!(d)

function clear!(d::OnceDifferentiable)
    _clear_f!(d)
    _clear_df!(d)
    nothing
end

function clear!(d::TwiceDifferentiable)
    _clear_f!(d)
    _clear_df!(d)
    _clear_h!(d)
    nothing
end

function clear!(d::TwiceDifferentiableHV)
    _clear_f!(d)
    _clear_df!(d)
    _clear_hv!(d)
    nothing
end

g_calls(d::NonDifferentiable) = 0
h_calls(d::Union{NonDifferentiable, OnceDifferentiable}) = 0
f_calls(d) = first(d.f_calls)
g_calls(d) = first(d.df_calls)
h_calls(d) = first(d.h_calls)
hv_calls(d) = 0
h_calls(d::TwiceDifferentiableHV) = 0
hv_calls(d::TwiceDifferentiableHV) = first(d.hv_calls)


text:
using Distributed

function main()
    if length(ARGS) != 2 || nprocs() == 1
        println("Usage: reduce -p n <# iterations> <vector_length>")
        exit(1)
    end

    argv = map(x -> parse(Int, x), ARGS)
    iterations, payload = argv

    ranks = nworkers()

    @sync for proc in workers()
        @spawnat proc begin
            global inbox = RemoteChannel(() -> Channel{Vector{Float64}}(0))
        end
    end

    futures = Vector{Future}(undef, ranks)
    @sync for proc in workers()
        futures[proc-1] = @spawnat proc begin
            reduce_channels = [(@fetchfrom proc get_channel()) for proc in sort(workers())]
            work(iterations, payload, inbox, reduce_channels)
        end
    end

    # Fetch results
    results = map(fetch, futures)
    value = results[1][1]
    time = maximum(map(x -> x[2], results))

    ground_truth = iterations+2.0+(iterations*iterations+5.0*iterations+4.0)*(ranks-1)/2;
    if abs(value - ground_truth) <= 1.e-8
        println("Solution validates")
    else
        println("Value provided: $value, expected: $ground_truth")
    end
    avgtime = time / iterations
    throughput = 1.e-6 * (2.0*ranks-1.0)*payload/avgtime
    println("Rate (MFlops/s): $throughput Avg time (s): $avgtime")
end

@everywhere function get_channel()
    return inbox
end

@everywhere function work(iterations, payload, local_channel, channels)
    local_sum = fill(1.0, payload)
    constant_vector = fill(1.0, payload)

    # Setup for tree topology
    n = nworkers()
    pow_2 = 2^Int(ceil(log(2,n))) # Smallest power of two ≤
    lowest_z = pow_2 - nworkers() + 1

    local t0, t1
    for it in 0:iterations
        if it == 1
            t0 = time_ns()
        end
        local_sum .+= constant_vector

        # Assign the accumulator
        if myid() == 2
            accumulator = local_sum
        else
            accumulator = copy(local_sum)
        end
        id_index = myid() - 1
        # z is the relative position variable in the reduction
        pos = z = pow_2 - myid() + 2
        k = 0 # offset for who to collect from next
        while z % 2 == 0
            if pos - 2^k >= lowest_z
                c = channels[id_index + 2^k]
                incoming = take!(c)
                accumulator .+= incoming
            end
            z ÷= 2
            k += 1
        end

        if myid() != 2
            put!(local_channel, accumulator)
        end
    end

    t1 = time_ns()
    time = (t1-t0)*1.e-9

    if myid() == 2
        return (local_sum[1], time)
    else
        return (nothing, time)
    end
end

main()


text:
using Test, Lasso, LassoPlot
using DataFrames, CSV, Pkg, Random, Plots

include("plots.jl")


text:
module ParameterJuMP

using SparseArrays

using JuMP
export index
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

export
ModelWithParams, ParameterRef, add_parameter, add_parameters, all_parameters, Param, parametrized_dual_objective_value

# types
# ------------------------------------------------------------------------------
struct ParameterRef <: JuMP.AbstractVariableRef
    ind::Int64 # local reference
    model::JuMP.Model
end

# Reference to a constraint in which the parameter has coefficient coef
struct ParametrizedConstraintRef{C}
    cref::JuMP.ConstraintRef{JuMP.Model, C}
    coef::Float64
end

const CtrRef{F, S} = ConstraintRef{JuMP.Model,MOI.ConstraintIndex{F,S}, JuMP.ScalarShape}
const SAF = MOI.ScalarAffineFunction{Float64}
const EQ = MOI.EqualTo{Float64}
const LE = MOI.LessThan{Float64}
const GE = MOI.GreaterThan{Float64}

mutable struct ParameterData
    inds::Vector{Int64}
    next_ind::Int64

    current_values::Vector{Float64}
    future_values::Vector{Float64}

    # Whether the JuMP model is in sync with the value of the parameters
    sync::Bool

    # constraints where it is a RHS
    constraints_map::Dict{Int64, Vector{ParametrizedConstraintRef}}

    parameters_map_saf_in_eq::Dict{CtrRef{SAF, EQ}, JuMP.GenericAffExpr{Float64,ParameterRef}}
    parameters_map_saf_in_le::Dict{CtrRef{SAF, LE}, JuMP.GenericAffExpr{Float64,ParameterRef}}
    parameters_map_saf_in_ge::Dict{CtrRef{SAF, GE}, JuMP.GenericAffExpr{Float64,ParameterRef}}

    names::Dict{ParameterRef, String}

    # overload addvariable
    # variables_map::Dict{Int64, Vector{Int64}}
    # variables_map_lb::Dict{Int64, Vector{Int64}}
    # variables_map_ub::Dict{Int64, Vector{Int64}}

    # solved::Bool

    has_deletion::Bool
    index_map::Dict{Int64, Int64}

    lazy::Bool

    no_duals::Bool
    dual_values::Vector{Float64}
    function ParameterData()
        new(Int64[],
            1,
            Float64[],
            Float64[],
            true,
            Dict{Int64, Vector{JuMP.ConstraintRef}}(),
            Dict{CtrRef{SAF, EQ}, JuMP.GenericAffExpr{Float64,ParameterRef}}(),
            Dict{CtrRef{SAF, LE}, JuMP.GenericAffExpr{Float64,ParameterRef}}(),
            Dict{CtrRef{SAF, GE}, JuMP.GenericAffExpr{Float64,ParameterRef}}(),
            Dict{ParameterRef, String}(),
            false,
            Dict{Int64, Int64}(),
            false,
            false,
            Float64[]
            )
    end
end

no_duals(data::ParameterData) = data.no_duals
no_duals(model::JuMP.Model) = no_duals(_getparamdata(model))

set_no_duals(model::JuMP.Model) = set_no_duals(_getparamdata(model))
function set_no_duals(data::ParameterData)
    if isempty(data.current_values)
        data.no_duals = true
    elseif no_duals(data)
        @warn "No duals mode is already activated"
    else
        error("Parameter JuMP's no duals mode can only be activated in empty models.")
    end
end

"""
    JuMP.index(p::ParameterRef)::Int64

Return the internal index of the parameter `p`.
"""
JuMP.index(p::ParameterRef) = index(_getparamdata(p.model), p)::Int64
JuMP.index(model::JuMP.Model, p::ParameterRef) = index(_getparamdata(model), p)::Int64
function JuMP.index(data::ParameterData, p::ParameterRef)::Int64
    if data.has_deletion
        return data.index_map[p.ind]
    else
        return p.ind
    end
end

lazy_duals(data::ParameterData) = data.lazy
lazy_duals(model::JuMP.Model) = lazy_duals(_getparamdata(model))

set_lazy_duals(model::JuMP.Model) = set_lazy_duals(_getparamdata(model))
function set_lazy_duals(data::ParameterData)
    if isempty(data.current_values)
        data.lazy = true
    elseif lazy_duals(data)
        @warn "Lazy mode is already activated"
    else
        error("Parameter JuMP's lazy mode can only be activated in empty models.")
    end
end

set_not_lazy_duals(model::JuMP.Model) = set_not_lazy_duals(_getparamdata(model))
function set_not_lazy_duals(data::ParameterData)
    if isempty(data.current_values)
        data.lazy = false
    elseif !lazy_duals(data)
        @warn "Lazy mode is already de-activated"
    else
        error("Parameter JuMP's lazy mode can only be de-activated in empty models.")
    end
end

_get_param_dict(data::ParameterData, ::Type{EQ}) = data.parameters_map_saf_in_eq
_get_param_dict(data::ParameterData, ::Type{LE}) = data.parameters_map_saf_in_le
_get_param_dict(data::ParameterData, ::Type{GE}) = data.parameters_map_saf_in_ge

_getmodel(p::ParameterRef) = p.model
_getparamdata(p::ParameterRef)::ParameterData = _getparamdata(_getmodel(p))::ParameterData
# _getparamdata(model::JuMP.Model) = model.ext[:params]::ParameterData
function _getparamdata(model::JuMP.Model)::ParameterData
    # TODO checks dict twice
    !haskey(model.ext, :params) && error("In order to use Parameters the model must be created with the ModelWithParams constructor")
    return model.ext[:params]
end

"""
    is_sync(model::JuMP.Model)

Test if the JuMP Model is updated to the latest value of all parameters.
"""
is_sync(data::ParameterData) = data.sync
is_sync(model::JuMP.Model) = is_sync(_getparamdata(model))

"""
    add_parameter(model::JuMP.Model, val::Real)::ParameterRef

Adds one parameter fixed at `val` to the model `model`.

    add_parameter(model::JuMP.Model)::ParameterRef

Adds one parameter fixed at zero to the model `model`.
"""
function add_parameter(model::JuMP.Model, val::Real)
    params = _getparamdata(model)::ParameterData

    ind = params.next_ind
    params.next_ind += 1

    push!(params.inds, ind)
    push!(params.current_values, 0.0)
    push!(params.future_values, val)
    if !iszero(val)
        params.sync = false
    end
    push!(params.dual_values, NaN)

    params.constraints_map[ind] = ParametrizedConstraintRef[]

    return ParameterRef(ind, model)
end
add_parameter(model) = add_parameter(model, 0.0)

_addsizehint!(vec::Vector, len::Integer) = sizehint!(vec, len+length(vec))

"""
    add_parameters(model::JuMP.Model, val::Vector{R})::Vector{ParameterRef}

Adds one parameter for each element of the vector `val`.

    add_parameters(model::JuMP.Model, N::Integer)::Vector{ParameterRef}

Adds `N` parameters to the model, all of them fixed in zero.
"""
function add_parameters(model::JuMP.Model, N::Integer)
    return add_parameters(model::JuMP.Model, zeros(N))
end
function add_parameters(model::JuMP.Model, val::AbstractArray{R,N}) where {R,N}
    params = _getparamdata(model)::ParameterData

    nparam = length(val)
    out = ParameterRef[]
    out = similar(val, ParameterRef)
    sizehint!(out, nparam)
    _addsizehint!(params.inds, nparam)
    _addsizehint!(params.current_values, nparam)
    _addsizehint!(params.future_values, nparam)
    _addsizehint!(params.dual_values, nparam)

    for i in eachindex(val)
        ind = params.next_ind
        params.next_ind += 1

        push!(params.inds, ind)
        push!(params.current_values, 0.0)
        push!(params.future_values, val[i])
        if !iszero(val[i])
            params.sync = false
        end
        push!(params.dual_values, NaN)

        params.constraints_map[ind] = ParametrizedConstraintRef[]

        out[i] = ParameterRef(ind, model)
    end

    return out
end

"""
    all_parameters(model::JuMP.AbstractModel)::Vector{ParameterRef}

Returns a list of all parameters currently in the model. The parameters are
ordered by creation time.
"""
function all_parameters(model::JuMP.AbstractModel)
    data = _getparamdata(model)
    return ParameterRef[ParameterRef(ind, model) for ind in data.inds]
end

# solve
# ------------------------------------------------------------------------------

"""
    ModelWithParams(args...; kwargs...)::JuMP.Model

Returns a JuMP.Model able to handle `Parameters`.

`args` and `kwargs` are the same parameters that would be passed
to the regular `Model` constructor.

Example using GLPK solver:

```julia
    model = ModelWithParams(with_optimizer(GLPK.Optimizer))
```
"""
function ModelWithParams(args...; kwargs...)
    m = JuMP.Model(args...; kwargs...)
    enable_parameters(m)
    return m
end

"""
    enable_parameters(m::JuMP.Model)::Nothing

Enables a `JuMP.Model` to handle `Parameters`.
"""
function enable_parameters(m::JuMP.Model)
    if haskey(m.ext, :params)
        error("Model already has parameter enabled")
    end
    # test and compare hook
    initialize_parameter_data(m)
    JuMP.set_optimize_hook(m, parameter_optimizehook)
    return nothing
end

function initialize_parameter_data(m::JuMP.Model)
    m.ext[:params] = ParameterData()
end

# The constraint has been deleted. We cannot keep `dict` in sync with deletions
# as we return a `JuMP.ConstraintRef`, not a custom type when the user create a
# parametrized constraint.
function _update_dicts_with_deletion(data::ParameterData, S::Type)
    dict = _get_param_dict(data, S)
    to_delete = eltype(keys(dict))[]
    for (cref, gaep) in dict
        if !JuMP.is_valid(cref.model, cref)
            push!(to_delete, cref)
        end
    end
    for cref in to_delete
        delete!(dict, cref)
    end
end
function _update_dicts_with_deletion(data::ParameterData)
    _update_dicts_with_deletion(data, EQ)
    _update_dicts_with_deletion(data, LE)
    _update_dicts_with_deletion(data, GE)
end
function parameter_optimizehook(m::JuMP.Model; kwargs...)
    data = _getparamdata(m)::ParameterData

    _update_dicts_with_deletion(data)

    # sync model rhs to newest parameter values
    sync(data)

    ret = JuMP.optimize!(m::JuMP.Model, ignore_optimize_hook = true, kwargs...)

    # update duals for later query
    if !no_duals(data) && JuMP.has_duals(m)
        _update_duals(data)
    end
    return ret
end

"""
    sync(model::JuMP.Model)::Nothing

Forces the model to update its constraints to the new values of the
Parameters.
"""
function sync(data::ParameterData)
    if !is_sync(data)
        _update_constraints(data)
        data.current_values .= data.future_values
        data.sync = true
    end
    return nothing
end
function sync(model::JuMP.Model)
    sync(_getparamdata(model))
end

function parametrized_dual_objective_value(model::JuMP.AbstractModel)
    params = all_parameters(model)
    linear = [param => dual(param) for param in params]
    obj = dual_objective_value(model)
    constant = obj - sum(value(term.first) * term.second for term in linear)
    return GenericAffExpr(constant, linear)
end

# constraints
# ------------------------------------------------------------------------------
include("constraints.jl")

# JuMP variable interface
# ------------------------------------------------------------------------------
include("variable_interface.jl")

# operators and mutable arithmetics
# ------------------------------------------------------------------------------
include("operators.jl")
include("mutable_arithmetics.jl")

# other
# ------------------------------------------------------------------------------
include("macros.jl")
include("print.jl")

end


text:
a=nothing
if a
    println("cool")
else
    println("looc")
end


text:
export device, GPU, CPU

abstract type AbstractDevice end
abstract type CPU <: AbstractDevice end
abstract type GPU{N} <: AbstractDevice end

device(x::T) where T = device(T)
device(::Type{T}) where T = CPU

import Base: show

show(io::IO, d::Type{CPU}) = print(io, "CPU")
show(io::IO, d::Type{GPU{N}}) where N = print(io, "GPU{$N}")


text:
using NIFTy
using Documenter

makedocs(;
    modules=[NIFTy],
    authors="Philipp Haim <philipp.haim@gmail.com> and contributors",
    repo="https://github.com/Philipp Haim/NIFTy.jl/blob/{commit}{path}#L{line}",
    sitename="NIFTy.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Philipp Haim.gitlab.io/NIFTy.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)