text:
using Test
using StaticArrays
using RootedTrees

using RootedTrees.Latexify: latexify

using Plots: Plots, plot
Plots.unicodeplots()

@testset "RootedTrees" begin

@testset "RootedTree" 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 butcher_representation(empty(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

  # 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

end # @testset "RootedTree"

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

    for (t1,t2,t3,t4,t5,t6,t7,t8,t9) 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 !(t1 == t6)
      @test !(t1 == t7)
      @test !(t1 == t8)

      @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 hash(t1) != hash(t6)
      @test hash(t1) != hash(t7)
      @test hash(t1) != hash(t8)

      @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)
      @test t1 < t6
      @test t1 < t7
      @test t1 < t8
      @test t6 < t7
      @test t6 < t8
      @test t7 < t8
      @test !(t1 < t9)
      @test !(t6 < t9)
      @test !(t7 < t9)
      @test !(t8 < t9)

      println(devnull, t1)
      println(devnull, t2)
      println(devnull, t3)
      println(devnull, t4)
      println(devnull, t5)
      println(devnull, t6)
      println(devnull, t7)
      println(devnull, t8)
    end

    @test rootedtree([1, 2]            ) > rootedtree([1]         )
    @test rootedtree([1, 2], Bool[0, 1]) > rootedtree([1], Bool[1])

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

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

    t = rootedtree([1, 2, 3, 2, 3, 4, 2, 3], [1, 0, 1, 1, 0, 0, 0, 1])
    @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], Bool[1, 0, 1, 1, 0, 0, 0, 1])
    @test t.level_sequence == [1, 2, 3, 3, 2, 3, 2, 3]
    @test !isempty(t)

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

    # misc
    @test_throws DimensionMismatch rootedtree([1, 2, 3], [1, 2])
  end

  @testset "hashing" begin
    hashes = [hash(rootedtree(Int[], Bool[]))]
    for o in 1:8
      for t in BicoloredRootedTreeIterator(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],
                   Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0])
    new_hash = @inferred hash(t)
    @test !(new_hash in hashes)
  end

  @testset "functions on trees" begin
    # See Araujo, Murua, and Sanz-Serna (1997), Table 1
    # https://doi.org/10.1137/S0036142995292128
    let t = rootedtree(Int[], Int[])
      @test order(t) == 0
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 1
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "∅"
    end

    let t = rootedtree([1], [1])
      @test order(t) == 1
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 1
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "τ₁"

      @test t ∘ t == rootedtree([1, 2], [1, 1])
    end

    let t = rootedtree([1], [2])
      @test order(t) == 1
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 1
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "τ₂"
    end

    let t = rootedtree([1], [3])
      @test order(t) == 1
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 1
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "τ₃"
    end

    let t = rootedtree([1], [typemax(Int)])
      @test order(t) == 1
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 1
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "τ₀"
    end

    let t = rootedtree([1, 2], [1, 1])
      @test order(t) == 2
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 2
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "[τ₁]₁"
    end

    let t = rootedtree([1, 2], [1, 2])
      @test order(t) == 2
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 2
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "[τ₂]₁"
    end

    let t = rootedtree([1, 2], [3, 1])
      @test order(t) == 2
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 2
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "[τ₁]₃"
    end

    let t = rootedtree([1, 2, 2], [2, 1, 1])
      @test order(t) == 3
      @test α(t) == 1
      @test σ(t) == 2
      @test γ(t) == 3
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "[τ₁²]₂"
    end

    let t = rootedtree([1, 2, 2], [2, 1, 2])
      @test order(t) == 3
      @test α(t) == 2
      @test σ(t) == 1
      @test γ(t) == 3
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "[τ₂τ₁]₂"
    end

    let t = rootedtree([1, 2, 3], [3, 2, 1])
      @test order(t) == 3
      @test α(t) == 1
      @test σ(t) == 1
      @test γ(t) == 6
      @test_nowarn println(devnull, t)
      @test butcher_representation(t) == "[[τ₁]₂]₃"
    end
  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 BicoloredRootedTreeIterator(order)
        num += 1
      end
      # number of plain rooted trees times number of possible color sequences
      # <= since not all possible color sequences are in canonical representation
      @test num <= number_of_rooted_trees[order] * 2^order
    end
  end

  # https://github.com/SciML/RootedTrees.jl/issues/72
  @testset "BicoloredRootedTreeIterator is canonical" begin
    for o in 1:10
      for t_iterator in BicoloredRootedTreeIterator(o)
        t_canonical = RootedTrees.canonical_representation(t_iterator)
        @test t_iterator == t_canonical
      end
    end
  end

  # See Sections 2.3 & 6.1 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
    # Example in Section 6.1
    let t = rootedtree([1, 2, 3, 3], Bool[0, 1, 0, 0])
      edge_set = [false, true, false]
      reference_forest = [rootedtree([3], Bool[0]),
                          rootedtree([2, 3], Bool[1, 0]),
                          rootedtree([1], Bool[0])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    # Other examples for single-colored trees
    let t = rootedtree([1, 2, 3, 4, 3], Bool[1, 1, 0, 1, 1])
      edge_set = [true, true, false, false]
      reference_forest = [rootedtree([1, 2, 3], Bool[1, 1, 0]),
                          rootedtree([4], Bool[1]),
                          rootedtree([3], Bool[1])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 4, 3], rand(Bool, 5))
      edge_set = [false, true, true, false]
      reference_forest = [rootedtree([3], t.color_sequence[5:5]),
                          rootedtree([2, 3, 4], t.color_sequence[2:4]),
                          rootedtree([1], t.color_sequence[1:1])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 4, 3], rand(Bool, 5))
      edge_set = [false, true, false, false]
      reference_forest = [rootedtree([4], t.color_sequence[4:4]),
                          rootedtree([3], t.color_sequence[5:5]),
                          rootedtree([2, 3], t.color_sequence[2:3]),
                          rootedtree([1], t.color_sequence[1:1])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 2, 2, 2], rand(Bool, 5))
      edge_set = [false, false, true, true]
      reference_forest = [rootedtree([2], t.color_sequence[2:2]),
                          rootedtree([2], t.color_sequence[3:3]),
                          rootedtree([1, 2, 2], t.color_sequence[[1,4,5]])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 2, 2], rand(Bool, 5))
      edge_set = [false, false, false, true]
      reference_forest = [rootedtree([3], t.color_sequence[3:3]),
                          rootedtree([2], t.color_sequence[2:2]),
                          rootedtree([2], t.color_sequence[4:4]),
                          rootedtree([1, 2], t.color_sequence[[1, 5]])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 2, 2], rand(Bool, 5))
      edge_set = [true, true, true, true]
      reference_forest = [rootedtree([1, 2, 3, 2, 2], t.color_sequence[:])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 2, 3], rand(Bool, 5))
      edge_set = [true, true, false, false]
      reference_forest = [rootedtree([3], t.color_sequence[5:5]),
                          rootedtree([2], t.color_sequence[4:4]),
                          rootedtree([1, 2, 3], t.color_sequence[1:3])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 2, 3], rand(Bool, 5))
      edge_set = [false, true, false, false]
      reference_forest = [rootedtree([2, 3], t.color_sequence[2:3]),
                          rootedtree([3], t.color_sequence[5:5]),
                          rootedtree([2], t.color_sequence[4:4]),
                          rootedtree([1], t.color_sequence[1:1])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

    let t = rootedtree([1, 2, 3, 3, 3], rand(Bool, 5))
      edge_set = [false, true, true, false]
      reference_forest = [rootedtree([3], t.color_sequence[5:5]),
                          rootedtree([2, 3, 3], t.color_sequence[2:4]),
                          rootedtree([1], t.color_sequence[1:1])]
      sort!(reference_forest)
      @test sort!(collect(PartitionForestIterator(t, edge_set))) == reference_forest

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

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

      reference_skeleton = rootedtree([1, 2])
      @test reference_skeleton.level_sequence == partition_skeleton(t, edge_set).level_sequence
    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 "PartitionIterator" begin
    t = rootedtree([1, 2, 3, 3], Bool[1, 0, 1, 0])
    partitions = collect(PartitionIterator(t))
    forests = map(first, partitions)
    skeletons = map(last, partitions)
    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], Bool[1, 0, 1, 0]),],
      [rootedtree([1], Bool[1]), rootedtree([1, 2, 2], Bool[0, 1, 0])],
      [rootedtree([1], Bool[1]), rootedtree([1, 2, 3], Bool[1, 0, 0])],
      [rootedtree([1], Bool[0]), rootedtree([1, 2, 3], Bool[1, 0, 1]),],
      [rootedtree([1], Bool[1]), rootedtree([1], Bool[1]), rootedtree([1, 2], Bool[0, 0])],
      [rootedtree([1], Bool[0]), rootedtree([1], Bool[1]), rootedtree([1, 2], Bool[1, 0])],
      [rootedtree([1], Bool[0]), rootedtree([1], Bool[1]), rootedtree([1, 2], Bool[0, 1])],
      [rootedtree([1], Bool[0]), rootedtree([1], Bool[0]), rootedtree([1], Bool[1]), rootedtree([1], Bool[1])],
    ]
    reference_skeletons = [
      rootedtree([1], Bool[1]),
      rootedtree([1, 2], Bool[1, 0]),
      rootedtree([1, 2], Bool[1, 0]),
      rootedtree([1, 2], Bool[1, 1]),
      rootedtree([1, 2, 2], Bool[1, 1, 0]),
      rootedtree([1, 2, 3], Bool[1, 0, 0]),
      rootedtree([1, 2, 3], Bool[1, 0, 1]),
      rootedtree([1, 2, 3, 3], Bool[1, 0, 1, 0]),
    ]
    for forest in reference_forests
      sort!(forest)
    end
    sort!(reference_forests)
    sort!(reference_skeletons)

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

    level_sequence = zeros(Int, RootedTrees.BUFFER_LENGTH + 1)
    level_sequence[1] -= 1
    color_sequence = rand(Bool, length(level_sequence))
    t = rootedtree(level_sequence, color_sequence)
    @inferred PartitionIterator(t)
    t = @inferred rootedtree!(view(level_sequence, :), view(color_sequence, :))
    @inferred PartitionIterator(t)
  end
end # @testset "ColoredRootedTree"

@testset "Order conditions" begin
  # Runge-Kutta method SSPRK33 of order 3
  @testset "RungeKuttaMethod, SSPRK33" begin
    A = [0 0 0; 1 0 0; 1/4 1/4 0]
    b = [1/6, 1/6, 2/3]
    rk = RungeKuttaMethod(A, b)
    show(IOContext(stdout, :compact=>false), rk)

    for order in 1:3
      for t in RootedTreeIterator(order)
        @test residual_order_condition(t, rk) ≈ 0 atol=eps()
      end
    end

    let order=4
      res = 0.0
      for t in RootedTreeIterator(order)
        res += abs(residual_order_condition(t, rk))
      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]
    rk = RungeKuttaMethod(A, b)
    show(IOContext(stdout, :compact=>true), rk)
    for order in 1:3
      @test all(RootedTreeIterator(order)) do t
        abs(residual_order_condition(t, rk)) < eps()
      end
    end

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

    # deprecations
    let order=4
      for t in RootedTreeIterator(order)
        @test elementary_weight(t, rk.A, rk.b, rk.c) ≈ elementary_weight(t, rk)
        @test derivative_weight(t, rk.A, rk.b, rk.c) ≈ derivative_weight(t, rk)
        @test residual_order_condition(t, rk.A, rk.b, rk.c) ≈ residual_order_condition(t, rk)
      end
    end
  end

  @testset "AdditiveRungeKuttaMethod, IMEX Euler" begin
    ex_euler = RungeKuttaMethod(
      @SMatrix([0//1]), @SVector [1]
    )
    im_euler = RungeKuttaMethod(
      @SMatrix([1//1]), @SVector [1]
    )
    ark = AdditiveRungeKuttaMethod([ex_euler, im_euler])
    show(IOContext(stdout, :compact=>true), ark)
    show(IOContext(stdout, :compact=>false), ark)

    @test elementary_weight(rootedtree(Int[], Bool[]), ark) ≈ 1

    for order in 1:1
      for t in BicoloredRootedTreeIterator(order)
        @test residual_order_condition(t, ark) ≈ 0 atol=eps()
      end
    end

    let order=2
      res = 0.0
      for t in BicoloredRootedTreeIterator(order)
        res += abs(residual_order_condition(t, ark))
      end
      @test res > 10*eps()
    end
  end

  @testset "AdditiveRungeKuttaMethod, Störmer-Verlet" begin
    # Hairer, Lubich, Wanner (2002)
    # Geometric numerical integration
    # Table II.2.1
    As = [
      [0 0; 1//2 1//2],
      [1//2 0; 1//2 0]
    ]
    bs = [
      [1//2, 1//2],
      [1//2, 1//2]
    ]
    ark = AdditiveRungeKuttaMethod(As, bs)

    for order in 1:2
      for t in BicoloredRootedTreeIterator(order)
        @test residual_order_condition(t, ark) ≈ 0 atol=eps()
      end
    end

    let order=3
      res = 0.0
      for t in BicoloredRootedTreeIterator(order)
        res += abs(residual_order_condition(t, ark))
      end
      @test res > 10*eps()
    end
  end

  @testset "AdditiveRungeKuttaMethod, Lobatto IIIA-IIIB pair (s = 3)" begin
    # Hairer, Lubich, Wanner (2002)
    # Geometric numerical integration
    # Table II.2.2
    As = [
      [0 0 0; 5//24 1//3 -1//24; 1//6 2//3 1//6],
      [1//6 -1//6 0; 1//6 1//3 0; 1//6 5//6 0]
    ]
    bs = [
      [1//6, 2//3, 1//6],
      [1//6, 2//3, 1//6]
    ]
    ark = AdditiveRungeKuttaMethod(As, bs)

    for order in 1:4
      for t in BicoloredRootedTreeIterator(order)
        @test residual_order_condition(t, ark) ≈ 0 atol=eps()
      end
    end

    let order=5
      res = 0.0
      for t in BicoloredRootedTreeIterator(order)
        res += abs(residual_order_condition(t, ark))
      end
      @test res > 10*eps()
    end
  end

  @testset "AdditiveRungeKuttaMethod, Griep3" begin
    # Oswald Knoth and J. Wensch (2014)
    # "Generalized Split-Explicit Runge-Kutta Methods for the
    # Compressible Euler Equations".
    # Monthly Weather Review, 142, 2067-2081
    A_explicit = @SArray [0 0 0; 1//2 0 0; -1 2 0]
    b_explicit = @SArray [1//6, 2//3, 1//6]
    rk_explicit = RungeKuttaMethod(A_explicit, b_explicit)
    β = sqrt(3) / 3
    A_implicit = @SArray [0 0 0; -β/2 (1+β)/2 0; (3+5β)/2 -1-3β (1+β)/2]
    b_implicit = @SArray [1//6, 2//3, 1//6]
    rk_implicit = RungeKuttaMethod(A_implicit, b_implicit)
    ark = AdditiveRungeKuttaMethod([rk_explicit, rk_implicit])

    for order in 1:3
      for t in BicoloredRootedTreeIterator(order)
        @test residual_order_condition(t, ark) ≈ 0 atol=eps()
      end
    end

    let order=4
      res = 0.0
      for t in BicoloredRootedTreeIterator(order)
        res += abs(residual_order_condition(t, ark))
      end
      @test res > 10*eps()
    end
  end

  @testset "AdditiveRungeKuttaMethod, ARK3(2)4L[2]SA" begin
    # Kennedy, Christopher A., and Mark H. Carpenter.
    # "Additive Runge-Kutta schemes for convection-diffusion-reaction equations."
    # Applied Numerical Mathematics 44, no. 1-2 (2003): 139-181.
    # https://doi.org/10.1016/S0168-9274(02)00138-1
    A_explicit = @SArray [
      0 0 0 0
      1767732205903/2027836641118 0 0 0
      5535828885825/10492691773637 788022342437/10882634858940 0 0
      6485989280629/16251701735622 -4246266847089/9704473918619 10755448449292/10357097424841 0
    ]
    b_explicit = @SArray [
      1471266399579/7840856788654,
      -4482444167858/7529755066697,
      11266239266428/11593286722821,
      1767732205903/4055673282236,
    ]
    rk_explicit = RungeKuttaMethod(A_explicit, b_explicit)
    A_implicit = @SArray [
      0 0 0 0
      1767732205903/4055673282236 1767732205903/4055673282236 0 0
      2746238789719/10658868560708 -640167445237/6845629431997 1767732205903/4055673282236 0
      1471266399579/7840856788654 -4482444167858/7529755066697 11266239266428/11593286722821 1767732205903/4055673282236
    ]
    b_implicit = @SArray [
      1471266399579/7840856788654,
      -4482444167858/7529755066697,
      11266239266428/11593286722821,
      1767732205903/4055673282236,
    ]
    rk_implicit = RungeKuttaMethod(A_implicit, b_implicit)
    ark = AdditiveRungeKuttaMethod([rk_explicit, rk_implicit])

    for order in 1:3
      for t in BicoloredRootedTreeIterator(order)
        @test residual_order_condition(t, ark) ≈ 0 atol=eps()
      end
    end

    let order=4
      res = 0.0
      for t in BicoloredRootedTreeIterator(order)
        res += abs(residual_order_condition(t, ark))
      end
      @test res > 10*eps()
    end
  end

  @testset "AdditiveRungeKuttaMethod, SSPRK33 three times" begin
    # Using the same method multiple times is equivalent to using a plain RK
    # method without any splitting/partitioning/decomposition
    A = @SArray [0 0 0; 1 0 0; 1/4 1/4 0]
    b = @SArray [1/6, 1/6, 2/3]
    rk = RungeKuttaMethod(A, b)
    ark = AdditiveRungeKuttaMethod([rk, rk, rk])

    let t = rootedtree([1, 2, 2], [1, 2, 3])
      @test residual_order_condition(t, ark) ≈ 0 atol=eps()
    end

    let t = rootedtree([1, 2, 3, 2], [1, 2, 3, 1])
      @test abs(residual_order_condition(t, ark)) > 10*eps()
    end
  end
end # @testset "Order conditions"

@testset "plots" begin
  @testset "RootedTree" begin
    plot(rootedtree(Int[]))

    for order in 1:4
      for t in RootedTreeIterator(order)
        plot(t)
      end
    end
  end

  @testset "ColoredRootedTree" begin
    let t = rootedtree(Int[], Bool[])
      plot(t)
    end

    let t = rootedtree([1], [1])
      plot(t)
    end

    let t = rootedtree([1], [2])
      plot(t)
    end

    let t = rootedtree([1], [3])
      plot(t)
    end

    let t = rootedtree([1, 2], [1, 1])
      plot(t)
    end

    let t = rootedtree([1, 2], [1, 2])
      plot(t)
    end

    let t = rootedtree([1, 2], [3, 1])
      plot(t)
    end

    let t = rootedtree([1, 2, 2], [2, 1, 1])
      plot(t)
    end

    let t = rootedtree([1, 2, 2], [2, 1, 2])
      plot(t)
    end

    let t = rootedtree([1, 2, 3], [3, 2, 1])
      plot(t)
    end
  end
end # @testset "plots"

end # @testset "RootedTrees"


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of default printing.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Default" begin

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data)
    @test result == expected

    # Without a newline at end.
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘"""

    result = pretty_table(String, data, newline_at_end = false)
    @test result == expected
end

@testset "Pre-defined formats" begin

    # ascii_dots
    # ==========================================================================
    expected = """
.....................................
: Col. 1 : Col. 2 : Col. 3 : Col. 4 :
:........:........:........:........:
:      1 :  false :    1.0 :      1 :
:      2 :   true :    2.0 :      2 :
:      3 :  false :    3.0 :      3 :
:      4 :   true :    4.0 :      4 :
:      5 :  false :    5.0 :      5 :
:      6 :   true :    6.0 :      6 :
:........:........:........:........:
"""
    result = pretty_table(String, data, tf = tf_ascii_dots)
    @test result == expected

    # ascii_rounded
    # ==========================================================================
    expected = """
.--------.--------.--------.--------.
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
:--------+--------+--------+--------:
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
'--------'--------'--------'--------'
"""
    result = pretty_table(String, data, tf = tf_ascii_rounded)
    @test result == expected

    # borderless
    # ==========================================================================
    expected = """
  Col. 1   Col. 2   Col. 3   Col. 4

       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
"""
    result = pretty_table(String, data, tf = tf_borderless)
    @test result == expected

    # compact
    # ==========================================================================
    expected = """
 -------- -------- -------- --------
  Col. 1   Col. 2   Col. 3   Col. 4
 -------- -------- -------- --------
       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
 -------- -------- -------- --------
"""
    result = pretty_table(String, data, tf = tf_compact)
    @test result == expected

    # markdown
    # ==========================================================================
    expected = """
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
|--------|--------|--------|--------|
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
"""
    result = pretty_table(String, data, tf = tf_markdown)
    @test result == expected

    # matrix
    # ==========================================================================

    expected = """
┌                  ┐
│ 1  false  1.0  1 │
│ 2   true  2.0  2 │
│ 3  false  3.0  3 │
│ 4   true  4.0  4 │
│ 5  false  5.0  5 │
│ 6   true  6.0  6 │
└                  ┘
"""
    result = pretty_table(String, data, tf = tf_matrix, noheader = true)
    @test result == expected

    # mysql
    # ==========================================================================
    expected = """
+--------+--------+--------+--------+
| Col. 1 | Col. 2 | Col. 3 | Col. 4 |
+--------+--------+--------+--------+
|      1 |  false |    1.0 |      1 |
|      2 |   true |    2.0 |      2 |
|      3 |  false |    3.0 |      3 |
|      4 |   true |    4.0 |      4 |
|      5 |  false |    5.0 |      5 |
|      6 |   true |    6.0 |      6 |
+--------+--------+--------+--------+
"""
    result = pretty_table(String, data, tf = tf_mysql)
    @test result == expected

    # simple
    # ==========================================================================
    expected = """
========= ======== ======== =========
  Col. 1   Col. 2   Col. 3   Col. 4
========= ======== ======== =========
       1    false      1.0        1
       2     true      2.0        2
       3    false      3.0        3
       4     true      4.0        4
       5    false      5.0        5
       6     true      6.0        6
========= ======== ======== =========
"""
    result = pretty_table(String, data, tf = tf_simple)
    @test result == expected

    # unicode_rounded
    # ==========================================================================
    expected = """
╭────────┬────────┬────────┬────────╮
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
╰────────┴────────┴────────┴────────╯
"""
    result = pretty_table(String, data, tf = tf_unicode_rounded)
    @test result == expected

    # Custom formats
    # ==========================================================================

    expected = """
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
"""

    tf = TextFormat(hlines = [:header])
    result = pretty_table(String, data, tf = tf)
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    tf = TextFormat(hlines = [:begin,:end])
    result = pretty_table(String, data, tf = tf)
    @test result == expected
end

@testset "Dictionaries" begin
    dict = Dict{Int64,String}(
        1 => "Jan",
        2 => "Feb",
        3 => "Mar",
        4 => "Apr",
        5 => "May",
        6 => "Jun"
    )

    expected = """
┌───────┬────────┐
│  Keys │ Values │
│ Int64 │ String │
├───────┼────────┤
│     1 │    Jan │
│     2 │    Feb │
│     3 │    Mar │
│     4 │    Apr │
│     5 │    May │
│     6 │    Jun │
└───────┴────────┘
"""

    result = pretty_table(String, dict, sortkeys = true)
    @test result == expected
end

@testset "Vectors" begin

    vec = 0:1:10

    expected = """
┌────────┐
│ Col. 1 │
├────────┤
│      0 │
│      1 │
│      2 │
│      3 │
│      4 │
│      5 │
│      6 │
│      7 │
│      8 │
│      9 │
│     10 │
└────────┘
"""

    result = pretty_table(String, vec)
    @test result == expected

    expected = """
┌─────┬────────┐
│ Row │ Col. 1 │
├─────┼────────┤
│   1 │   0    │
│   2 │   1    │
│   3 │   2    │
│   4 │   3    │
│   5 │   4    │
│   6 │   5    │
│   7 │   6    │
│   8 │   7    │
│   9 │   8    │
│  10 │   9    │
│  11 │   10   │
└─────┴────────┘
"""

    result = pretty_table(String, vec; alignment = :c, show_row_number = true)
    @test result == expected

    expected = """
┌────────────┐
│     Header │
│ Sub-header │
├────────────┤
│          0 │
│          1 │
│          2 │
│          3 │
│          4 │
│          5 │
│          6 │
│          7 │
│          8 │
│          9 │
│         10 │
└────────────┘
"""

    result = pretty_table(
        String,
        vec;
        header = (["Header"], ["Sub-header"])
    )
    @test result == expected

    @test_throws Exception pretty_table(
        vec;
        header = ["1", "1"]
    )
end

@testset "Print missing, nothing, and #undef" begin

    matrix = Matrix{Any}(undef,3,3)
    matrix[1,1:2] .= missing
    matrix[2,1:2] .= nothing
    matrix[3,1]   = missing
    matrix[3,2]   = nothing

    expected = """
┌─────────┬─────────┬────────┐
│  Col. 1 │  Col. 2 │ Col. 3 │
├─────────┼─────────┼────────┤
│ missing │ missing │ #undef │
│ nothing │ nothing │ #undef │
│ missing │ nothing │ #undef │
└─────────┴─────────┴────────┘
"""

    result = pretty_table(String, matrix)
    @test result == expected
end

@testset "Overwrite" begin
    result = pretty_table(
        String,
        data;
        header = (["A", "B", "C", "D"], ["E", "F", "G", "H"]),
        body_hlines = collect(1:1:6)
    )

    num_lines = length(findall(x->x == '\n', result))

    io = IOBuffer()
    pretty_table(
        io,
        data;
        header = (["A", "B", "C", "D"], ["E", "F", "G", "H"]),
        body_hlines = collect(1:1:6),
        overwrite = true
    )

    io_result = String(take!(io))

    @test io_result == ("\e[1F\e[2K"^(num_lines) * result)
end

ext: jl


text:
module FiveThreeOne

include("TextTable.jl")
import .TextTable
using .TextTable
export TextTable

using Base: Float64
using Printf:@printf

import YAML

export print_main_lift_table, MainLiftSet

struct Order531 end
struct Order351 end
struct Week1 end
struct Week2 end
struct Week3 end

struct Reps
    number::Union{Int, UnitRange{Int}, AbstractString}
    is_pr_set::Bool
end

Reps(number) = Reps(number, false)

function Base.show(io::IO, reps::Reps)
    if reps.number isa UnitRange{Int}
        print(io, first(reps.number), "-", last(reps.number))
    else
        print(io, reps.number)
    end
    reps.is_pr_set && print(io, "+")
end

struct MainLift
    percentage::Float64
    weight::Float64
    sets::Int
    reps::Reps
    e1rm_to_beat::Union{Nothing, Float64}
end

MainLift(percentage, weight, sets, reps::Int) = MainLift(percentage, weight, sets, Reps(reps), nothing)
MainLift(percentage, weight, sets, reps::UnitRange{Int}) = MainLift(percentage, weight, sets, Reps(reps), nothing)
MainLift(percentage, weight, sets, reps::Reps) = MainLift(percentage, weight, sets, reps, nothing)
MainLift(percentage, weight, sets, reps::Int, e1rm_to_beat) = MainLift(percentage, weight, sets, Reps(reps), e1rm_to_beat)
MainLift(percentage, weight, sets, reps::UnitRange{Int}, e1rm_to_beat) = MainLift(percentage, weight, sets, Reps(reps), e1rm_to_beat)

struct AssistanceLift
    name::AbstractString
    weight::Union{Float64, Int, String}
    sets::Int
    reps::Reps
end

AssistanceLift(percentage, weight, sets, reps::Union{Int, UnitRange{Int}, AbstractString}) = AssistanceLift(percentage, weight, sets, Reps(reps, false))

function routine_from_file(; days, main_lifts, maxes_file, assistance_file, order, print_config=nothing)
    function _routine_for_week(week)
        active_print_config = Dict(
            :number_columns => 4,
        )
        if print_config !== nothing
            for k in keys(active_print_config)
                if haskey(print_config, k)
                    active_print_config[k] = print_config[k]
                end
            end
        end
        maxes_data = YAML.load_file(maxes_file)
        assistance_data = YAML.load_file(assistance_file)

        main_sets = Vector{FiveThreeOne.MainLift}[]
        assistance_sets = Vector{FiveThreeOne.AssistanceLift}[]
        for day in days
            day_mains = Vector{MainLift}[]
            push!(day_mains, warmup_sets(maxes_data["training"][day]))
            for lift in main_lifts
                push!(day_mains, lift(maxes_data, day, week, order))
            end
            day_mains = vcat(day_mains...)
            @show day_mains
            push!(main_sets, day_mains)
            day_assistance = AssistanceLift[]
            for lift in assistance_data[day]
                push!(day_assistance, AssistanceLift(lift...))
            end
            push!(assistance_sets, day_assistance)
        end
        print_routine(
            days,
            main_sets,
            assistance_sets,
            n_columns=active_print_config[:number_columns])
    end
    return _routine_for_week
end

function make_routine(training_maxes, supplemental, order, assistance, e1rm_for_pr=nothing; print_config=nothing)
    active_print_config = Dict(
        :number_columns => 4,
    )
    if print_config !== nothing
        for k in keys(active_print_config)
            if haskey(print_config, k)
                active_print_config[k] = print_config[k]
            end
        end
    end
    names = Vector{AbstractString}
    lifts = Vector{MainLift}
    assistance_lifts = Vector{AssistanceLift}

    if e1rm_for_pr === nothing
        e1rm_for_pr = repeat([nothing], length(training_maxes))
        has_pr_sets = false
    else
        has_pr_sets = true
    end

    function _print_routine(week; fives_progression=false)
        names = [n for (n, _) in training_maxes]
        if fives_progression
            _main = ((tm, week) -> main_lifts_5pro(tm, week, order))
        else
            _main = ((tm, week) -> main_lifts(tm, week, order, has_pr_sets, e1rm_to_beat))
        end
        main_sets = [
            vcat(warmup_sets(tm),
                 _main(tm, week),
                 supplemental(tm, week, order))
            for ((_, tm), lift_e1rm) in zip(training_maxes, e1rm_for_pr)
        ]
        assistance_sets = [
            [AssistanceLift(each_lift...) for each_lift in daily_assistance]
            for daily_assistance in assistance
        ]
        print_routine(
            names,
            main_sets,
            assistance_sets,
            n_columns=active_print_config[:number_columns])
    end

    return _print_routine
end

function make_main_lift_cell(name, lifts::Vector{MainLift})
    header_cell = TextCell([name])
    percentages = TextCell(vcat(
        ["%"],
        [string(Int(round(lift.percentage))) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    set_width!(percentages, 4)
    weights = TextCell(vcat(
        ["Weight"],
        [string(lift.weight) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    sets = TextCell(vcat(
        ["Sets"],
        [string(lift.sets) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    reps = TextCell(vcat(
        ["Reps"],
        [string(lift.reps) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    lift_details = hmerge([percentages, weights, sets, reps]; sep="  ")

    # append pr details if there's a pr set
    pr_content = AbstractString[]
    for lift in lifts
        if lift.reps.is_pr_set
            pr_goal = pr_set_goal_reps(lift.weight, lift.e1rm_to_beat)
            push!(pr_content, "$(pr_goal.reps) reps for E1RM $(pr_goal.new_e1rm)")
        end
    end
    if length(pr_content) > 0
        pr_cell = TextCell(pr_content)
        lift_details = vmerge(lift_details, pr_cell)
    end

    # add header to complete
    return vmerge(header_cell, lift_details; sep="-")
end

function make_assistance_cell(lifts::Vector{AssistanceLift})
    header_cell = TextCell(["Assistance"])
    weights = TextCell(vcat(
        ["Weight"],
        [string(lift.weight) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    sets = TextCell(vcat(
        ["Sets"],
        [string(lift.sets) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    reps = TextCell(vcat(
        ["Reps"],
        [string(lift.reps) for lift in lifts]),
        horizontal_alignment=align_right,
        )
    lift_details = hmerge([weights, sets, reps]; horizontal_alignment=align_right, sep="  ")

    names = TextCell([lift.name for lift in lifts])
    assistance = interleave(lift_details, names; horizontal_alignment=align_right)
    return vmerge(header_cell, assistance; sep="-")
end

function print_routine(
    names::Vector{<:AbstractString},
    main_lifts::Vector{Vector{MainLift}},
    assistance::Vector{Vector{AssistanceLift}};
    n_columns::Int=length(names),
    row_sep=" ",
    )
    main_lift_cells = [make_main_lift_cell(name, lifts) for (name, lifts) in zip(names, main_lifts)]
    assistance_cells = [make_assistance_cell(lifts) for lifts in assistance]
    daily_cells = [vmerge(main, assistance; sep="=") for (main, assistance) in zip(main_lift_cells, assistance_cells)]
    n_rows = Int(ceil(length(names) / n_columns))
    row_cells = Vector{TextCell}(undef, n_rows)
    groups = [(1+(i-1)*n_columns):i*n_columns for i in 1:n_rows]
    for i in 1:n_rows
        j1 = first(groups[i])
        j2 = min(last(groups[i]), length(daily_cells))
        row_cells[i] = hmerge(daily_cells[j1:j2], sep=" | ")
    end
    if length(row_cells) > 1
        routine = vmerge(row_cells; sep=row_sep)
    else
        routine = row_cells[1]
    end
    println(format_text(routine))
end

make_single_sets(percentages, weights, reps) = collect(map(MainLift, percentages, weights, repeat([1], length(percentages)), reps))
function make_single_sets(percentages, weights, reps, e1rm_to_beat)
    collect(map(MainLift, percentages, weights, repeat([1], length(percentages)), reps, repeat([e1rm_to_beat], length(percentages))))
end
round_to_half_pound(x) = round(x * 2, RoundNearest) / 2
make_weights(percentages, training_max) = round_to_half_pound.(percentages / 100 .* training_max)

function warmup_sets(training_max)
    percentages = [40, 50, 60]
    reps = [5,5,3]
    weights = make_weights(percentages, training_max)
    return make_single_sets(percentages, weights, reps)
end

function main_lifts(training_max, week, order, pr_sets, e1rm_to_beat)
    make_pr_reps(reps) = Reps(reps, pr_sets)

    if order === Order351 && week === Week1
        week = Week2
    elseif order === Order351 && week === Week2
        week = Week1
    end
    if week === Week1
        percentages = [65, 75, 85]
        reps = [5,5, make_pr_reps(5)]
    elseif week === Week2
        percentages = [70, 80, 90]
        reps = [3,3, make_pr_reps(3)]
    elseif week === Week3
        percentages = [75, 85, 95]
        reps = [5,3, make_pr_reps(1)]
    else
        throw(DomainError("Argument `week` must be Week1, Week2, or Week3."))
    end
    weights = make_weights(percentages, training_max)
    return make_single_sets(percentages, weights, reps, e1rm_to_beat)
end

function main_lifts_5pro(training_max, week, order)
    if order === Order351 && week === Week1
        week = Week2
    elseif order === Order351 && week === Week2
        week = Week1
    end
    if week === Week1
        percentages = [65, 75, 85]
        reps = [5,5, 5]
    elseif week === Week2
        percentages = [70, 80, 90]
        reps = [5,5, 5]
    elseif week === Week3
        percentages = [75, 85, 95]
        reps = [5,5, 5]
    else
        throw(DomainError("Argument `week` must be Week1, Week2, or Week3."))
    end
    weights = make_weights(percentages, training_max)
    return make_single_sets(percentages, weights, reps)
end

function deload_lifts(training_max)
    percentages = [70, 80, 90, 100]
    reps = [5, 3, 1, 1]
    weights = make_weights(percentages, training_max)
    return make_single_sets(percentages, weights, reps)
end

function training_max_test_lifts(training_max)
    percentages = [70, 80, 90, 100]
    reps = [5, 5, 5, 3:5]
    weights = make_weights(percentages, training_max)
    return make_single_sets(percentages, weights, reps)
end

function widowmaker(training_max, week, order)
    if order === Order351 && week === Week1
        week = Week2
    elseif order === Order351 && week === Week2
        week = Week1
    end
    if week === Week1
        percentages = [65]
        reps = [15:20]
    elseif week === Week2
        percentages = [70]
        reps = [15:20]
    elseif week === Week3
        percentages = [75]
        reps = [15:20]
    else
        throw(DomainError("Argument `week` must be Week1, Week2, or Week3."))
    end
    weights = make_weights(percentages, training_max)
    return make_single_sets(percentages, weights, reps)
end

function boringbutbig_light(training_max, week, order)
    if order === Order351 && week === Week1
        week = Week2
    elseif order === Order351 && week === Week2
        week = Week1
    end
    if week === Week1
        percentages = 40
    elseif week === Week2
        percentages = 50
    elseif week === Week3
        percentages = 60
    else
        throw(DomainError("Argument `week` must be Week1, Week2, or Week3."))
    end
    weights = make_weights(percentages, training_max)
    return [MainLift(percentages, weights, 5, 10)]
end

e1rm(weight, reps::Int) = weight * reps * 0.0333 + weight
e1rm(weight, reps::UnitRange{Int}) = collect(map(x -> e1rm(weight, x), reps))
e1rm(weight, reps::Reps) = e1rm(weight, reps.number)

struct PRGoal
    reps::Reps
    new_e1rm::Float64
end

function pr_set_goal_reps(pr_set_weight, e1rm_to_beat)
    rep_range = 1:50
    e1rm_above_goal = e1rm(pr_set_weight, rep_range) .>= e1rm_to_beat
    goal_reps = collect(rep_range)[e1rm_above_goal][1]
    return PRGoal(Reps(goal_reps), e1rm(pr_set_weight, goal_reps))
end

end


text:
s = open("input.txt") do file
    readlines(file);
end

function EvalExpression(expr)
    parts = split(expr, " ");
    curVal = parse(Int, parts[1])

    i = 2;
    while i <= length(parts)
        op = parts[i];
        num = parse(Int, parts[i+1]);

        if op == "+"
            curVal += num;
        else
            curVal *= num;
        end

        i += 2;
    end

    return curVal
end

function DoTheMath(expression)
    expressionStack = [];
    currentExpr = "";
    for i = 1:length(expression)
        if expression[i] == '('
            push!(expressionStack, currentExpr);
            currentExpr = "";
        elseif expression[i] == ')'
            val = EvalExpression(currentExpr);
            currentExpr = pop!(expressionStack) * "$(val)";
        else
            currentExpr *= string(expression[i])
        end
    end

    val = EvalExpression(currentExpr);
    return val;
end

answer = sum(map(DoTheMath, s))


text:
__precompile__(true)
module LightGraphs

using GZip
using Distributions: Binomial
using Base.Collections
using LightXML
using ParserCombinator: Parsers.DOT, Parsers.GML
using StatsBase: fit, Histogram

import Base: write, ==, <, *, ≈, isless, issubset, union, intersect,
            reverse, reverse!, blkdiag, getindex, setindex!, show, print, copy, in,
            sum, size, sparse, eltype, length, ndims, issymmetric, transpose,
            ctranspose, join, start, next, done, eltype, get, A_mul_B!

# core
export SimpleGraph, Edge, Graph, DiGraph, vertices, edges, src, dst,
fadj, badj, in_edges, out_edges, has_vertex, has_edge, is_directed,
nv, ne, add_edge!, rem_edge!, add_vertex!, add_vertices!,
indegree, outdegree, degree, degree_histogram, density, Δ, δ,
Δout, Δin, δout, δin, neighbors, in_neighbors, out_neighbors,
common_neighbors, all_neighbors, has_self_loop, rem_vertex!,

# distance
eccentricity, diameter, periphery, radius, center,

# distance between graphs
spectral_distance, edit_distance,

# edit path cost functions
MinkowskiCost, BoundedMinkowskiCost,

# operators
complement, reverse, reverse!, blkdiag, union, intersect,
difference, symmetric_difference,
induced_subgraph, join, tensor_product, cartesian_product,
crosspath,

# graph visit
SimpleGraphVisitor, TrivialGraphVisitor, LogGraphVisitor,
discover_vertex!, open_vertex!, close_vertex!,
examine_neighbor!, visited_vertices, traverse_graph!, traverse_graph_withlog,

# bfs
BreadthFirst, gdistances, gdistances!, bfs_tree, is_bipartite, bipartite_map,

# dfs
DepthFirst, is_cyclic, topological_sort_by_dfs, dfs_tree,

# random
randomwalk, saw, non_backtracking_randomwalk,

# connectivity
connected_components, strongly_connected_components, weakly_connected_components,
is_connected, is_strongly_connected, is_weakly_connected, period,
condensation, attracting_components, neighborhood, egonet, isgraphical,

# maximum_adjacency_visit
MaximumAdjacency, AbstractMASVisitor, mincut, maximum_adjacency_visit,

# a-star, dijkstra, bellman-ford, floyd-warshall
a_star, dijkstra_shortest_paths,
bellman_ford_shortest_paths, has_negative_edge_cycle, enumerate_paths,
floyd_warshall_shortest_paths,

# centrality
betweenness_centrality, closeness_centrality, degree_centrality,
indegree_centrality, outdegree_centrality, katz_centrality, pagerank,

# spectral
adjacency_matrix,laplacian_matrix, adjacency_spectrum, laplacian_spectrum,
CombinatorialAdjacency, non_backtracking_matrix, incidence_matrix,
nonbacktrack_embedding, Nonbacktracking,
contract,

# astar
a_star,

# persistence
# readgraph, readgraphml, readgml, writegraphml, writegexf, readdot,
load, save,

# flow
maximum_flow, EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm,
multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm,

# randgraphs
erdos_renyi, watts_strogatz, random_regular_graph, random_regular_digraph, random_configuration_model,
StochasticBlockModel, make_edgestream, nearbipartiteSBM, blockcounts, blockfractions,
stochastic_block_model, barabasi_albert, barabasi_albert!, static_fitness_model, static_scale_free,

#community
modularity, core_periphery_deg,
local_clustering,local_clustering_coefficient, global_clustering_coefficient, triangles,
label_propagation, maximal_cliques,

#generators
CompleteGraph, StarGraph, PathGraph, WheelGraph, CycleGraph,
CompleteBipartiteGraph, CompleteDiGraph, StarDiGraph, PathDiGraph, Grid,
WheelDiGraph, CycleDiGraph, BinaryTree, DoubleBinaryTree, RoachGraph, CliqueGraph,

#smallgraphs
smallgraph,

# Euclidean graphs
euclidean_graph

"""An optimized graphs package.

Simple graphs (not multi- or hypergraphs) are represented in a memory- and
time-efficient manner with adjacency lists and edge sets. Both directed and
undirected graphs are supported via separate types, and conversion is available
from directed to undirected.

The project goal is to mirror the functionality of robust network and graph
analysis libraries such as NetworkX while being simpler to use and more
efficient than existing Julian graph libraries such as Graphs.jl. It is an
explicit design decision that any data not required for graph manipulation
(attributes and other information, for example) is expected to be stored
outside of the graph structure itself. Such data lends itself to storage in
more traditional and better-optimized mechanisms.
"""
LightGraphs

include("core.jl")
    include("digraph.jl")
    include("graph.jl")
        include("edgeiter.jl")
        include("traversals/graphvisit.jl")
            include("traversals/bfs.jl")
            include("traversals/dfs.jl")
            include("traversals/maxadjvisit.jl")
            include("traversals/randomwalks.jl")
        include("connectivity.jl")
        include("distance.jl")
        include("edit_distance.jl")
        include("shortestpaths/astar.jl")
        include("shortestpaths/bellman-ford.jl")
        include("shortestpaths/dijkstra.jl")
        include("shortestpaths/floyd-warshall.jl")
        include("spectral.jl")
        include("operators.jl")
        include("persistence/common.jl")
            include("persistence/lg.jl")
            include("persistence/dot.jl")
            include("persistence/gexf.jl")
            include("persistence/gml.jl")
            include("persistence/graphml.jl")
            include("persistence/net.jl")
            include("persistence/jld.jl")
        include("generators/staticgraphs.jl")
            include("generators/randgraphs.jl")
            include("generators/euclideangraphs.jl")
            include("generators/smallgraphs.jl")
        include("centrality/betweenness.jl")
        include("centrality/closeness.jl")
        include("centrality/degree.jl")
        include("centrality/katz.jl")
        include("centrality/pagerank.jl")
        include("community/modularity.jl")
        include("community/label_propagation.jl")
        include("community/core-periphery.jl")
        include("community/clustering.jl")
        include("community/cliques.jl")
        include("flow/maximum_flow.jl")
            include("flow/edmonds_karp.jl")
            include("flow/dinic.jl")
            include("flow/boykov_kolmogorov.jl")
            include("flow/push_relabel.jl")
            include("flow/multiroute_flow.jl")
                include("flow/kishimoto.jl")
                include("flow/ext_multiroute_flow.jl")
        include("utils.jl")

end # module


text:
module FrostHelperSpeedJournalTrigger

using ..Ahorn, Maple

@mapdef Trigger "FrostHelper/SpeedChallengeJournal" JournalTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    challengeNames::String="sc2020_beg_rainbowBerryRush", autoAddSid::Bool=true)

const placements = Ahorn.PlacementDict(
    "Speed Challenge Journal (Frost Helper)" => Ahorn.EntityPlacement(
        JournalTrigger,
        "rectangle",
    ),
)

function Ahorn.nodeLimits(trigger::JournalTrigger)
    return 0, 0
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using Base.Iterators: product
import Base: vec, getindex

export TensorProductBasis, TensorProductBasisState, basis_states

"""
A simple covenience wrapper for a tuple of subsystem dimensions for a tensor product space.
"""
struct TensorProductBasis{N}
    dims::Tuple{Vararg{Int, N}}
end

# We may want to paramaterize by the dims. So, hide the implementation
dims(b::TensorProductBasis) = b.dims

# helper contructor to convert from dimensions given in AbstractVector form
TensorProductBasis(dims::AbstractVector{Int}) = TensorProductBasis(Tuple(dims))

abstract type BasisState end

"""
A basis element of a tensor product space.
"""
struct TensorProductBasisState{N} <: BasisState
    basis::TensorProductBasis{N}
    states::Tuple{Vararg{Int, N}}

    function TensorProductBasisState(b::TensorProductBasis{N}, states::Tuple{Vararg{Int, N}}) where N
        # check that given subsystem basis states are compatible with TensorProductBasis dimensions
        @assert all(states .>= 0) && all(states .< dims(b))
        return new{N}(b, states)
    end
end

# helper contructor to convert from dimensions given in AbstractVector form
TensorProductBasisState(b::TensorProductBasis, states::AbstractVector{Int}) = TensorProductBasisState(b, Tuple(states))

"""
    vec(state::TensorProductBasisState)

Create the complex basis vector correspoding to a given TensorProductBasisState.

## args
* `state`: a TensorProductBasisState to convert to a complex state vector.

## returns
The complex basis state vector.
"""
function vec(state::TensorProductBasisState)
    state_vector = zeros(ComplexF64, prod(dims(state.basis)))
    state_vector[index(state)] = 1.0
    return state_vector
end

"""
    basis_states(b::TensorProductBasis)

Enumerate all basis states in a tensor product space in the canonical order used by the Kronecker
product.

## args
* `b`: a TensorProductBasis.

## returns
A vector of TensorProductBasisState.

## example
`[bs.states for bs in basis_states(TensorProductBasis((2,2)))] == [(0,0), (0,1), (1,0), (1,1)]]`.
"""
function basis_states(b::TensorProductBasis)
    return [TensorProductBasisState(b, reverse(x .- 1))
            for x in vec(collect(product([1:dim for dim in reverse(dims(b))]...)))]
end

"""
    index(bs::TensorProductBasisState)

Convert basis state of a tensor product space to an index into an enumerate of all basis states from
the canonical order used by the Kronecker product.

## args
* `bs`: a TensorProductBasisState.

## returns
The linear index of the state into the canonical tensor product basis states.

## example
Consider a tensor product of two qubits `b = TensorProductBasis((2,2))`. Then states `(0,0), (0,1),
(1,0), (1,1)` are numbered `1, 2, 3, 4` so that `index(TensorProductBasisState(b, (1,0)) == 3`.
"""
function index(bs::TensorProductBasisState)
    return LinearIndices(tuple(reverse(dims(bs.basis))...))[reverse(bs.states.+1)...]
end

"""
    getindex(b::TensorProductBasis, i::Int)

Convert an index into a tensor product space to the basis state. Note that
this function can be called with the notation `b[i]`.

## args
* `b`: the TensorProductBasis into which to index.
* `i`: the index of the basis state in the tensor product space.

## returns
The indexed TensorProductBasisState.

## example
Consider a tensor product of two qubits so that `b = TensorProductBasis((2,2))`. Then the basis
states `(0,0), (0,1), (1,0), (1,1)` are numbered `1, 2, 3, 4` so that
`b[3] == TensorProductBasisState(b, (1,0))`.
"""
function getindex(b::TensorProductBasis, i::Int)
    states = reverse(Tuple(CartesianIndices(tuple(reverse(dims(b))...))[i])).-1
    return TensorProductBasisState(b, states)
end


text:
@testset "Linear" begin
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[1.0 0 0; 0 0 0; 0 0 2.0], [2.89241 0; 0 1.89241], [2.35008 0; 0 2.35008], [3.62629 0; 0 1.62629]],
        [0.0, -2.0, 0.0, 3.0, 1.0]))
    MOIT.lin1ftest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 2.0]],
        [3.0, 1.0]))
    MOIT.lin1vtest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[16.0], [3.0], [13.9279 0.0; 0.0 17.9279], [14.3814 0.0; 0.0 17.3814], [25.3959 0.0; 0.0 9.3959], [15.8211 0.0; 0.0 15.8211]],
        [-7.0, -2.0, 4.0, 0.0, 0.0, -7.0]))
    MOIT.lin2ftest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[0.0], [16.0], [3.0], [75.9147 0.0; 0.0 79.9147]],
        [-7.0, -2.0, 4.0, 70.2306]))
    MOIT.lin2vtest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, [-0.5, 0.5]))
    MOIT.lin3test(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        MOI.INFEASIBLE, MOI.INFEASIBLE_POINT, [-1.0]))
    MOIT.lin4test(bridged, config)
end
@testset "SOC" begin
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[1.0 0.7071 0.7071; 0.7071 1.0 -0.0; 0.7071 -0.0 1.0], [2.2375 0.0; 0.0 1.2375], [2.0662 0.0; 0.0 1.3591], [2.0662 0.0; 0.0 1.3591]],
        [1.4142, -0.7071, 1.0, -0.3536, 1.0, -0.7071, -0.3536]))
    MOIT.soc1ftest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[1.0 0.7071 0.7071; 0.7071 1.0 -0.0; 0.7071 -0.0 1.0], [2.2375 0.0; 0.0 1.2375], [2.0662 0.0; 0.0 1.3591], [2.0662 0.0; 0.0 1.3591]],
        [1.4142, -0.7071, 1.0, -0.3536, 1.0, -0.7071, -0.3536]))
    MOIT.soc1vtest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[0.0], [1.0 -0.7071 0.7071; -0.7071 1.0 0.0; 0.7071 0.0 1.0], [0.4046 0.0; 0.0 1.1117], [1.0895 0.0; 0.0 0.3824], [1.3079 0.0; 0.0 0.3079]],
        [-1, -√2, -√2/2, -1, -√2/4, 1, √2/2, -√2/4]))
    MOIT.soc2ntest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[0.0], [1.0 -0.7071 0.7071; -0.7071 1.0 0.0; 0.7071 0.0 1.0], [0.4046 0.0; 0.0 1.1117], [1.0895 0.0; 0.0 0.3824], [1.3079 0.0; 0.0 0.3079]],
        [1, -√2, -√2/2, -1, -√2/4, 1, √2/2, -√2/4]))
    MOIT.soc2ptest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                MOI.INFEASIBLE,
                                                                MOI.INFEASIBLE_POINT,
                                                                [-1.0, 1.0, -0.5, 1.0, -0.5]))
    MOIT.soc3test(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                [[1.0 0.8944 0.4472; 0.8944 1.0 0.0; 0.4472 0.0 1.0], [1.9332 0.0; 0.0 0.9332], [1.868 0.0; 0.0 0.9736], [1.5795 0.0; 0.0 1.1323], [1.868 0.0; 0.0 0.9736], [1.5795 0.0; 0.0 1.1323]],
                                                                [2.2361, 2.0, 1.0, -1.118, 2.0, -0.8944, 1.0, -0.8944, -0.2236]))
    MOIT.soc4test(bridged, config)
end
@testset "RSOC" begin
        MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                    [[0.5 0.7071 0.7071; 0.7071 2.0 0.0; 0.7071 0.0 2.0], [0.8667 0.0; 0.0 0.1596], [0.8667 0.0; 0.0 0.1596]],
                                                                    [-1.4142, 1.0, -0.1768, 1.0, -0.3536, -0.1768]))
        MOIT.rotatedsoc1ftest(bridged, config)
        MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                    [[0.0], [0.0], [0.5 0.7071 0.7071; 0.7071 2.0 0.0; 0.7071 0.0 2.0], [2.6164 0.0; 0.0 1.9093], [2.6164 0.0; 0.0 1.9093]],
                                                                    [5183.15, 5182.44, -1.4142, 1.0, -0.1768, 1.0, -0.3536, -0.1768]))
        MOIT.rotatedsoc1vtest(bridged, config)
        MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                    MOI.INFEASIBLE,
                                                                    tuple(),
                                                                    [141.088, -47.8864, 47.5533, -46.2201]))
        MOIT.rotatedsoc2test(bridged, config)
        # FIXME u >= 0.0 followed by u <= ub. We need to drop support for
        #       SingleVariable-in-LessThan but we need variable bridge otherwise,
        #       it creates a slack
#            MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
#                                                                        [[0.0], [0.0], [3.0], [0.0], [1.0 0.0; 0.0 0.0], [0.7071 1.0 0.0001; 1.0 1.4142 -0.0; 0.0001 -0.0 1.4142], [0.7071 1.7321; 1.7321 4.2426], [4.5908 0.0; 0.0 2.8588]],
#                                                                        [0.0, 0.0, 0.2887, -0.6124, 0.866, -0.3062, 0.0, -0.0001, -0.0, -1.2247, 1.0, -0.2041]))
#            MOIT.rotatedsoc3test(bridged, config)
end
@testset "GeoMean" begin
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                [[0.0], [0.0], [1.4142 2.0; 2.0 2.8284], [1.0 1.4142; 1.4142 2.0], [1.0 1.4142; 1.4142 2.0], [1.8324 0.0; 0.0 0.8324], [1.7758 0.0; 0.0 0.7758], [1.7772 0.0; 0.0 0.7772], [1.7728 0.0; 0.0 0.7728], [2.5017 0.0; 0.0 0.5017], [2.0744 0.0; 0.0 0.6602], [2.0671 0.0; 0.0 0.6529]],
                                                                [1.0, 0.3333, -0.4714, 0.6667, -0.2357, -0.3333, 0.4714, -0.1667, -0.3333, 0.4714, -0.1667]))
    MOIT.geomean1ftest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                                [[0.0], [0.0], [1.4142 2.0; 2.0 2.8284], [1.0 1.4142; 1.4142 2.0], [1.0 1.4142; 1.4142 2.0], [1.8324 0.0; 0.0 0.8324], [1.7758 0.0; 0.0 0.7758], [1.7772 0.0; 0.0 0.7772], [1.7728 0.0; 0.0 0.7728], [2.5017 0.0; 0.0 0.5017], [2.0744 0.0; 0.0 0.6602], [2.0671 0.0; 0.0 0.6529]],
                                                                [1.0, 0.3333, -0.4714, 0.6667, -0.2357, -0.3333, 0.4714, -0.1667, -0.3333, 0.4714, -0.1667]))
    MOIT.geomean1vtest(bridged, config)
end
@testset "PSD" begin
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[1.0 1.0; 1.0 1.0], [2.2633 0.0; 0.0 1.2633], [2.2669 0.0; 0.0 1.2669], [2.2633 0.0; 0.0 1.2633]],
        [-1.0, 2.0, -1.0, -2.0]))
    MOIT.psdt0ftest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[1.0 1.0; 1.0 1.0]],
        [-2.0]))
    MOIT.psdt0vtest(bridged, config)
    # PSD1: see comments in MOI/src/Test/contconic.jl to see where these constants come from
    δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
    ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
    y2 = 1 - ε*δ
    y1 = 1 - √2*y2
    obj = y1 + y2/2
    k = -2*δ/ε
    x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
    α = √(3-2obj-4x2)/2
    β = k*α
    Xv = [α^2, α*β, β^2, α^2, α*β, α^2]
    xv = [√2*x2, x2, x2]
    cX0 = 1+(√2-1)*y2
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[Xv[1] Xv[2] Xv[1]; Xv[2] Xv[3] Xv[2]; Xv[1] Xv[2] Xv[1]],
        [0.254409 0.179895 0.179895; 0.179895 0.254409 0.0; 0.179895 0.0 0.254409],
        [0.349485 0; 0 0.132235],
        [0.114062 0; 0 0.374032],
        [0.419535 0; 0 0.108444],
        [0.313002 0; 0 0.0957514],
        [0.114062 0; 0 0.374032],
        [0.349485 0; 0 0.132235],
        [0.354357 0; 0 0.0999473],
        [0.332125 0; 0 0.152231],
        [0.332125 0; 0 0.152231]],
        [-cX0, -1.35619, -cX0, 2y2, -1.35619, -cX0, -y2/√2, y2, -y2/√2/2, y2, -y2/√2, -y2/√2/2, y2*√2-1, -y2]))
    MOIT.psdt1ftest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[Xv[1] Xv[2] Xv[1]; Xv[2] Xv[3] Xv[2]; Xv[1] Xv[2] Xv[1]],
         [xv[1] xv[2] xv[2]; xv[2] xv[1] 0.0; xv[2] 0.0 xv[1]],
         [0.486623 0; 0 0.486623-xv[1]],
         [0.4475 0; 0 0.4475-xv[2]],
         [0.4475 0; 0 0.4475-xv[2]]],
        [-y2/√2, y2, -y2/√2/2, y2, -y2/√2, -y2/√2/2, y2*√2-1, -y2]))
    MOIT.psdt1vtest(bridged, config)
    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
        [[0.0], [20/3 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 10/3 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0], [4.0981 -2.1213; -2.1213 1.0981], [7.3275 0.0; 0.0 0.6608], [1.7945 0.0; 0.0 1.7945], [4.1881 0.0; 0.0 0.8547], [1.7945 0.0; 0.0 1.7945], [1.7945 0.0; 0.0 1.7945], [1.7945 0.0; 0.0 1.7945], [2.8504 0.0; 0.0 0.9485]],
        [-0.190192, 0, 0.125977, 0, 0.142644, 0.142644, 0.0127405, -0.211325, -0.816497, -0.788675, 0]))
    MOIT.psdt2test(bridged, config)
end


text:
module ParameterPaths

using ..Platforms, BasisFunctions
import Base: getindex, unsafe_getindex, first
import ..Platforms: unsafe_dictionary, correctparamformat, param_first, param_double,
    param_inbetween, hasparam_inbetween, param_increment, dualdictionary, measure,
    SamplingStyle, SolverStyle, ProblemStyle
using BasisFunctions: AbstractMeasure
# For the moment only size is in parameter space. (no spline degree, extensionparameter T, ...)
export ParameterPath
"""
    abstract type ParameterPath{OUT} end

An indexable type that returns something of type `OUT`.
The index can have a lot of complexity, e.g.:
A dictionary expansion and a function together determine a platform parameter
that leads to a expansion with a smaller error, but more degrees of freedom.
"""
abstract type ParameterPath{OUT} end

correctparamformat(::ParameterPath, param) = false

export outtype
outtype(::ParameterPath{OUT}) where OUT = OUT

function getindex(path::ParameterPath, index)
    if !correctparamformat(path, index)
        return ArgumentError("Index $index not compatible with parameterpath $path")
    end
    unsafe_getindex(path, index)
end

export default_param_path
default_param_path(platform::Platform) =
    LinearParameterPath1D()

first(path::ParameterPath) = unsafe_getindex(path,  param_first(path))
function hasparam_inbetween(path::ParameterPath, param1, param2, stoptolerance::Real)
    p = param_inbetween(path, param1, param2)
    norm(param2 - p) >= stoptolerance && norm(param1 - p) >= stoptolerance
end

export LinearParameterPath1D
"""
    struct LinearParameterPath1D{T<:Real} <: ParameterPath{Int}

"""
struct LinearParameterPath1D{T<:Real} <: ParameterPath{Int}
    constant    ::  T
end

LinearParameterPath1D() = LinearParameterPath1D(1)
param_first(path::LinearParameterPath1D) = 1
param_double(path::LinearParameterPath1D, index::Real) =
    2index
param_inbetween(path::LinearParameterPath1D, index1::Real, index2::Real) =
    (index1+index2)/2

hasparam_inbetween(path::LinearParameterPath1D, index1::Real, index2::Real, stoptolerance::Real) =
    index2 - index1 >= 2stoptolerance

param_increment(path::LinearParameterPath1D, param::Real) =
    param + 1/path.constant

unsafe_getindex(path::LinearParameterPath1D, index::Real) =
    round(Int, index*path.constant)

correctparamformat(::LinearParameterPath1D, ::Real) = true

export CartesianParameterPath
"""
    struct CartesianParameterPath{TNTuple{N,Int}} <: ParameterPath{NTuple{N,Int}}

"""
struct CartesianParameterPath{N,T<:NTuple{N,<:Real}} <: ParameterPath{NTuple{N,Int}}
    constant    ::  T
end

CartesianParameterPath{N}() where N = CartesianParameterPath(CartesianIndex{N}(1).I)

param_first(path::CartesianParameterPath) = 1
param_double(path::CartesianParameterPath, param::Real) =
    2param
param_inbetween(path::CartesianParameterPath{N}, param1::Real, param2::Real) where N =
    max(param1, min((param1+param2)/2, minimum(((((param2).^(1/N)).*path.constant.-1)./path.constant).^N)))
param_double(path::CartesianParameterPath, param::NTuple) =
    (2^(1/N)).*param
param_inbetween(path::CartesianParameterPath{N}, param1::NTuple, param2::NTuple) where N =
    max.(param1, min.((param1 .+ param2)./2, (((((param2).^(1/N)).*path.constant.-1)./path.constant).^N)))
function hasparam_inbetween(path::CartesianParameterPath{N}, param1::NTuple, param2::NTuple, stoptolerance::Real) where {N}
    p = param_inbetween(path, param1, param2)
    p != param1 && p != param2
end
param_increment(path::CartesianParameterPath{N}, param::Real) where N =
    minimum((round.(Int,((param).^(1/N)).*path.constant.+1)./path.constant).^N)
param_increment(path::CartesianParameterPath{N}, param::NTuple{N}) where N =
    ((round.(Int,((param).^(1/N)).*path.constant.+1)./path.constant).^N)

unsafe_getindex(path::CartesianParameterPath{N}, index::Real) where N =
    round.(Int, (index^(1/N)).*path.constant)
unsafe_getindex(path::CartesianParameterPath{N}, index::Vararg{<:Real,N}) where N =
    unsafe_getindex(path, index)
unsafe_getindex(path::CartesianParameterPath{N}, index::CartesianIndex{N}) where N =
    unsafe_getindex(path, index.I)
unsafe_getindex(path::CartesianParameterPath{N}, index::NTuple{N,<:Real}) where N =
    round.(Int, (index).*path.constant)
correctparamformat(::CartesianParameterPath, ::Real) = true
correctparamformat(::CartesianParameterPath{N}, ::NTuple{N,<:Real}) where N = true
correctparamformat(::CartesianParameterPath{N}, ::CartesianIndex{N}) where N = true
correctparamformat(::CartesianParameterPath{N}, ::Vararg{<:Real,N}) where N = true

export IncrementalCartesianParameterPath
"""
    struct IncrementalCartesianParameterPath{TNTuple{N,Int}} <: ParameterPath{NTuple{N,Int}}

"""
struct IncrementalCartesianParameterPath{N} <: ParameterPath{NTuple{N,Int}}
end

param_first(path::IncrementalCartesianParameterPath) = 1
param_double(path::IncrementalCartesianParameterPath, param::Real) =
    2param
param_inbetween(path::IncrementalCartesianParameterPath{N}, param1::Real, param2::Real) where N =
    max(param1, min((param1+param2)/2, param2-1))
param_increment(path::IncrementalCartesianParameterPath{N}, param::Real) where N =
    param + 1
function hasparam_inbetween(path::IncrementalCartesianParameterPath{N}, param1::NTuple, param2::NTuple, stoptolerance::Real) where {N}
    p = param_inbetween(path, param1, param2)
    p != param1 && p != param2
end

function unsafe_getindex(path::IncrementalCartesianParameterPath{N}, index::Real) where N
    q, r = divrem(round(Int, index-1), N)
    ntuple(k-> k<=r ? q + 2 : q +1, Val(N))
end

correctparamformat(::IncrementalCartesianParameterPath, ::Real) = true

export NTupleParameterPath
struct NTupleParameterPath{N,OUT} <: ParameterPath{NTuple{N,OUT}}
    first::OUT
    NTupleParameterPath{N}(out::OUT) where {N,OUT} = new{N,OUT}(out)
end

correctparamformat(::NTupleParameterPath{N,OUT}, param) where {N,OUT} = param isa OUT

unsafe_getindex(path::NTupleParameterPath{N}, param) where N =
    ntuple(k->param, Val(N))

export HierarchyPath
"""
"""
struct HierarchyPath{OUT,PATH1<:ParameterPath,PATH2<:ParameterPath{OUT}} <: ParameterPath{OUT}
    pathfirst :: PATH1
    pathsecond :: PATH2
    function HierarchyPath(pathfirst::ParameterPath, pathsecond::ParameterPath)
        if !correctparamformat(pathsecond,first(pathfirst))
            error("parameters produced by the first path are not compatible with the input of the second path.")
        end
        new{outtype(pathsecond),typeof(pathfirst),typeof(pathsecond)}(pathfirst,pathsecond)
    end
end

correctparamformat(path::HierarchyPath, param) =
    correctparamformat(path.pathfirst, param)

unsafe_getindex(path::HierarchyPath, param) =
    unsafe_getindex(path.pathsecond,path.pathfirst[param])

param_first(path::HierarchyPath) =
    param_first(path.pathfirst)

param_double(path::HierarchyPath, param) =
    param_double(path.pathfirst, param)

param_inbetween(path::HierarchyPath, param1, param2) =
    param_inbetween(path.pathfirst, param1, param2)

hasparam_inbetween(path::HierarchyPath, param1, param2, stoptolerance::Real) =
    hasparam_inbetween(path.pathfirst, param1, param2, stoptolerance)

param_increment(path::HierarchyPath, param1) =
    param_increment(path.pathfirst, param1)

export ProductPath
"""
"""
struct ProductPath{N,PATHS,OUT} <: ParameterPath{OUT}
    paths :: PATHS
    function ProductPath(paths::ParameterPath...)
        new{length(paths),typeof(paths),typeof(map(first, paths))}(paths)
    end
end

correctparamformat(path::ProductPath, param) =
    all(map(correctparamformat,path.paths, param))

unsafe_getindex(path::ProductPath, param) =
    map(unsafe_getindex, path.paths, param)

param_first(path::ProductPath) =
    map(param_first, path.paths)

param_double(path::ProductPath, param) =
    map(param_double, path.paths, param)

param_inbetween(path::ProductPath, param1, param2) =
    map(param_inbetween, path.paths, param1, param2)

hasparam_inbetween(path::ProductPath, param1, param2, stoptolerance::Real) =
    any(map(hasparam_inbetween, path.paths, param1, param2, Ref(stoptolerance))...)

default_param_path(platform::ProductPlatform) where N =
    ProductPath(map(default_param_path, elements(platform))...)

export ParametrizedPlatform
struct ParametrizedPlatform{PLATFORM<:Platform, PATH<:ParameterPath} <: Platform
    platform :: PLATFORM
    path     :: PATH

    function ParametrizedPlatform(platform::Platform, path::ParameterPath)
        if !correctparamformat(platform, first(path))
            throw(ArgumentError("The given path returns no parameters that are compatible with the given platform"))
        end
        new{typeof(platform),typeof(path)}(platform, path)
    end
end

export parametrizedplatform
parametrizedplatform(platform::Platform) =
    parametrizedplatform(platform, default_param_path(platform))

parametrizedplatform(platform::ParametrizedPlatform) =
    platform

parametrizedplatform(platform::Platform, path::ParameterPath) =
    ParametrizedPlatform(platform, path)

correctparamformat(platform::ParametrizedPlatform, param) =
    correctparamformat(platform.path, param)

unsafe_dictionary(platform::ParametrizedPlatform, param) =
    unsafe_dictionary(platform.platform, platform.path[param])

default_param_path(platform::ParametrizedPlatform) = platform.path

param_first(platform::ParametrizedPlatform) =
    param_first(platform.path)

param_double(platform::ParametrizedPlatform, param) =
    param_double(platform.path, param)

param_inbetween(platform::ParametrizedPlatform, param1, param2) =
    param_inbetween(platform.path, param1, param2)

hasparam_inbetween(platform::ParametrizedPlatform, param1, param2, stoptolerance::Real) =
    hasparam_inbetween(platform.path, param1, param2, stoptolerance)

param_increment(platform::ParametrizedPlatform, param1) =
    param_increment(platform.path, param1)

dualdictionary(platform::ParametrizedPlatform, param, measure::AbstractMeasure; options...) =
    dualdictionary(platform.platform, param, measure; options...)

measure(platform::ParametrizedPlatform) =
    measure(platform.platform)

SamplingStyle(platform::ParametrizedPlatform) =
    SamplingStyle(platform.platform)

SolverStyle(platform::ParametrizedPlatform, ss::SamplingStyle) =
    SolverStyle(platform.platform, ss)

ProblemStyle(platform::ParametrizedPlatform) =
    ProblemStyle(platform.platform)

end # module


text:
abstract type DiscreteCountingProcess <: Sampleable{Multivariate, Discrete} end

struct BernouliProcess<:DiscreteCountingProcess
    n::Int64
    p::Float64
    # domain::Tuple

    inter_arrival_distribution::Distribution
    success_distribution::Distribution

    function BernoulliProcess(n::Int64, p::Float64)
        new(n,p, Geometric(p), Bernoulli(p))
    end
end

rand(bp::BernouliProcess) = rand(bp.inter_arrival_distribution)
rand(bp::BernouliProcess, t::Int64) = rand(bp.success_distribution)
params(bp::BernouliProcess) = (bp.n, bp.p)


text:
# this uses catalyst reaction networks -- probably not optimised for performance
# function calculate_steady_state_of_ode(repression_threshold,
#                                        hill_coefficient,
#                                        mRNA_degradation_rate,
#                                        protein_degradation_rate,
#                                        basal_transcription_rate,
#                                        translation_rate)
#    # define Hes5 reaction network
#    hes5_network = @reaction_network begin
#       hillr(protein,basal_transcription_rate,repression_threshold,hill_coefficient), ∅ --> mRNA
#       mRNA_degradation_rate, mRNA → ∅
#       protein_degradation_rate, protein → ∅
#       translation_rate, mRNA → mRNA + protein
#    end repression_threshold hill_coefficient mRNA_degradation_rate protein_degradation_rate basal_transcription_rate translation_rate
#
#    # parameters [α,K,n,δ,γ,β,μ]
#    parameters = (repression_threshold,
#                  hill_coefficient,
#                  mRNA_degradation_rate,
#                  protein_degradation_rate,
#                  basal_transcription_rate,
#                  translation_rate)
#    # initial condition [m,p]
#    u₀ = [40.,50000.]
#    tspan = (0.0,100.0)
#    # create the ODEProblem we want to solve
#    prob = ODEProblem(hes5_network, u₀,tspan,parameters)
#    sol = solve(prob)
#
# end #function

"""
Calculate the steady state of the Hes5 ODE system, for a specific set of parameters.

# Arguments

- `repression_threshold::Float64`

- `hill_coefficient::Float64`

- `mRNA_degradation_rate::Float64`

- `protein_degradation_rate::Float64`

- `basal_transcription_rate::Float64`

- `translation_rate::Float64`

# Returns

- `steady_state_solution::Array{Float64,1}`: A 2-element array, giving the steady state for the mRNA and protein respectively.
"""
function calculate_steady_state_of_ode(repression_threshold,
                                       hill_coefficient,
                                       mRNA_degradation_rate,
                                       protein_degradation_rate,
                                       basal_transcription_rate,
                                       translation_rate)

   hill_function(protein,repression_threshold,hill_coefficient) = 1/(1 + (protein/repression_threshold)^hill_coefficient)

   function hes_ode!(du,u,p,t)
      du[1] = p[5]*hill_function(u[2],p[1],p[2]) - p[3]*u[1]
      du[2] = p[6]*u[1] - p[4]*u[2]
   end #function

   parameters = (repression_threshold,
                 hill_coefficient,
                 mRNA_degradation_rate,
                 protein_degradation_rate,
                 basal_transcription_rate,
                 translation_rate)

   # initial condition [m,p]
   initial_guess = [40.,50000.]
   # create the SteadyStateProblem we want to solve
   prob = SteadyStateProblem(hes_ode!, initial_guess, parameters)
   # solve the problem
   steady_state_solution = solve(prob,SSRootfind())
end #function