text:
module RootedTrees

@doc read(joinpath(dirname(@__DIR__), "README.md"), String) RootedTrees

using LinearAlgebra: dot

using Latexify: Latexify
using RecipesBase: RecipesBase
using Requires: @require

export RootedTree, rootedtree, rootedtree!, RootedTreeIterator,
       ColoredRootedTree, BicoloredRootedTree, BicoloredRootedTreeIterator

export butcher_representation

export α, β, γ, density, σ, symmetry, order, root_color

export residual_order_condition, elementary_weight, derivative_weight

export count_trees

export partition_forest, PartitionForestIterator,
       partition_skeleton,
       all_partitions, PartitionIterator

export all_splittings, SplittingIterator

export RungeKuttaMethod, AdditiveRungeKuttaMethod

abstract type AbstractRootedTree end

"""
    RootedTree(level_sequence, is_canonical::Bool=false)

Represents a rooted tree using its level sequence.

# References

- Terry Beyer and Sandra Mitchell Hedetniemi.
  "Constant time generation of rooted trees".
  SIAM Journal on Computing 9.4 (1980): 706-712.
  [DOI: 10.1137/0209055](https://doi.org/10.1137/0209055)
"""
struct RootedTree{T<:Integer, V<:AbstractVector{T}} <: AbstractRootedTree
  level_sequence::V
  iscanonical::Bool

  function RootedTree(level_sequence::V, iscanonical::Bool=false) where {T<:Integer, V<:AbstractVector{T}}
    new{T, V}(level_sequence, iscanonical)
  end
end

"""
    rootedtree(level_sequence)

Construct a canonical [`RootedTree`](@ref) object from a `level_sequence`, i.e.,
a vector of integers representing the levels of each node of the tree.

# References

- Terry Beyer and Sandra Mitchell Hedetniemi.
  "Constant time generation of rooted trees".
  SIAM Journal on Computing 9.4 (1980): 706-712.
  [DOI: 10.1137/0209055](https://doi.org/10.1137/0209055)
"""
rootedtree(level_sequence::AbstractVector) = canonical_representation(RootedTree(level_sequence))

"""
    rootedtree!(level_sequence)

Construct a canonical `RootedTree` object from a `level_sequence` which may be
modified in this process. See also [`rootedtree`](@ref).

# References

- Terry Beyer and Sandra Mitchell Hedetniemi.
  "Constant time generation of rooted trees".
  SIAM Journal on Computing 9.4 (1980): 706-712.
  [DOI: 10.1137/0209055](https://doi.org/10.1137/0209055)
"""
rootedtree!(level_sequence::AbstractVector) = canonical_representation!(RootedTree(level_sequence))

iscanonical(t::RootedTree) = t.iscanonical
#TODO: Validate rooted tree in constructor?

Base.copy(t::RootedTree) = RootedTree(copy(t.level_sequence), t.iscanonical)
Base.similar(t::RootedTree) = RootedTree(similar(t.level_sequence), true)
Base.isempty(t::RootedTree) = isempty(t.level_sequence)
Base.empty(t::RootedTree) = RootedTree(empty(t.level_sequence), iscanonical(t))

@inline function Base.copy!(t_dst::RootedTree, t_src::RootedTree)
  copy!(t_dst.level_sequence, t_src.level_sequence)
  return t_dst
end

"""
    unsafe_deleteat!(t::AbstractRootedTree, i)

Delete the node `i` from the rooted tree `t`. This is an unsafe operation
since the rooted tree will not necessarily be in canonical representation
afterwards, even if the corresponding flag of `t` is set. Use with caution!

!!! warning "Internal interface"
    This function is considered to be an internal implementation detail and
    will not necessarily be stable.
"""
@inline function unsafe_deleteat!(t::RootedTree, i)
  deleteat!(t.level_sequence, i)
  return t
end

"""
    unsafe_resize!(t::AbstractRootedTree, n::Integer)

Resize the rooted tree `t` to `n` nodes. This is an unsafe operation
since the rooted tree will not necessarily be in canonical representation
afterwards, even if the corresponding flag of `t` is set. Use with caution!

!!! warn "Internal interface"
    This function is considered to be an internal implementation detail and
    will not necessarily be stable.
"""
@inline function unsafe_resize!(t::RootedTree, n::Integer)
  resize!(t.level_sequence, n)
  return t
end

"""
    unsafe_copyto!(t_dst::AbstractRootedTree, dst_offset,
                   t_src::AbstractRootedTree, src_offset, N)

Copy `N` nodes from `t_src` starting at offset `src_offset` to `t_dst`
starting at offset `dst_offset`. The types of the rooted trees must match.
For example, you cannot copy a [`ColoredRootedTree`](@ref) to a
[`RootedTree`](@ref).

This is an unsafe operation since the rooted tree `t_dst` will not necessarily
be in canonical representation afterwards, even if the corresponding flag
of `t_dst` is set. Use with caution!

!!! warning "Internal interface"
    This function is considered to be an internal implementation detail and
    will not necessarily be stable.
"""
@inline function unsafe_copyto!(t_dst::RootedTree, dst_offset,
                                t_src::RootedTree, src_offset, N)
  copyto!(t_dst.level_sequence, dst_offset, t_src.level_sequence, src_offset, N)
  return t_dst
end

#  #function RootedTree(sequence::Vector{T}, valid::Bool)
#  function RootedTree(sequence::Array{T,1})
#    length(sequence) < 1 && throw(ArgumentError("Rooted trees must have a root, in particular at least one element!"))#
#
#    ## If there is only one element, the sequence must be valid.
#    #if !valid && length(sequence) > 1
#    #  # Test, whether there is exactly one root element at the beginning of sequence, if necessary.
#    #  root = sequence[1]
#    #  for level in sequence[2:end]
#    #    level <= root && throw(ArgumentError("Rooted trees must have exactly one element at root level at the beginning."))
#    #  end
#    #end
#    # If there is only one element, the sequence must be valid.
#    if length(sequence) > 1
#      # Test, whether there is exactly one root element at the beginning of sequence, if necessary.
#      root = sequence[1]
#      for level in sequence[2:end]
#        level <= root && throw(ArgumentError("Rooted trees must have exactly one element at root level at the beginning."))
#      end
#    end
#
#    new(sequence)
#  end
#RootedTree{T<:Integer}(sequence::Vector{T}) = RootedTree{T}(sequence, false)

function Base.show(io::IO, t::RootedTree{T}) where {T}
  print(io, "RootedTree{", T, "}: ")
  show(io, t.level_sequence)
end

# comparison

"""
    isless(t1::RootedTree, t2::RootedTree)

Compares two rooted trees using a lexicographical comparison of their level
sequences while considering equivalence classes given by different root indices.
"""
function Base.isless(t1::RootedTree, t2::RootedTree)
  if isempty(t1.level_sequence)
    if isempty(t2.level_sequence)
      # empty trees are equal
      return false
    else
      # the empty tree `isless` than any other tree
      return true
    end
  elseif isempty(t2.level_sequence)
    # the empty tree `isless` than any other tree
    return false
  end

  root1_minus_root2 = first(t1.level_sequence) - first(t2.level_sequence)
  for (e1, e2) in zip(t1.level_sequence, t2.level_sequence)
    v1 = e1
    v2 = e2 + root1_minus_root2
    v1 == v2 || return isless(v1, v2)
  end
  return isless(length(t1.level_sequence), length(t2.level_sequence))
end

"""
    ==(t1::RootedTree, t2::RootedTree)

Compares two rooted trees based on their level sequences while considering
equivalence classes given by different root indices.

# Examples

```jldoctest
julia> t1 = rootedtree([1, 2, 3]);

julia> t2 = rootedtree([2, 3, 4]);

julia> t3 = rootedtree([1, 2, 2]);

julia> t1 == t2
true

julia> t1 == t3
false
```
"""
function Base.:(==)(t1::RootedTree, t2::RootedTree)
  length(t1.level_sequence) == length(t2.level_sequence) || return false

  if isempty(t1.level_sequence)
    # empty trees are equal
    return true
  end

  root1_minus_root2 = first(t1.level_sequence) - first(t2.level_sequence)
  for (e1, e2) in zip(t1.level_sequence, t2.level_sequence)
    e1 == e2 + root1_minus_root2 || return false
  end

  return true
end

# Factor out equivalence classes given by different roots
function Base.hash(t::RootedTree, h::UInt)
  # Use a fast path if possible
  if UInt == UInt64 && length(t.level_sequence) <= 16
    return simple_hash(t, h)
  end

  isempty(t.level_sequence) && return h
  root = first(t.level_sequence)
  for l in t.level_sequence
    h = hash(l - root, h)
  end
  return h
end

# Map the level sequence to an unsigned integer by concatenating the bit
# representations of level sequence differences. If the level sequence increases
# from one vertex to the next, it can increase at most by unity. Since we want
# to use simple bits representations, we measure the decrease compared to the
# maximal possible increase.
# The maximal drop in the level sequence is
#   maximal_drop = length(t.level_sequence) - 3
# We need at most
#   number_of_bits = trunc(Int, log2(maximal_drop)) + 1
# bits to represent this. Thus, 64 bit allow us to compute unique hashes for
# level sequence up to length 16 in the following simple way; 64 bit result
# in `number_of_bits = 4` for `maximal_drop = 16 - 3 = 13`.
# For 32 bits, we could use a maximal length of 10 with `number_of_bits = 3`.
# However, most user systems should use 64 bit by default, so we only implement
# this option for simplicity.
@inline function simple_hash(t::RootedTree, h_base::UInt64)
  isempty(t.level_sequence) && return h_base
  h = zero(h_base)
  l_prev = first(t.level_sequence)
  for l in t.level_sequence
    h = (h << 4) | (l_prev + 1 - l)
    l_prev = l
  end
  return hash(h, h_base)
end

# generation and canonical representation
"""
    canonical_representation(t::AbstractRootedTree)

Returns a new tree using the canonical representation of the rooted tree `t`,
i.e., the one with lexicographically biggest level sequence.

See also [`canonical_representation!`](@ref).
"""
function canonical_representation(t::AbstractRootedTree)
  canonical_representation!(copy(t))
end

# A very simple implementation of `canonical_representation!` could read as
# follows.
#   function canonical_representation!(t::RootedTree)
#     subtr = subtrees(t)
#     for i in eachindex(subtr)
#       canonical_representation!(subtr[i])
#     end
#     sort!(subtr, rev=true)
#
#     i = 2
#     for τ in subtr
#       t.level_sequence[i:i+order(τ)-1] = τ.level_sequence
#       i += order(τ)
#     end
#
#     RootedTree(t.level_sequence, true)
#   end
# However, this would create a lot of intermediate allocations, which make it
# rather slow. Since most trees in use are relatively small, we can use a
# non-allocating sorting algorithm instead - although bubble sort is slower in
# general when comparing the complexity with quicksort etc., it will be faster
# here since we can avoid allocations.
"""
    canonical_representation!(t::AbstractRootedTree)

Change the representation of the rooted tree `t` to the canonical one, i.e., the
one with lexicographically biggest level sequence.

See also [`canonical_representation`](@ref).
"""
function canonical_representation!(t::RootedTree, buffer=similar(t.level_sequence))
  # Since we use a recursive implementation, it is useful to exit early for
  # small trees. If there are at most 3 vertices in a valid rooted tree, its
  # level sequence must already be in canonical representation.
  if order(t) <= 3
    return RootedTree(t.level_sequence, true)
  end

  # First, sort all subtrees recursively. Here, we use `view`s to avoid memory
  # allocations.
  # TODO: Assume 1-based indexing in the following
  subtree_root_index = 2
  number_of_subtrees = 0

  while subtree_root_index <= order(t)
    subtree_last_index = _subtree_last_index(subtree_root_index, t.level_sequence)

    # We found a complete subtree
    subtree = RootedTree(view(t.level_sequence, subtree_root_index:subtree_last_index))
    canonical_representation!(subtree, view(buffer, subtree_root_index:subtree_last_index))

    subtree_root_index = subtree_last_index + 1
    number_of_subtrees += 1
  end

  # Next, we need to sort the subtrees of `t` (in lexicographically decreasing
  # order of the level sequences).
  if number_of_subtrees > 1
    # Simple bubble sort that can act in-place, avoiding allocations
    # We keep track of the last index of the last subtree that we need to sort
    # since we know that the last `n` subtrees are already sorted after `n`
    # iterations.
    subtree_last_index_to_sort = order(t)
    swapped = true
    while swapped
      swapped = false

      # Search the first complete subtree
      subtree1_root_index = 2
      subtree1_last_index = 0
      subtree2_last_index = 0
      while subtree1_root_index <= subtree_last_index_to_sort
        subtree1_last_index = _subtree_last_index(subtree1_root_index, t.level_sequence)
        subtree2_last_index = subtree1_last_index

        # Search the next complete subtree
        subtree1_last_index == subtree_last_index_to_sort && break

        subtree2_root_index = subtree1_last_index + 1
        subtree2_last_index = _subtree_last_index(subtree2_root_index, t.level_sequence)

        # Swap the subtrees if they are not sorted correctly
        subtree1 = RootedTree(view(t.level_sequence, subtree1_root_index:subtree1_last_index))
        subtree2 = RootedTree(view(t.level_sequence, subtree2_root_index:subtree2_last_index))
        if isless(subtree1, subtree2)
          copyto!(buffer, 1, t.level_sequence, subtree1_root_index, order(subtree1) + order(subtree2))
          copyto!(t.level_sequence, subtree1_root_index, buffer, order(subtree1) + 1, order(subtree2))
          copyto!(t.level_sequence, subtree1_root_index + order(subtree2), buffer, 1, order(subtree1))
          # `subtree1_root_index` will be updated below using `subtree1_last_index`.
          # Thus, we need to adapt this variable here.
          subtree1_last_index = subtree1_root_index + order(subtree2) - 1
          swapped = true
        end

        # Move on to the next pair of subtrees
        subtree2_last_index == subtree_last_index_to_sort && break
        subtree1_root_index = subtree1_last_index + 1
      end

      # Update the last subtree we need to look at
      subtree_last_index_to_sort = min(subtree1_last_index, subtree2_last_index)
    end
  end

  RootedTree(t.level_sequence, true)
end

@inline function _subtree_last_index(subtree_root_index, level_sequence)
  # The subtree goes up to the next node that has the same (or lower)
  # rank as its root.
  subtree_last_index = subtree_root_index
  subtree_root_level = level_sequence[subtree_root_index]
  while subtree_last_index < length(level_sequence)
    if level_sequence[subtree_last_index + 1] > subtree_root_level
      subtree_last_index += 1
    else
      break
    end
  end
  return subtree_last_index
end

# Allocate global buffer for `canonical_representation!` for each thread
const BUFFER_LENGTH = 64
const CANONICAL_REPRESENTATION_BUFFER = Vector{Vector{Int}}()

function canonical_representation!(t::RootedTree{Int, Vector{Int}})
  if order(t) <= BUFFER_LENGTH
    buffer = CANONICAL_REPRESENTATION_BUFFER[Threads.threadid()]
  else
    buffer = similar(t.level_sequence)
  end
  canonical_representation!(t, buffer)
end

"""
    check_canonical(t::AbstractRootedTree)

Check whether `t` is in canonical representation.

!!! warning "Internal interface"
    This function is considered to be an internal implementation detail and
    will not necessarily be stable.
"""
function check_canonical(t::AbstractRootedTree)
  for subtree in SubtreeIterator(t)
    if !check_canonical(subtree)
      return false
    end
  end

  return issorted(SubtreeIterator(t), rev=true)
end

"""
    normalize_root!(t::AbstractRootedTree, root=one(eltype(t.level_sequence)))

Normalize the level sequence of the rooted tree `t` such that the root is
set to `root`.
"""
function normalize_root!(t::AbstractRootedTree, root=one(eltype(t.level_sequence)))
  t.level_sequence .+= root - first(t.level_sequence)
  t
end

"""
    RootedTreeIterator(order::Integer)

Iterator over all rooted trees of given `order`. The returned trees are views to
an internal tree modified during the iteration. If the returned trees shall be
stored or modified during the iteration, a `copy` has to be made.
"""
struct RootedTreeIterator{T<:Integer}
  order::T
  t::RootedTree{T,Vector{T}}

  function RootedTreeIterator(order::T) where {T<:Integer}
    new{T}(order, RootedTree(Vector{T}(one(T):order), true))
  end
end

Base.IteratorSize(::Type{<:RootedTreeIterator}) = Base.SizeUnknown()
Base.eltype(::Type{RootedTreeIterator{T}}) where T = RootedTree{T,Vector{T}}

@inline function Base.iterate(iter::RootedTreeIterator{T}) where {T}
  iter.t.level_sequence[:] = one(T):iter.order
  (iter.t, false)
end

@inline function Base.iterate(iter::RootedTreeIterator{T}, state) where {T}
  state && return nothing

  two = iter.t.level_sequence[1] + one(T)
  p = 1
  q = 1
  @inbounds for i in 2:length(iter.t.level_sequence)
    if iter.t.level_sequence[i] > two
      p = i
    end
  end
  p == 1 && return nothing

  level_q = iter.t.level_sequence[p] - one(T)
  @inbounds for i in 1:p
    if iter.t.level_sequence[i] == level_q
      q = i
    end
  end

  @inbounds for i in p:length(iter.t.level_sequence)
    iter.t.level_sequence[i] = iter.t.level_sequence[i - (p-q)]
  end

  (iter.t, false)
end

"""
    count_trees(order)

Counts all rooted trees of `order`.
"""
function count_trees(order)
  order < 1 && throw(ArgumentError("The `order` must be at least one."))

  num = 0
  for _ in RootedTreeIterator(order)
    num += 1
  end
  num
end

# subtrees
struct SubtreeIterator{Tree<:AbstractRootedTree}
  t::Tree
end

# optional: define some interface methods such as
# Base.IteratorSize(::Type{<:SubtreeIterator}) = Base.SizeUnknown()
# Base.eltype(::Type{SubtreeIterator})

@inline function Base.iterate(subtrees::SubtreeIterator{<:RootedTree})
  subtree_root_index = firstindex(subtrees.t.level_sequence) + 1
  iterate(subtrees, subtree_root_index)
end

@inline function Base.iterate(subtrees::SubtreeIterator{<:RootedTree}, subtree_root_index)
  level_sequence = subtrees.t.level_sequence

  # terminate the iteration if there are no further subtrees
  if subtree_root_index > lastindex(level_sequence)
    return nothing
  end

  # find the next complete subtree
  subtree_last_index = _subtree_last_index(subtree_root_index, level_sequence)
  subtree = RootedTree(
    view(level_sequence, subtree_root_index:subtree_last_index),
    # if t is in canonical representation, its subtrees are, too
    iscanonical(subtrees.t))

  return (subtree, subtree_last_index + 1)
end

"""
    subtrees(t::RootedTree)

Returns a vector of all subtrees of `t`.
"""
function subtrees(t::RootedTree)
  subtr = typeof(t)[]

  if length(t.level_sequence) < 2
    return subtr
  end

  start = 2
  i = 3
  while i <= length(t.level_sequence)
    if t.level_sequence[i] <= t.level_sequence[start]
      push!(subtr, RootedTree(t.level_sequence[start:i-1]))
      start = i
    end
    i += 1
  end
  push!(subtr, RootedTree(t.level_sequence[start:end]))
end

# partitions
"""
    partition_forest(t::RootedTree, edge_set)

Form the partition forest of the rooted tree `t` where edges marked with `false`
in the `edge_set` are removed. The ith value in the Boolean iterable `edge_set`
corresponds to the edge connecting node `i+1` in the level sequence to its parent.

See also [`partition_skeleton`](@ref), [`PartitionIterator`](@ref), and
[`PartitionForestIterator`](@ref).

# References

Section 2.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)
"""
function partition_forest(t::RootedTree, edge_set)
  @boundscheck begin
    @assert length(t.level_sequence) == length(edge_set) + 1
  end

  level_sequence = copy(t.level_sequence)
  edge_set_copy = copy(edge_set)
  forest = Vector{RootedTree{eltype(level_sequence), typeof(level_sequence)}}()
  partition_forest!(forest, level_sequence, edge_set_copy)
  return forest
end

# Internal implementation that `push!`es to `forest` and modifies the vectors
# `level_sequence` and `edge_set`.
function partition_forest!(forest, level_sequence, edge_set)
  # Iterate over all edges that shall be removed.
  edge_to_remove = findlast(==(false), edge_set)
  while edge_to_remove !== nothing
    # Remember the convention node = edge + 1
    subtree_root_index = edge_to_remove + 1
    subtree_last_index = _subtree_last_index(subtree_root_index, level_sequence)
    subtree = rootedtree!(level_sequence[subtree_root_index:subtree_last_index])

    # Since we search from the end, there is no additional edge that needs to
    # be removed in the current subtree. Thus, we can `push!` it to the `forest`
    # and remove it from the active `level_sequence` and `edge_set`.
    push!(forest, subtree)
    deleteat!(level_sequence, subtree_root_index:subtree_last_index)
    deleteat!(edge_set, subtree_root_index-1:subtree_last_index-1)

    edge_to_remove = findlast(==(false), edge_set)
  end
  push!(forest, rootedtree(level_sequence))
end

"""
    PartitionForestIterator(t::AbstractRootedTree, edge_set)

Lazy iterator representation of the [`partition_forest`](@ref) of the rooted
tree `t`.
Similar to [`RootedTreeIterator`](@ref), you should `copy` the iterates
if you want to store or modify them during the iteration since they may be
views to internal caches.

See also [`partition_forest`](@ref), [`partition_skeleton`](@ref), and
[`PartitionIterator`](@ref).

# References

Section 2.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)
"""
struct PartitionForestIterator{Tree<:AbstractRootedTree}
  t_iter::Tree # return value from `iterate`
  t_temp::Tree # internal temporary buffer
  edge_set::Vector{Bool}
end

function PartitionForestIterator(t::AbstractRootedTree, edge_set)
  t_iter = copy(t)
  t_temp = copy(t)
  PartitionForestIterator(t_iter, t_temp, copy(edge_set))
end

Base.IteratorSize(::Type{<:PartitionForestIterator}) = Base.HasLength()
Base.length(forest::PartitionForestIterator) = count(==(false), forest.edge_set) + 1
Base.eltype(::Type{PartitionForestIterator{Tree}}) where {Tree} = Tree

@inline function Base.iterate(forest::PartitionForestIterator)
  iterate(forest, lastindex(forest.edge_set))
end

@inline function Base.iterate(forest::PartitionForestIterator, search_start)
  t_iter = forest.t_iter
  t_temp = forest.t_temp
  edge_set = forest.edge_set

  # We use `search_start = typemin(Int)` to indicate that we have already
  # returned the final tree in the previous call.
  if search_start == typemin(Int)
    return nothing
  end

  edge_to_remove = findprev(==(false), edge_set, search_start)

  # There are no further edges to remove and we can return the final tree.
  if edge_to_remove === nothing
    unsafe_resize!(t_iter, order(t_temp))
    copy!(t_iter, t_temp)
    canonical_representation!(t_iter)
    return (t_iter, typemin(Int))
  end

  # On to the next subtree
  # Remember the convention node = edge + 1
  subtree_root_index = edge_to_remove + 1
  subtree_last_index = _subtree_last_index(subtree_root_index, t_temp.level_sequence)
  subtree_length = subtree_last_index - subtree_root_index + 1

  # Since we search from the end, there is no additional edge that needs to
  # be removed in the current subtree. Thus, we can return it as the next
  # iterate of the partition forest
  unsafe_resize!(t_iter, subtree_length)
  unsafe_copyto!(t_iter, 1, t_temp, subtree_root_index, subtree_length)
  canonical_representation!(t_iter)

  # Now, we can remove the next subtree iterate from the active
  # level sequence in `t_temp` and the `edge_set`.
  unsafe_deleteat!(t_temp, subtree_root_index:subtree_last_index)
  deleteat!(edge_set, subtree_root_index-1:subtree_last_index-1)

  return (t_iter, edge_to_remove - 1)
end

# necessary for simple and convenient use since the iterates may be modified
function Base.collect(forest::PartitionForestIterator)
  iterates = Vector{eltype(forest)}()
  sizehint!(iterates, length(forest))
  for t in forest
    push!(iterates, copy(t))
  end
  return iterates
end

"""
    partition_skeleton(t::AbstractRootedTree, edge_set)

Form the partition skeleton of the rooted tree `t`, i.e., the rooted tree
obtained by contracting each tree of the partition forest to a single vertex
and re-establishing the edges removed to obtain the partition forest.

See also [`partition_forest`](@ref) and [`PartitionIterator`](@ref).

# References

Section 2.3 (and Section 6.1 for colored trees) 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)
"""
function partition_skeleton(t::AbstractRootedTree, edge_set)
  @boundscheck begin
    @assert order(t) == length(edge_set) + 1
  end

  edge_set_copy = copy(edge_set)
  skeleton = copy(t)
  return partition_skeleton!(skeleton, edge_set_copy)
end

# internal in-place version of partition_skeleton modifying the inputs
function partition_skeleton!(skeleton::AbstractRootedTree, edge_set)
  level_sequence = skeleton.level_sequence

  # Iterate over all edges that shall be kept/contracted.
  # We start the iteration at the end since this will result in less memory
  # moves because we have already reduced the size of the vectors when reaching
  # the beginning.
  edge_to_contract = findlast(edge_set)
  while edge_to_contract !== nothing
    # Contract the corresponding edge by removing the subtree root and promoting
    # the rest of the subtree.
    # Remember the convention node = edge + 1
    subtree_root_index = edge_to_contract + 1
    subtree_last_index = subtree_root_index + 1
    while subtree_last_index <= length(level_sequence)
      if level_sequence[subtree_last_index] > level_sequence[subtree_root_index]
        level_sequence[subtree_last_index] -= 1
        subtree_last_index += 1
      else
        break
      end
    end

    # Remove the root node
    unsafe_deleteat!(skeleton, subtree_root_index)
    deleteat!(edge_set, edge_to_contract)

    edge_to_contract = findprev(edge_set, edge_to_contract - 1)
  end

  # The level sequence `level_sequence` will not automatically be a canonical
  # representation.
  canonical_representation!(skeleton)
end

"""
    all_partitions(t::RootedTree)

Create all partition forests and skeletons of a rooted tree `t`. This returns
vectors of the return values of [`partition_forest`](@ref) and
[`partition_skeleton`](@ref) when looping over all possible edge sets.

See also [`PartitionIterator`](@ref).

# References

Section 2.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)
"""
function all_partitions(t::RootedTree)
  edge_set = zeros(Bool, order(t) - 1)
  forests   = [partition_forest(t, edge_set)]
  skeletons = [partition_skeleton(t, edge_set)]

  for edge_set_value in 1:(2^length(edge_set) - 1)
    binary_digits!(edge_set, edge_set_value)
    push!(forests,   partition_forest(t, edge_set))
    push!(skeletons, partition_skeleton(t, edge_set))
  end

  return (; forests, skeletons)
end

# A helper function to comute the binary representation of an integer `n` as
# a vector of `Bool`s. This is a more efficient version of
#   binary_digits!(digits, n) = digits!(digits, n, base=2)
function binary_digits!(digits::Vector{Bool}, n::Int)
  bit = 1
  for i in eachindex(digits)
    digits[i] = n & bit > 0
    bit = bit << 1
  end
  digits
end

"""
    PartitionIterator(t::AbstractRootedTree)

Iterator over all partition forests and skeletons of the rooted tree `t`.
This is basically a pure iterator version of [`all_partitions`](@ref).
In particular, the partition forest may only be realized as an iterator.
Similar to [`RootedTreeIterator`](@ref), you should `copy` the iterates
if you want to store or modify them during the iteration since they may be
views to internal caches.

See also [`partition_forest`](@ref), [`partition_skeleton`](@ref),
and [`PartitionForestIterator`](@ref).

# References

Section 2.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)
"""
struct PartitionIterator{TreeInput<:AbstractRootedTree, TreeOutput<:AbstractRootedTree}
  t::TreeInput
  forest::PartitionForestIterator{TreeOutput}
  skeleton::TreeOutput
  edge_set::Vector{Bool}
  edge_set_tmp::Vector{Bool}
end

function PartitionIterator(t::AbstractRootedTree)
  skeleton = similar(t)
  edge_set = Vector{Bool}(undef, order(t) - 1)
  edge_set_tmp = similar(edge_set)

  t_forest = similar(t)
  t_temp_forest = similar(t)
  forest = PartitionForestIterator(t_forest, t_temp_forest, edge_set_tmp)
  PartitionIterator{typeof(t), typeof(skeleton)}(t, forest, skeleton, edge_set, edge_set_tmp)
end

# Allocate global buffer for `PartitionIterator` for each thread
const PARTITION_ITERATOR_BUFFER_FOREST_T = Vector{Vector{Int}}()
const PARTITION_ITERATOR_BUFFER_FOREST_T_COLORS = Vector{Vector{Bool}}()
const PARTITION_ITERATOR_BUFFER_FOREST_LEVEL_SEQUENCE = Vector{Vector{Int}}()
const PARTITION_ITERATOR_BUFFER_FOREST_COLOR_SEQUENCE = Vector{Vector{Bool}}()
const PARTITION_ITERATOR_BUFFER_SKELETON = Vector{Vector{Int}}()
const PARTITION_ITERATOR_BUFFER_SKELETON_COLORS = Vector{Vector{Bool}}()
const PARTITION_ITERATOR_BUFFER_EDGE_SET = Vector{Vector{Bool}}()
const PARTITION_ITERATOR_BUFFER_EDGE_SET_TMP = Vector{Vector{Bool}}()

function PartitionIterator(t::RootedTree{Int, Vector{Int}})
  order_t = order(t)

  if order_t <= BUFFER_LENGTH
    id = Threads.threadid()

    buffer_forest_t = PARTITION_ITERATOR_BUFFER_FOREST_T[id]
    resize!(buffer_forest_t, order_t)
    level_sequence  = PARTITION_ITERATOR_BUFFER_FOREST_LEVEL_SEQUENCE[id]
    resize!(level_sequence, order_t)
    buffer_skeleton = PARTITION_ITERATOR_BUFFER_SKELETON[id]
    resize!(buffer_skeleton, order_t)
    edge_set        = PARTITION_ITERATOR_BUFFER_EDGE_SET[id]
    resize!(edge_set, order_t - 1)
    edge_set_tmp    = PARTITION_ITERATOR_BUFFER_EDGE_SET_TMP[id]
    resize!(edge_set_tmp, order_t - 1)
  else
    buffer_forest_t = Vector{Int}(undef, order_t)
    level_sequence  = similar(buffer_forest_t)
    buffer_skeleton = similar(buffer_forest_t)
    edge_set        = Vector{Bool}(undef, order_t - 1)
    edge_set_tmp    = similar(edge_set)
  end

  skeleton = RootedTree(buffer_skeleton, true)
  t_forest = RootedTree(buffer_forest_t, true)
  t_temp_forest = RootedTree(level_sequence, true)
  forest = PartitionForestIterator(t_forest, t_temp_forest, edge_set_tmp)
  PartitionIterator{typeof(t), RootedTree{Int, Vector{Int}}}(
    t, forest, skeleton, edge_set, edge_set_tmp)
end

Base.IteratorSize(::Type{<:PartitionIterator}) = Base.HasLength()
Base.length(partitions::PartitionIterator) = 2^length(partitions.edge_set)
Base.eltype(::Type{PartitionIterator{TreeInput, TreeOutput}}) where {TreeInput, TreeOutput} = Tuple{PartitionForestIterator{TreeOutput}, TreeOutput}

@inline function Base.iterate(partitions::PartitionIterator)
  edge_set_value = 0
  iterate(partitions, edge_set_value)
end

@inline function Base.iterate(partitions::PartitionIterator, edge_set_value)
  edge_set_value >= length(partitions) && return nothing

  t            = partitions.t
  forest       = partitions.forest
  skeleton     = partitions.skeleton
  edge_set     = partitions.edge_set
  edge_set_tmp = partitions.edge_set_tmp

  binary_digits!(edge_set, edge_set_value)

  # Compute the partition skeleton.
  # The following is a more efficient version of
  #   skeleton = partition_skeleton(t, edge_set)
  # avoiding some allocations.
  resize!(edge_set_tmp, length(edge_set))
  copy!(edge_set_tmp, edge_set)
  unsafe_resize!(skeleton, order(t))
  copy!(skeleton, t)
  partition_skeleton!(skeleton, edge_set_tmp)

  # Compute the partition forest.
  # The following is a more efficient version of
  #   forest = partition_forest(t, edge_set)
  # avoiding some allocations and using a lazy iterator.
  resize!(edge_set_tmp, length(edge_set))
  copy!(edge_set_tmp, edge_set)
  unsafe_resize!(forest.t_temp, order(t))
  copy!(forest.t_temp, t)

  ((forest, skeleton), edge_set_value + 1)
end

# necessary for simple and convenient use since the iterates may be modified
function Base.collect(partitions::PartitionIterator{TreeInput, TreeOutput}) where {TreeInput, TreeOutput}
  iterates = Vector{Tuple{Vector{TreeOutput}, TreeOutput}}()
  sizehint!(iterates, length(partitions))
  for (forest, skeleton) in partitions
    push!(iterates, (collect(forest), copy(skeleton)))
  end
  return iterates
end

# splittings
"""
    all_splittings(t::RootedTree)

Create all splitting forests and subtrees associated to ordered subtrees of a
rooted tree `t`.

Seee also [`SplittingIterator`](@ref).

# References

Section 2.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)
"""
function all_splittings(t::RootedTree)
  node_set = zeros(Bool, order(t))
  ls = t.level_sequence
  T = eltype(ls)
  forests  = Vector{Vector{RootedTree{T, Vector{T}}}}()
  subtrees = Vector{RootedTree{T, Vector{T}}}() # ordered subtrees

  for node_set_value in 0:(2^order(t) - 1)
    binary_digits!(node_set, node_set_value)

    # Check that if a node is removed then all of its descendants are removed
    subtree_root_index = 1
    forest = Vector{RootedTree{T, Vector{T}}}()
    while subtree_root_index <= order(t)
      if node_set[subtree_root_index] == false # This node is removed
        # Find complete subtree
        subtree_last_index = _subtree_last_index(subtree_root_index, ls)

        # Check that subtree is all removed
        if !any(@view node_set[subtree_root_index:subtree_last_index])
          new_tree = rootedtree(@view ls[subtree_root_index:subtree_last_index])
          push!(forest, new_tree)
          subtree_root_index = subtree_last_index + 1
        else
          break
        end
      else
        subtree_root_index += 1
      end
    end

    if subtree_root_index == order(t) + 1
      # This is a valid ordered subtree.
      # The `level_sequence` will not automatically be a canonical representation.
      level_sequence = ls[node_set]
      subtree = rootedtree!(level_sequence)
      push!(subtrees, subtree)
      push!(forests, forest)
    end
  end

  return (; forests, subtrees)
end

"""
    SplittingIterator(t::RootedTree)

Iterator over all splitting forests and subtrees of the rooted tree `t`.
This is basically an iterator version of [`all_splittings`](@ref).

See also [`partition_forest`](@ref) and [`partition_skeleton`](@ref).

# References

Section 2.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)
"""
struct SplittingIterator{T<:RootedTree}
  t::T
  node_set::Vector{Bool}
  max_node_set_value::Int

  function SplittingIterator(t::T) where {T<:RootedTree}
    node_set = zeros(Bool, order(t))
    new{T}(t, node_set, 2^order(t) - 1)
  end
end

Base.IteratorSize(::Type{<:SplittingIterator}) = Base.SizeUnknown()
Base.eltype(::Type{SplittingIterator{T}}) where {T} = Tuple{Vector{T}, T}

@inline function Base.iterate(splittings::SplittingIterator)
  node_set_value = 0
  iterate(splittings, node_set_value)
end

@inline function Base.iterate(splittings::SplittingIterator, node_set_value)
  node_set_value > splittings.max_node_set_value && return nothing

  node_set = splittings.node_set
  t = splittings.t
  ls = t.level_sequence
  T = eltype(ls)
  forest = Vector{RootedTree{T, Vector{T}}}()

  while node_set_value <= splittings.max_node_set_value
    binary_digits!(node_set, node_set_value)

    # Check that if a node is removed then all of its descendants are removed
    subtree_root_index = 1
    empty!(forest)
    while subtree_root_index <= order(t)
      if node_set[subtree_root_index] == false # This node is removed
        subtree_last_index = _subtree_last_index(subtree_root_index, ls)

        # Check that subtree is all removed
        if !any(@view node_set[subtree_root_index:subtree_last_index])
          # If `iscanonical(t)`, the subtree starting at the root of `t`
          # is also in canonical representation. Thus, we don't need to
          # use the more expensive version
          #   push!(forest, rootedtree!(level_sequence))
          # but can use the cheaper version below.
          level_sequence = ls[subtree_root_index:subtree_last_index]
          push!(forest, RootedTree(level_sequence, iscanonical(t)))
          subtree_root_index = subtree_last_index + 1
        else
          break
        end
      else
        subtree_root_index += 1
      end
    end

    if subtree_root_index == order(t) + 1
      # This is a valid ordered subtree.
      # The `level_sequence` will not automatically be a canonical representation.
      # TODO: splittings;
      #       Decide whether canonical representations should be used. Disabling
      #       them will increase the performance.
      level_sequence = ls[node_set]
      subtree = rootedtree!(level_sequence)
      return ((forest, subtree), node_set_value + 1)
    else
      node_set_value = node_set_value + 1
    end
  end

  return nothing
end

# functions on trees

"""
    order(t::AbstractRootedTree)

The `order` of a rooted tree `t`, i.e., the length of its level sequence.
"""
order(t::AbstractRootedTree) = length(t.level_sequence)

"""
    σ(t::AbstractRootedTree)
    symmetry(t::AbstractRootedTree)

The symmetry `σ` of a rooted tree `t`, i.e., the order of the group of automorphisms
on a particular labelling (of the vertices) of `t`.

Reference: Section 301 of
- Butcher, John Charles.
  Numerical methods for ordinary differential equations.
  John Wiley & Sons, 2008.
"""
function symmetry(t::AbstractRootedTree)
  if order(t) <= 2
    return 1
  end

  # Rely on the canonical ordering to guarantee that identical subtrees are
  # next to each other.
  if !iscanonical(t)
    return symmetry(canonical_representation(t))
  end

  # Iterate over all subtrees. Since we know that the `order(t)` is at least 3,
  # there must be at least one subtree.
  subtrees = SubtreeIterator(t)

  # Unroll the `for` loop manually to be able to compare the next iterate with
  # the previous one.
  previous_subtree, state = iterate(subtrees)

  result = 1
  num_same_subtrees = 1
  iter = iterate(subtrees, state)
  while iter !== nothing
    subtree, state = iter
    if subtree == previous_subtree
      num_same_subtrees += 1
    else
      result *= factorial(num_same_subtrees) * symmetry(previous_subtree)^num_same_subtrees
      num_same_subtrees = 1
    end

    previous_subtree = subtree
    iter = iterate(subtrees, state)
  end

  result *= factorial(num_same_subtrees) * symmetry(previous_subtree)^num_same_subtrees
  return result
end

const σ = symmetry

"""
    γ(t::AbstractRootedTree)
    density(t::AbstractRootedTree)

The density `γ(t)` of a rooted tree, i.e., the product over all vertices of `t`
of the order of the subtree rooted at that vertex.

Reference: Section 301 of
- Butcher, John Charles.
  Numerical methods for ordinary differential equations.
  John Wiley & Sons, 2008.
"""
function density(t::AbstractRootedTree)
  isempty(t) && return 1

  result = order(t)
  for subtree in SubtreeIterator(t)
    result *= density(subtree)
  end
  result
end

const γ = density

"""
    α(t::AbstractRootedTree)

The number of monotonic labelings of `t` not equivalent under the symmetry group.

Reference: Section 302 of
- Butcher, John Charles.
  Numerical methods for ordinary differential equations.
  John Wiley & Sons, 2008.
"""
function α(t::AbstractRootedTree)
  div(factorial(order(t)), σ(t) * γ(t))
end

"""
    β(t::AbstractRootedTree)

The total number of labelings of `t` not equivalent under the symmetry group.

Reference: Section 302 of
- Butcher, John Charles.
  Numerical methods for ordinary differential equations.
  John Wiley & Sons, 2008.
"""
function β(t::AbstractRootedTree)
  div(factorial(order(t)), σ(t))
end

# additional representation and construction methods
"""
    t1 ∘ t2

The non-associative Butcher product of rooted trees. It is formed
by adding an edge from the root of `t1` to the root of `t2`.

Reference: Section 301 of
- Butcher, John Charles.
  Numerical methods for ordinary differential equations.
  John Wiley & Sons, 2016.
"""
function Base.:∘(t1::RootedTree, t2::RootedTree)
  offset = first(t1.level_sequence) - first(t2.level_sequence) + 1
  level_sequence = vcat(t1.level_sequence, t2.level_sequence .+ offset)
  rootedtree(level_sequence)
end

"""
    butcher_representation(t::RootedTree)

Returns the representation of `t::RootedTree` introduced by Butcher as a string.
Thus, the rooted tree consisting whose only vertex is the root itself is
represented as `τ`. The representation of other trees is defined recursively;
if `t₁, t₂, ... tₙ` are the [`subtrees`](@ref) of the rooted tree `t`, it is
represented as `t = [t₁ t₂ ... tₙ]`. If multiple subtrees are the same, their
number of occurences is written as a power.

# Examples

```jldoctest
julia> rootedtree([1, 2, 3, 2]) |> butcher_representation
"[[τ]τ]"

julia> rootedtree([1, 2, 3, 3, 2]) |> butcher_representation
"[[τ²]τ]"
```

# References

Section 300 of
- Butcher, John Charles.
  Numerical methods for ordinary differential equations.
  John Wiley & Sons, 2008.
"""
function butcher_representation(t::RootedTree, normalize::Bool=true)
  if order(t) == 0
    return "∅"
  elseif order(t) == 1
    return "τ"
  end

  result = ""
  for subtree in SubtreeIterator(t)
    result = result * butcher_representation(subtree, normalize)
  end
  result = "[" * result * "]"

  if normalize
    # normalize the result by grouping repeated occurrences of τ
    # TODO: Decide whether powers should also be used for subtrees,
    #       e.g., "[[τ]²]" instead of "[[τ][τ]]"
    #       for rootedtree([1, 2, 3, 2, 3]).
    #       Currently, powers are only used for τ.
    for n in order(t):-1:2
      n_str = string(n)
      n_str = replace(n_str, "1" => "¹")
      n_str = replace(n_str, "2" => "²")
      n_str = replace(n_str, "3" => "³")
      n_str = replace(n_str, "4" => "⁴")
      n_str = replace(n_str, "5" => "⁵")
      n_str = replace(n_str, "6" => "⁶")
      n_str = replace(n_str, "7" => "⁷")
      n_str = replace(n_str, "8" => "⁸")
      n_str = replace(n_str, "9" => "⁹")
      n_str = replace(n_str, "0" => "⁰")
      result = replace(result, "τ"^n => "τ" * n_str)
    end
  end

  return result
end

include("colored_trees.jl")
include("latexify.jl")
include("plot_recipes.jl")
include("time_integration_methods.jl")

function __init__()
  # canonical_representation!
  Threads.resize_nthreads!(CANONICAL_REPRESENTATION_BUFFER,
                           Vector{Int}(undef, BUFFER_LENGTH))

  # PartitionIterator
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_FOREST_T,
                           Vector{Int}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_FOREST_T_COLORS,
                           Vector{Bool}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_FOREST_LEVEL_SEQUENCE,
                           Vector{Int}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_FOREST_COLOR_SEQUENCE,
                           Vector{Bool}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_SKELETON,
                           Vector{Int}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_SKELETON_COLORS,
                           Vector{Bool}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_EDGE_SET,
                           Vector{Bool}(undef, BUFFER_LENGTH))
  Threads.resize_nthreads!(PARTITION_ITERATOR_BUFFER_EDGE_SET_TMP,
                           Vector{Bool}(undef, BUFFER_LENGTH))

  @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    using .Plots: Plots
    include("plots.jl")
  end

  return nothing
end

end # module


text:
export UnconstrainedOEDParam, getNumberOfConstraints

import jInv.ForwardShare.ForwardProbType

"""
type UnconstrainedOEDParam <: ForwardProbType

Parameter for OED with unconstrained lower level problem

	f(p) = argmin_f 0.5*| M(p)*f - d(p)|^2 + 0.5*|L*f|^2

STATUS: this is under development

Fields:

	Design  - design type (OEDTypeAParam)
	L       - regularization Operator
	Xr      - field for reconstructions
	Ainv    - field for linear solvers
"""
type UnconstrainedOEDParam <: ForwardProbType
   Design::AbstractOEDParam
   L
   Xr
   Ainv
end

getNumberOfConstraints(pFor::UnconstrainedOEDParam) = 0

function Base.display(pFor::UnconstrainedOEDParam)
	(m,n,nrhs) = getProblemSize(pFor.Design)
	println("--- OED.UnconstrainedOEDParam <: ForwardProbType ---")
	println("Type of Design:                    $(typeof(pFor.Design))")
	println("Size of linear system:             $m x $n")
	println("Number of right hand sides:        $nrhs")
	println("Size of regularization operator:   $(size(pFor.L))")
end

import jInv.ForwardShare.getData
function getData(p::Vector{Float64},pFor::UnconstrainedOEDParam,doClear::Bool=false)
	# get design
	M  = getM(p,pFor.Design, true)
	D  = getD(p,pFor.Design)
	L  = pFor.L

	# build and factorize normal equations
	H   = M'*M + L'*L
	pFor.Ainv = cholfact(H)

	Xr = pFor.Ainv\(M'*D)

    pFor.Xr = Xr
    return Xr,pFor
end

ext: jl


text:
#=
Problem 102
Three distinct points are plotted at random on a Cartesian plane, for which -1000 ≤ x, y ≤ 1000, such that a triangle is formed.

Consider the following two triangles:

A(-340,495), B(-153,-910), C(835,-947)
X(-175,41), Y(-421,-714), Z(574,-645)

It can be verified that triangle ABC contains the origin, whereas triangle XYZ does not.

Using triangles.txt (right click and 'Save Link/Target As...'), a 27K text file containing the co-ordinates of one thousand "random" triangles, find the number of triangles for which the interior contains the origin.

NOTE: The first two examples in the file represent the triangles in the example given above.
=#

using DataFrames

function containsOrigin(p1, p2, p3)
    # Change of basis approach...
    # Shift origin to p1, then write original origin as a linear combination. Problem gets very easy
    lin_comb = inv([(p2 - p1) (p3 - p1)]) * (-p1)
    if (lin_comb[1] > 0 && lin_comb[2] > 0 && (lin_comb[1] + lin_comb[2] < 1))
        return 1
    end
    return 0
end

df = readtable("appendices/p102_triangles.txt", header = false)
numTriangles = size(df)[1]

result = zeros(numTriangles)

for i = 1:numTriangles
    result[i] = containsOrigin(Array(df[i, 1:2])', Array(df[i, 3:4])', Array(df[i, 5:6])')
end

@show sum(result)

max_stars_repo_name: mebrunet/project-euler

max_issues_repo_name: mebrunet/project-euler

max_forks_repo_name: mebrunet/project-euler


text:
# run this on the RCI cluster in the `data/experiment_count` dir
using DrWatson
@quickactivate
using PrettyTables
using Statistics
using FileIO, BSON, DataFrames
using GenerativeAD
using CSV

# just count models, score files and training time
# number of models
data_folders = [datadir("experiments/tabular"), datadir("experiments/images"),
	datadir("experiments/images_leave-one-in")]
#allfiles = vcat(map(f->GenerativeAD.Evaluation.collect_files(f), data_folders)...)

# tabular data
df = load(datadir("evaluation/tabular_eval.bson"))[:df];
ddir = data_folders[1]
models = unique(df.modelname)
models = models[models .!= "vae+ocsvm"]
datasets = unique(df.dataset)

mcount_df_tabular = DataFrame(
	:dataset => datasets)
for model in models
	mcount_df_tabular[Symbol(model)] = zeros(Float32,length(datasets))
end
scount_df_tabular = DataFrame(
	:dataset => datasets)
for model in models
	scount_df_tabular[Symbol(model)] = zeros(Float32,length(datasets))
end

for (di,dataset) in enumerate(datasets)
	for (mi,model) in enumerate(models)
		datapath = joinpath(ddir, model, dataset)
		nseeds = length(readdir(datapath))
		files = GenerativeAD.Evaluation.collect_files(datapath)
		sfiles = filter(x->!occursin("model", x), files)
		mfiles = filter(x->occursin("model", x), files)
		mcount_df_tabular[di,mi+1] = length(mfiles)/nseeds
		scount_df_tabular[di,mi+1] = length(sfiles)/nseeds
	end
end

CSV.write("model_count_df_tabular_raw.csv", mcount_df_tabular)
CSV.write("score_count_df_tabular_raw.csv", scount_df_tabular)

mcount_df_tabular = CSV.read("model_count_df_tabular_raw.csv")
scount_df_tabular = CSV.read("score_count_df_tabular_raw.csv")

# some models dont save modelfiles
nosave_models = ["MO_GAAL", "abod", "hbos", "if", "knn", "loda", "lof", "ocsvm",
	"ocsvm_nu", "ocsvm_rbf", "pidforest"]

for model in nosave_models
	mcount_df_tabular[!,Symbol(model)] .= scount_df_tabular[!,Symbol(model)]
end
mcount_df_tabular.knn .= mcount_df_tabular.knn/3
mcount_df_tabular.vae_knn .= mcount_df_tabular.vae_knn/3

# some dont save score files
nosave_scores = ["vae_knn", "vae_ocsvm"]
for model in nosave_scores
	scount_df_tabular[!,Symbol(model)] .= mcount_df_tabular[!,Symbol(model)]
end
scount_df_tabular.vae_knn .= scount_df_tabular.vae_knn*3

# now save the fit and predict times
fit_ts = []
pred_ts = []
for model in models
	subdf = filter(r->r.modelname==model, df)
	push!(fit_ts, sum(subdf.fit_t) + sum(subdf.fs_fit_t))
	push!(pred_ts, sum(subdf.tr_eval_t) + sum(subdf.tst_eval_t) + sum(subdf.val_eval_t)
		+ sum(subdf.fs_eval_t))
end
push!(mcount_df_tabular, vcat(["fit_t"], fit_ts))
push!(mcount_df_tabular, vcat(["pred_t"], pred_ts))
#CSV.write("model_count_df_tabular.csv", mcount_df_tabular)
push!(scount_df_tabular, vcat(["fit_t"], fit_ts))
push!(scount_df_tabular, vcat(["pred_t"], pred_ts))
#CSV.write("score_count_df_tabular.csv", scount_df_tabular)
mcount_df_tabular.pidf[1:end-2] .= mcount_df_tabular.pidf[1:end-2]/3

mcount_df_tabular = CSV.read("model_count_df_tabular.csv")
scount_df_tabular = CSV.read("score_count_df_tabular.csv")

# rename datasets
include("aliases.jl")
mcount_df_tabular.dataset .= map(x->get(DATASET_ALIAS, x, x), mcount_df_tabular.dataset)
scount_df_tabular.dataset .= map(x->get(DATASET_ALIAS, x, x), scount_df_tabular.dataset)

# rename columns
sum_models = Dict(
	"aae" => ["aae", "aae_full", "aae_vamp"],
	"ocsvm" => ["ocsvm", "ocsvm_nu", "ocsvm_rbf"],
	"vae" => ["vae", "vae_full", "vae_simple"],
	"wae" => ["wae", "wae_full", "wae_vamp"]
	)
models_ordered = ["aae", "avae", "gano", "vae", "wae", "abod", "hbos", "if",
	"knn", "loda", "lof", "orbf", "osvm", "pidf", "maf", "rnvp", "sptn", "fmgn",
	"gan", "mgal", "dsvd", "vaek", "vaeo"]
mc_final_tabular = DataFrame(:dataset => copy(mcount_df_tabular.dataset))
for model in models_ordered
	ind = findfirst(values(MODEL_ALIAS) .== model)
	mc_final_tabular[!,Symbol(model)] .=
		mcount_df_tabular[!,Symbol(collect(keys(MODEL_ALIAS))[ind])]
	if model in keys(sum_models)
		for m in sum_models[model][2:end]
			mc_final_tabular[!,Symbol(model)] .+= mcount_df_tabular[!,Symbol(m)]
		end
	end
end
sc_final_tabular = DataFrame(:dataset => copy(scount_df_tabular.dataset))
for model in models_ordered
	ind = findfirst(values(MODEL_ALIAS) .== model)
	sc_final_tabular[!,Symbol(model)] .=
		scount_df_tabular[!,Symbol(collect(keys(MODEL_ALIAS))[ind])]
	if model in keys(sum_models)
		for m in sum_models[model][2:end]
			sc_final_tabular[!,Symbol(model)] .+= scount_df_tabular[!,Symbol(m)]
		end
	end
end

CSV.write("model_count_df_tabular.csv", mc_final_tabular)
CSV.write("score_count_df_tabular.csv", sc_final_tabular)

mc_final_tabular = CSV.read("model_count_df_tabular.csv")
sc_final_tabular = CSV.write("score_count_df_tabular.csv")

# now finalize them
mc_final_tabular[1:end-2, 2:end] .= round.(mc_final_tabular[1:end-2, 2:end])
smc_tabular = pretty_table(mc_final_tabular, backend = :latex)
vals=Array(mc_final_tabular[end-1,2:end])/3600
@printf("%.1e &", vals[1])
for val in vals
	@printf("%.1e & ", val)
end
for val in Array(mc_final_tabular[end,2:end])/3600
	@printf("%.1e & ", val)
end

sc_final_tabular[1:end-2, 2:end] .= round.(sc_final_tabular[1:end-2, 2:end])
ssc_tabular = pretty_table(sc_final_tabular, backend = :latex)

#################
# now the same for image models
df = load(datadir("evaluation/images_eval.bson"))[:df];
ddir = data_folders[2]
models = unique(df.modelname)
datasets = unique(df.dataset)

mcount_df_images = DataFrame(
	:dataset => datasets)
for model in models
	mcount_df_images[Symbol(model)] = zeros(Float32,length(datasets))
end
scount_df_images = DataFrame(
	:dataset => datasets)
for model in models
	scount_df_images[Symbol(model)] = zeros(Float32,length(datasets))
end

for (di,dataset) in enumerate(datasets)
	for (mi,model) in enumerate(models)
		datapath = joinpath(ddir, model, dataset)
		nac = length(readdir(datapath))
		files = GenerativeAD.Evaluation.collect_files(datapath)
		files = filter(x->occursin("seed=1", x), files)
		sfiles = filter(x->!occursin("model", x), files)
		mfiles = filter(x->occursin("model", x), files)
		mcount_df_images[di,mi+1] = length(mfiles)/nac
		scount_df_images[di,mi+1] = length(sfiles)/nac
	end
end

mcount_df_images.knn .= scount_df_images.knn/3
mcount_df_images.ocsvm .= scount_df_images.ocsvm

nosave_scores = ["vae_knn", "vae_ocsvm", "aae_ocsvm"]
for model in nosave_scores
	scount_df_images[!,Symbol(model)] .= mcount_df_images[!,Symbol(model)]
end
scount_df_tabular.vae_knn .= scount_df_tabular.vae_knn*3

CSV.write("model_count_df_images_raw.csv", mcount_df_images)
CSV.write("score_count_df_images_raw.csv", scount_df_images)

mcount_df_images = CSV.read("model_count_df_images_raw.csv")
scount_df_images = CSV.read("score_count_df_images_raw.csv")

models_ordered = ["aae", "gano", "skip", "vae", "wae", "knn", "osvm", "fano",
	"fmgn", "dsvd", "vaek", "vaeo"]
mc_final_images = DataFrame(:dataset => copy(mcount_df_images.dataset))
for model in models_ordered
	ind = findlast(values(MODEL_ALIAS) .== model)
	mc_final_images[!,Symbol(model)] .=
		mcount_df_images[!,Symbol(collect(keys(MODEL_ALIAS))[ind])]
end
sc_final_images = DataFrame(:dataset => copy(scount_df_images.dataset))
for model in models_ordered
	ind = findlast(values(MODEL_ALIAS) .== model)
	sc_final_images[!,Symbol(model)] .=
		scount_df_images[!,Symbol(collect(keys(MODEL_ALIAS))[ind])]
end

CSV.write("model_count_df_images.csv", mc_final_images)
CSV.write("score_count_df_images.csv", sc_final_images)

mc_final_images[1:end, 2:end] .= round.(mc_final_images[1:end, 2:end])
smc_images = pretty_table(mc_final_images, backend = :latex, nosubheader = true)

sc_final_images[1:end, 2:end] .= round.(sc_final_images[1:end, 2:end])
ssc_images = pretty_table(sc_final_images, backend = :latex, nosubheader = true)

fit_t_tabular = sum(df.fit_t)/3600
pred_t_tabular = (sum(df.tr_eval_t) + sum(df.tst_eval_t) + sum(df.val_eval_t))/3600
nexp_images = round(Int,sum(Array(CSV.read("model_count_df_images.csv")[1:end,2:end]))*10)
nscore_images = round(Int,sum(Array(CSV.read("score_count_df_images.csv")[1:end,2:end]))*10)

######## LOI ##############
df = load(datadir("evaluation/images_leave-one-in_eval.bson"))[:df];
ddir = data_folders[3]
models = unique(df.modelname)
datasets = unique(df.dataset)

mcount_df_images = DataFrame(
	:dataset => datasets)
for model in models
	mcount_df_images[Symbol(model)] = zeros(Float32,length(datasets))
end
scount_df_images = DataFrame(
	:dataset => datasets)
for model in models
	scount_df_images[Symbol(model)] = zeros(Float32,length(datasets))
end

for (di,dataset) in enumerate(datasets)
	for (mi,model) in enumerate(models)
		datapath = joinpath(ddir, model, dataset)
		nac = length(readdir(datapath))
		files = GenerativeAD.Evaluation.collect_files(datapath)
		files = filter(x->occursin("seed=1", x), files)
		sfiles = filter(x->!occursin("model", x), files)
		mfiles = filter(x->occursin("model", x), files)
		mcount_df_images[di,mi+1] = length(mfiles)/nac
		scount_df_images[di,mi+1] = length(sfiles)/nac
	end
end

mcount_df_images.knn .= scount_df_images.knn/3
mcount_df_images.ocsvm .= scount_df_images.ocsvm

nosave_scores = ["vae_knn", "vae_ocsvm", "aae_ocsvm"]
for model in nosave_scores
	scount_df_images[!,Symbol(model)] .= mcount_df_images[!,Symbol(model)]
end
scount_df_tabular.vae_knn .= scount_df_tabular.vae_knn*3

CSV.write("model_count_df_images_loi_raw.csv", mcount_df_images)
CSV.write("score_count_df_images_loi_raw.csv", scount_df_images)

mcount_df_images = CSV.read("model_count_df_images_loi_raw.csv")
scount_df_images = CSV.read("score_count_df_images_loi_raw.csv")

models_ordered = ["aae", "gano", "skip", "vae", "wae", "knn", "osvm", "fano",
	"fmgn", "dsvd", "vaek", "vaeo"]
mc_final_images = DataFrame(:dataset => copy(mcount_df_images.dataset))
for model in models_ordered
	ind = findlast(values(MODEL_ALIAS) .== model)
	mc_final_images[!,Symbol(model)] .=
		mcount_df_images[!,Symbol(collect(keys(MODEL_ALIAS))[ind])]
end
sc_final_images = DataFrame(:dataset => copy(scount_df_images.dataset))
for model in models_ordered
	ind = findlast(values(MODEL_ALIAS) .== model)
	sc_final_images[!,Symbol(model)] .=
		scount_df_images[!,Symbol(collect(keys(MODEL_ALIAS))[ind])]
end

CSV.write("model_count_df_images_loi.csv", mc_final_images)
CSV.write("score_count_df_images_loi.csv", sc_final_images)

mc_final_images[1:end, 2:end] .= round.(mc_final_images[1:end, 2:end])
smc_images = pretty_table(mc_final_images, backend = :latex, nosubheader = true)

sc_final_images[1:end, 2:end] .= round.(sc_final_images[1:end, 2:end])
ssc_images = pretty_table(sc_final_images, backend = :latex, nosubheader = true)

fit_t_tabular = sum(df.fit_t)/3600
pred_t_tabular = (sum(df.tr_eval_t) + sum(df.tst_eval_t) + sum(df.val_eval_t))/3600
nexp_images_loi = round(Int,sum(Array(CSV.read("model_count_df_images_loi.csv")[1:end,2:end]))*10)
nscore_images_loi = round(Int,sum(Array(CSV.read("score_count_df_images_loi.csv")[1:end,2:end]))*10)

# total times and counts
# tabular
df = load(datadir("evaluation/tabular_eval.bson"))[:df];
nexp_tabular = round(Int,sum(Array(CSV.read("model_count_df_tabular.csv")[1:end-2,2:end]))*5)
nscore_tabular = round(Int,sum(Array(CSV.read("score_count_df_tabular.csv")[1:end-2,2:end]))*5)

cdf = CSV.read("model_count_df_tabular.csv")
sdf = CSV.read("score_count_df_tabular.csv")
multi_models = ["aae", "wae", "vae", "knn", "pidf"]
scales = []
for m in multi_models
	push!(scales, sum(sdf[1:end-2,Symbol(m)])/sum(cdf[1:end-2,Symbol(m)]))
end
fit_t_tabular = 0
for model in unique(df.modelname)
	ft = sum(filter(r->r.modelname==model, df)[!,:fit_t])
	i = findfirst(map(x->occursin(x,model),multi_models))
	if i != nothing && !(occursin("ocsvm", model))
		global fit_t_tabular += ft/scales[i]
	else
		global fit_t_tabular += ft
	end
end
fit_t_tabular = fit_t_tabular/3600/24
pred_t_tabular = (sum(df.tr_eval_t) + sum(df.tst_eval_t) + sum(df.val_eval_t))/3600/24

# images
df = load(datadir("evaluation/images_eval.bson"))[:df];
nexp_images = round(Int,sum(Array(CSV.read("model_count_df_images.csv")[1:end,2:end]))*10)
nscore_images = round(Int,sum(Array(CSV.read("score_count_df_images.csv")[1:end,2:end]))*10)

cdf = CSV.read("model_count_df_images.csv")
sdf = CSV.read("score_count_df_images.csv")
multi_models = ["aae", "gano", "skip", "wae", "vae", "knn"]
scales = []
for m in multi_models
	push!(scales, sum(sdf[1:end-2,Symbol(m)])/sum(cdf[1:end-2,Symbol(m)]))
end
multi_models2 = ["aae", "GANo", "skip", "wae", "vae", "knn"]
fit_t_images = 0
for model in unique(df.modelname)
	ft = sum(filter(r->r.modelname==model, df)[!,:fit_t])
	i = findfirst(map(x->occursin(x,model),multi_models2))
	if i != nothing && !(occursin("vae_knn", model))
		global fit_t_images += ft/scales[i]
	else
		global fit_t_images += ft
	end
end
fit_t_images = fit_t_images/3600/24
pred_t_images = (sum(df.tr_eval_t) + sum(df.tst_eval_t) + sum(df.val_eval_t))/3600/24

#images LOI
df = load(datadir("evaluation/images_leave-one-in_eval.bson"))[:df];
nexp_images_loi = round(Int,sum(Array(CSV.read("model_count_df_images_loi.csv")[1:end,2:end]))*10)
nscore_images_loi = round(Int,sum(Array(CSV.read("score_count_df_images_loi.csv")[1:end,2:end]))*10)

cdf = CSV.read("model_count_df_images_loi.csv")
sdf = CSV.read("score_count_df_images_loi.csv")
multi_models = ["aae", "gano", "skip", "wae", "vae", "knn"]
scales = []
for m in multi_models
	push!(scales, sum(sdf[1:end-2,Symbol(m)])/sum(cdf[1:end-2,Symbol(m)]))
end
multi_models2 = ["aae", "GANo", "skip", "wae", "vae", "knn"]
fit_t_images_loi = 0
for model in unique(df.modelname)
	ft = sum(filter(r->r.modelname==model, df)[!,:fit_t])
	i = findfirst(map(x->occursin(x,model),multi_models2))
	if i != nothing && !(occursin("vae_knn", model))
		global fit_t_images_loi += ft/scales[i]
	else
		global fit_t_images_loi += ft
	end
end
fit_t_images_loi = fit_t_images_loi/3600/24
pred_t_images_loi = (sum(df.tr_eval_t) + sum(df.tst_eval_t) + sum(df.val_eval_t))/3600/24

total_nmodels = nexp_images + nexp_images_loi + nexp_tabular
total_nscores = nscore_images + nscore_images_loi + nscore_tabular
total_fit_t = fit_t_tabular + fit_t_images + fit_t_images_loi
total_pred_t = pred_t_tabular + pred_t_images + pred_t_images_loi


text:
### A Pluto.jl notebook ###
# v0.16.0

using Markdown
using InteractiveUtils

# ╔═╡ 43a2b752-e6d2-11eb-0dbb-059487ed50e4
begin
    using Pkg
    Pkg.activate("..")
    using PlutoUI, Colors, Javis, Animations, StatsBase, Fontconfig
    import Luxor
end

# ╔═╡ 21ed956a-2b0e-435a-977d-202ed102d2d0
begin
    red = colorant"#ea5153"
    blue = colorant"#212c52"
    # grey = colorant"#323232"
    grey = colorant"#7a8097"

    n_frames = 200
    left_margin = -230
    right_margin = 230
    top_margin = -125
    bottom_margin = 125
    step_len = Int(n_frames / 5)
    right_end = right_margin + 50
    bottom_end = bottom_margin + 50
    global_margin = 15
    sizeof4 = 150
    sizeofletters = 200

    function ground(args...)
        background("white")
        sethue("black")
    end

    function convex_hull(p1, p2, p3, frame, n_frames, idx)

        if frame < idx * step_len
            x1, x2 = p1[1], p2[1]
            y1, y2 = p1[2], p2[2]
            Point(
                x1 * (1 - frame % step_len / step_len) + x2 * (frame % step_len / step_len),
                y1 * (1 - frame % step_len / step_len) + y2 * (frame % step_len / step_len),
            )
        else
            x1, x2 = p2[1], p3[1]
            y1, y2 = p2[2], p3[2]
            Point(
                # x1 * (1 - frame%step_len/step_len) + x2 * (frame%step_len/step_len),
                x1 * (1 - (frame - idx * step_len) / n_frames) +
                x2 * ((frame - idx * step_len) / n_frames),
                # y1 * (1 - frame%step_len/step_len) + y2 * (frame%step_len/step_len),
                y1 * (1 - (frame - idx * step_len) / n_frames) +
                y2 * ((frame - idx * step_len) / n_frames),
            )
        end
    end

    function from_to_end(idx, step_len)
        return (((idx - 1) * step_len) + 1):n_frames
    end

    function gen_text(
        args...;
        size,
        body,
        x = 0,
        y = 0,
        color = red,
        valign = :middle,
        halign = :left,
    )

        fontface("Barlow")
        fontsize(size)
        setfont("Barlow", size)
        sethue(color)
        # settext(body, Point(x, y); valign=string(valign), halign=string(halign))
        text(body, Point(x, y); valign = valign, halign = halign)
    end

    my_video = Video(2 * right_end, 2 * bottom_end)

    Background(1:n_frames, ground)

    converter = [1, 2, 3, 4]

    letters = ["B", "4", "D", "S"]
    locations = [
        left_margin - 16 + global_margin,
        left_margin - 16 + global_margin + (sizeofletters ÷ 5) * 3 + global_margin,
        (right_margin - left_margin) / 4 - global_margin,
        right_margin - global_margin,
    ]
    letters = map(zip(1:4, letters, locations)) do (idx, letter, loc)

        Object(
            from_to_end(converter[idx], step_len),
            (args...) -> gen_text(
                args...,
                size = size = letter != "4" ? sizeofletters : sizeof4,
                body = letter,
                x = loc,
                y = letter != "4" ? top_margin + global_margin :
                    top_margin + sizeofletters - sizeof4,
                color = (letter != "4" ? blue : red),
                halign = idx <= 2 ? :left : :right,
                valign = :top,
            ),
        )
    end

    map(enumerate(letters)) do (idx, l)
        act!(l, Action(0:step_len, linear(), appear(:fade)))
    end

    words = ["Business", "Engineering", "Data", "Science"]
    words = map(enumerate(words)) do (idx, word)
        ob1 = Object(
            from_to_end(converter[idx], step_len),
            (args...) -> gen_text(
                args...,
                size = 45,
                body = word,
                x = if idx == 1
                    left_margin - 3 + global_margin
                elseif idx == 3
                    left_margin + 150
                else
                    right_margin - global_margin
                end,
                y = idx <= 2 ? bottom_margin - 70 : bottom_margin - global_margin,
                color = idx != 2 ? red : blue,
                halign = idx % 2 == 1 ? :left : :right,
                valign = idx <= 2 ? :middle : :bottom,
            ),
            keep = true,
        )
        if idx == 2
            ob2 = Object(
                from_to_end(converter[idx], step_len),
                (args...) -> gen_text(
                    args...,
                    size = 50,
                    body = "for",
                    x = left_margin + global_margin,
                    y = bottom_margin - global_margin,
                    color = blue,
                    halign = :left,
                    valign = :bottom,
                ),
                keep = true,
            )
            [ob1, ob2]
        else
            ob1
        end
    end

    map(enumerate(words)) do (idx, w)
        act!(w, Action(0:step_len, linear(), appear(:fade)))
    end

    map(1:4) do idx
        if idx == 1
            p1 = Point(left_margin, bottom_margin)
            p2 = Point(left_margin, top_margin)
            p3 = Point(left_margin, -bottom_end)
        elseif idx == 2
            p1 = Point(left_margin, top_margin)
            p2 = Point(right_margin, top_margin)
            p3 = Point(right_end, top_margin)
        elseif idx == 3
            p1 = Point(right_margin, top_margin)
            p2 = Point(right_margin, bottom_margin)
            p3 = Point(right_margin, bottom_end)
        elseif idx == 4
            p1 = Point(right_margin, bottom_margin)
            p2 = Point(left_margin, bottom_margin)
            p3 = Point(-right_end, bottom_margin)
        end
        Object(
            from_to_end(converter[idx], step_len),
            (args...) -> line(p1, convex_hull(p1, p2, p3, args[3], step_len, idx), :stroke),
            keep = true,
        )
    end

    render(my_video, pathname = "../output/B4DS.gif")
end

# ╔═╡ baafcd46-e0a7-40a1-84d2-d51b82669794
let
    function ground(args...)
        background("white") # canvas background
        sethue("black") # pen color
    end

    function object(p = O, color = "black")
        sethue(color)
        circle(p, 25, :fill)
        return p
    end

    function path!(points, pos, color)
        sethue(color)
        push!(points, pos) # add pos to points
        circle.(points, 2, :fill) # draws a circle for each point using broadcasting
    end

    function connector(p1, p2, color)
        sethue(color)
        line(p1, p2, :stroke)
    end

    myvideo = Video(500, 500)
    path_of_red = []
    Background(1:70, ground)
    red_ball = Object(1:70, (args...) -> object(O, "red"), Point(100, 0))
    act!(red_ball, Action(anim_rotate_around(2π, O)))
    Object(1:70, (args...) -> path!(path_of_red, pos(red_ball), "red"))
    path_of_blue = []
    blue_ball = Object(1:70, (args...) -> object(O, "blue"), Point(200, 80))
    act!(blue_ball, Action(anim_rotate_around(2π, 0.0, red_ball)))
    Object(1:70, (args...) -> path!(path_of_blue, pos(blue_ball), "blue"))
    Object(1:70, (args...) -> connector(pos(red_ball), pos(blue_ball), "black"))

    render(myvideo; pathname = "../output/circle.gif")
end

# ╔═╡ d2e53571-5e30-4979-a1a6-39c33c98481a
md"## Utility"

# ╔═╡ e1f05461-f356-4f09-a02f-7a864acc3361
function egg(radius, action = :none)
    A, B = [Point(x, 0) for x in [-radius, radius]]
    nints, C, D = intersectionlinecircle(Point(0, -2radius), Point(0, 2radius), A, 2radius)

    flag, C1 = intersectionlinecircle(C, D, O, radius)
    nints, I3, I4 = intersectionlinecircle(A, C1, A, 2radius)
    nints, I1, I2 = intersectionlinecircle(B, C1, B, 2radius)

    if distance(C1, I1) < distance(C1, I2)
        ip1 = I1
    else
        ip1 = I2
    end
    if distance(C1, I3) < distance(C1, I4)
        ip2 = I3
    else
        ip2 = I4
    end

    newpath()
    arc2r(B, A, ip1, :path)
    arc2r(C1, ip1, ip2, :path)
    arc2r(A, ip2, B, :path)
    arc2r(O, B, A, :path)
    closepath()

    do_action(action)
end

# ╔═╡ Cell order:
# ╟─baafcd46-e0a7-40a1-84d2-d51b82669794
# ╟─21ed956a-2b0e-435a-977d-202ed102d2d0
# ╟─d2e53571-5e30-4979-a1a6-39c33c98481a
# ╟─43a2b752-e6d2-11eb-0dbb-059487ed50e4
# ╟─e1f05461-f356-4f09-a02f-7a864acc3361


text:
using CommonLH

# Must be outside of testset.
@common_fields test_fields begin
    x :: Float64
    # Comment
    y :: String
    # Parametric type
    z :: Vector{Int}
end

struct Foo
    @test_fields
    x1 :: Int
end

function struct_test()
    @testset "Structs" begin
        f = Foo(1.0, "y", [1,2], 3);
        @test f.x == 1.0;
        @test f.y == "y";
        @test f.z == [1,2];
        @test f.x1 == 3;
    end
end

@testset "Structs" begin
    struct_test();
end

# -----------


text:
"""
Wrapper module for [*tectonic*](https://github.com/tectonic-typesetting/tectonic).
"""
module Tectonic

export binary, tectonic, version, Biber

using Pkg.Artifacts

"""
    tectonic(f)

Call function `f` with the *tectonic* binary path as argument.

```julia
tectonic() do bin
    run(`\$bin file.tex`)
end
```

All required setup and clean-up to run the binary correctly is done automatically.

!!! note

    Currently nothing is needed to be done to run the binary, but this may
    change in future and this provides a future-proof interface to the package.
"""
tectonic(f) = f(binary())

"""
Version of *tectonic* provided by this package.
"""
version() = VersionNumber(lstrip(!isnumeric, tectonic(bin -> readchomp(`$bin --version`))))

binary() = joinpath(artifact"tectonic_bin", "tectonic$(Sys.iswindows() ? ".exe" : "")")

"""
Path to *tectonic* binary.

!!! warning

    This should not be used to run the program directly, please use the
    [`tectonic`](@ref) function instead since it performs the required setup
    and clean up steps needed.
"""
binary()

include("Biber.jl")

end # module


text:
# maximal
"""
    set_maximal_state(mechanism, z)

    set the maximal state of a mechanism

    mechanism: Mechanism
    z: state
"""
function set_maximal_state!(mechanism::Mechanism, z::AbstractVector)
    off = 0
    for body in mechanism.bodies
        x2, v15, q2, ϕ15 = unpack_data(z[off+1:end]); off += 13
        q2 = Quaternion(q2...)
        body.state.v15 = v15
        body.state.ϕ15 = ϕ15
        body.state.x2 = x2
        body.state.q2 = q2
		initialize_state!(mechanism) # set x1, q1 and zeroes out F2 τ2
    end
	# warm-start solver
	for body in mechanism.bodies
		set_velocity_solution!(body)
	end
end

function initialize_state!(mechanism::Mechanism)
    for body in mechanism.bodies initialize_state!(body, mechanism.timestep) end
end

# inputs
"""
    set_input(mechanism, u)

    set input for each joint in mechanism

    mechanism: Mechanism
    u: input
"""
function set_input!(mechanism::Mechanism{T}, u::AbstractVector) where T
	joints = mechanism.joints
	# set the controls in the equality constraints
	off = 0
	for joint in joints
		nu = input_dimension(joint)
		set_input!(joint, SVector{nu,T}(u[off .+ (1:nu)]))
		off += nu
	end
	# apply the controls to each body's state
	for joint in joints
		input_impulse!(joint, mechanism)
	end
end

function set_input!(mechanism::Mechanism, dict::Dict)
    for (id, joint) in pairs(mechanism.joints)
        set_input!(joint, dict[id])
    end
end

# minimal
"""
    set_minimal_state(mechanism, y)

    set the maximal state of a mechanism

    mechanism: Mechanism
    y: state
"""
function set_minimal_state!(mechanism::Mechanism, y::AbstractVector)
    z = minimal_to_maximal(mechanism, y)
    set_maximal_state!(mechanism, z)
end

function set_minimal_coordinates!(mechanism::Mechanism, dict::Dict)
    for (id, joint) in pairs(mechanism.joints)
        set_minimal_coordinates!(mechanism, joint, dict[id])
    end
end

function set_minimal_velocities!(mechanism::Mechanism, dict::Dict)
    for (id, joint) in pairs(mechanism.joints)
        set_minimal_velocities!(mechanism, joint, dict[id])
    end
end

# velocity
"""
    zero_velocity!(mechanism)

    set all mechanism body velocities to zero

    mechanism: Mechanism
"""
function zero_velocity!(mechanism::Mechanism)
    for (i, body) in enumerate(mechanism.bodies)
        try
            set_maximal_velocities!(body, v=zeros(3), ω=zeros(3))
            set_previous_configuration!(body, mechanism.timestep)
        catch
            nothing
        end
    end
end

# springs + dampers
function set_spring_offset!(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, x::AbstractVector) where {T,Nn,Ne,Nb,Ni}
	off = 0
	for id in mechanism.root_to_leaves
		(id > Ne) && continue # only treat joints
		joint = mechanism.joints[id]
        N̄ = 3 - length(joint)
        joint.spring_offset = x[off .+ (1:N̄)]
        off += 2N̄
	end
	return nothing
end


text:
using Gurobi

const GRB_OPTIMIZER = Gurobi.Optimizer(OutputFlag=0)
MOI.set(GRB_OPTIMIZER, MOI.Silent(), true)
const GRB_CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const GRB_CACHED = MOIU.CachingOptimizer(GRB_CACHE, GRB_OPTIMIZER)
const GRB_BRIDGED = MOI.Bridges.full_bridge_optimizer(GRB_CACHED, Float64)

push!(solvers, (opt = GRB_BRIDGED, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_quad, (opt = GRB_BRIDGED, mode = BilevelJuMP.SOS1Mode()))


text:
if Sys.iswindows()
    download("https://wj32.org/wp/download/1455/", joinpath(@__DIR__, "EmptyStandbyList.exe"))
end

max_stars_repo_name: davidanthoff/csv-comparison

max_issues_repo_name: davidanthoff/csv-comparison

max_forks_repo_name: davidanthoff/csv-comparison


text:
pars(d) = pars(d, true) + pars(d, false)
freepars(d) = pars(d, true)
fixedpars(d) = pars(d, false)
ispar(d, s::Symbol) = (s ∈ keys(pars(d)))
isfreepar(d, s::Symbol) = (s ∈ keys(freepars(d)))
#

abstract type AbstractFunctionWithParameters end
#
nfreepars(d::AbstractFunctionWithParameters) = length(freepars(d))
#
func(d::AbstractFunctionWithParameters, x::AbstractArray; p=freepars(d)) = func.(Ref(d), x; p=p+fixedpars(p))
func(d::AbstractFunctionWithParameters, x::AbstractRange; p=freepars(d)) = func.(Ref(d), x; p=p+fixedpars(p))
#
(d::AbstractFunctionWithParameters)(x; p=freepars(d)) = func(d,x;p)
const ArrayOrRange = Union{AbstractArray,AbstractRange}
(d::AbstractFunctionWithParameters)(x, v::ArrayOrRange) = d(x; p=v2p(v,d))

# methods that call `updatevalueorflag`

updateisfree(d::AbstractFunctionWithParameters, s::Symbol, isfree::Bool) =
    ispar(d,s) ? updatevalueorflag(d,s,isfree) : d
updatevalue( d::AbstractFunctionWithParameters, s::Symbol, v) =
    ispar(d,s) ? updatevalueorflag(d, s, isfreepar(d, s), v) : d
#
# singular
updatepar( d::AbstractFunctionWithParameters, s::Symbol, v) = updatevalue(d, s, v)
releasepar(d::AbstractFunctionWithParameters, s::Symbol) = updateisfree(d, s, true )
fixpar(    d::AbstractFunctionWithParameters, s::Symbol) = updateisfree(d, s, false)
fixpar(    d::AbstractFunctionWithParameters, s::Symbol, v) = updateisfree(updatepar(d, s, v), s, false)

# plural
const SymbolSequenceType = Union{AbstractVector{Symbol}, Tuple{Vararg{Symbol}}}
function fixpars(d::AbstractFunctionWithParameters, sequence::SymbolSequenceType)
    dnew = d
    for s in sequence
        dnew = fixpar(dnew, s)
    end
    return dnew
end

function releasepars(d::AbstractFunctionWithParameters, sequence::SymbolSequenceType)
    dnew = d
    for s in sequence
        dnew = releasepar(dnew, s)
    end
    return dnew
end

function updatepars(d::AbstractFunctionWithParameters, sequence::NamedTuple)
    dnew = d
    for (s,v) in zip(keys(sequence), sequence)
        dnew = updatepar(dnew, s, v)
    end
    return dnew
end
fixpars(d::AbstractFunctionWithParameters, sequence::NamedTuple) =
    fixpars(updatepars(d, sequence), keys(sequence))

# Number <: AbstractFunctionWithParameters
pars(d::Number, isfree::Bool) = ∅
func(d::Number, x::NumberOrTuple; p=∅) = d

# Function <: AbstractFunctionWithParameters
pars(f::Function, isfree::Bool) = ∅
func(f::Function, x::NumberOrTuple; p=∅) = f(x)

# # to be implemented for F <: AbstractFunctionWithParameters
# func(d::F, x::NumberOrTuple; p)
# pars(d::F, isfree::Bool)
# updatevalueorflag( d::F, s::Symbol, isfree::Bool, v=getproperty(pars(d),s))

# default method assumes that the d has d.p, and is a single arg
pars(d::AbstractFunctionWithParameters, isfree::Bool) = pars(d.p, isfree)
updatevalueorflag(d::AbstractFunctionWithParameters, s::Symbol, isfree::Bool, v=getproperty(pars(d),s)) =
    typeof(d)(updatevalueorflag(d.p, s, isfree, v))
#

# wrap julia function
@with_kw struct FunctionWithParameters{T} <: AbstractFunctionWithParameters
    f::Function
    p::T
end
FunctionWithParameters(f;p) = FunctionWithParameters(;f,p)

# two methods to be defined
func(d::FunctionWithParameters, x::NumberOrTuple; p=freepars(d)) = d.f(x; p=p+fixedpars(d))
pars(d::FunctionWithParameters, isfree::Bool) = pars(d.p, isfree)
updatevalueorflag(d::FunctionWithParameters, s::Symbol, isfree::Bool, v=getproperty(pars(d),s)) =
    FunctionWithParameters(;f=d.f, p=updatevalueorflag(d.p,s,isfree,v))
#

"""
    @makefuntype MyPDF(x;p) = unnormdensity(x, p.a, p.b)

    Expected form of the expression is `f(x;p)` on the left
"""
macro makefuntype(ex)
    #
    fpx = ex.args[1].args
    name = fpx[1]
    p = fpx[2].args[1]
    (p != :p) && error("expected format f(x;p) = ... " )
    x = fpx[3]
    body = ex.args[2]
    #
    quote
        struct $name{T} <: AbstractFunctionWithParameters
            p::T
        end
        $(esc(name))(;p) = $(esc(name))(p)
        import AlgebraPDF: func
        $(esc(:func))(d::$(esc(name)), $(esc(x))::NumberOrTuple; p=$(esc(:pars))(d)) = $(esc(body))
    end
end

# parameters conversion
v2p(v,d::AbstractFunctionWithParameters) = NamedTuple{keys(freepars(d))}(v)
p2v(p,d::AbstractFunctionWithParameters) = [getproperty(p, k) for k in keys(freepars(d))]
p2v(  d::AbstractFunctionWithParameters) = p2v(freepars(d), d)

# single-argument lambda-function with fixed parameters
noparsf(d::AbstractFunctionWithParameters; p=pars(d)) = (x;kw...)->func(d,x;p=p)


text:
notary = Notary()
cypher = CypherSuite(notary)

# Somewhere far far away

mixkey = newsigner(notary)

MIX_PORT = 2000
MIX_ID = id(mixkey)

mix = Mix(MIX_PORT, MIX_ID, cypher)
@async SynchronicBallot.serve(mix, mixkey)

sleep(1.)

# Using the mixer

guardiankey = newsigner(notary)
transactions = Trigger(TransactionVector(), t -> println(typeof(t.document)))
tlog = BraidChain(transactions, Guardian(id(guardiankey)), notary)

registrator = newsigner(notary)
reg = Registrator(id(registrator))
push!(tlog, Transaction(reg, guardiankey))

braider = newsigner(notary)

push!(tlog, Transaction(Braider(id(braider)), guardiankey))

gk = GateKeeper(2001, id(braider), cypher, mix)
sb = BraiderConfig(UInt8(3), UInt8(64), gk)
bo = BraidOfficer(sb, tlog)
botask = @async BraidChains.serve!(bo, braider)
sleep(1.)

### Registration of members ###

member1key = newsigner(notary)
member2key = newsigner(notary)
member3key = newsigner(notary)

member1 = Member(id(member1key))
member2 = Member(id(member2key))
member3 = Member(id(member3key))

# on the server
push!(tlog, Transaction(member1, registrator))
push!(tlog, Transaction(member2, registrator))
push!(tlog, Transaction(member3, registrator))

# back to user

pseudonym1 = newsigner(notary)
pseudonym2 = newsigner(notary)
pseudonym3 = newsigner(notary)

@show istaskstarted(botask)

@sync begin
    @async braid(member1key, pseudonym1, bo)
    @async braid(member2key, pseudonym2, bo)
    @async braid(member3key, pseudonym3, bo)
end

### Some tests

@test pseudonyms(tlog, 5) == Set([id(member1key), id(member2key), id(member3key)]) == tlog.state.members
@test pseudonyms(tlog) == Set([id(pseudonym1), id(pseudonym2), id(pseudonym3)])
@test tlog.state.braiders == Set([id(braider)])
@test tlog.state.registrators == Set([id(registrator)])

# ### TODO #####
# kc = KeyChain(notary, user)
# braid!(braidconf, tlog, kc)
# ###############

# ### TODO ######
# # perhaps kc[n] could select the key for nth row
# sign(data,n,kc) # signs data with coresponding pseudonym tagged by nth row in the TransactionLog
# #########


text:
import Base: FastMath

import SpecialFunctions

const MATH_INTRINSICS = GCNIntrinsic[]

for jltype in (Float16, Float32, Float64)
    for intrinsic in (
        :sin, :cos, :tan, :asin, :acos, :atan, :atan2,
        :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
        :sinpi, :cospi, :tanpi,
        :asinpi, :acospi, :atanpi, :atan2pi,
        :sqrt, :rsqrt, :cbrt, :rcbrt, :recip,
        :log, :log2, :log10, :log1p, :logb, :ilogb,
        :exp, :exp2, :exp10, :expm1,
        :erf, :erfinv, :erfc, :erfcinv, :erfcx,
        # TODO: :brev, :clz, :ffs, :byte_perm, :popc,
        :isnormal, :nearbyint, :nextafter,
        :lgamma)
        # FIXME: :lgamma_r segfaults on GPU

        if intrinsic == :sin && jltype == Float16
            continue # FIXME: https://github.com/JuliaGPU/AMDGPU.jl/issues/177
            push!(MATH_INTRINSICS, GCNIntrinsic(intrinsic, inp_args=(jltype,), out_arg=jltype, isbroken=true))
        end

        push!(MATH_INTRINSICS, GCNIntrinsic(intrinsic, inp_args=(jltype,), out_arg=jltype))
    end

    push!(MATH_INTRINSICS, GCNIntrinsic(:besselj0, :j0; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:besselj1, :j1; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:bessely0, :y0; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:bessely1, :y1; inp_args=(jltype,), out_arg=jltype))

    push!(MATH_INTRINSICS, GCNIntrinsic(:sin_fast, :native_sin; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:cos_fast, :native_cos; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:sqrt_fast, :native_sqrt; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:rsqrt_fast, :native_rsqrt; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:recip_fast, :native_recip; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log_fast, :native_log; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log2_fast, :native_log2; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:log10_fast, :native_log10; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp_fast, :native_exp; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp2_fast, :native_exp2; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:exp10_fast, :native_exp10; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:abs, :fabs; inp_args=(jltype,), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:gamma, :tgamma; inp_args=(jltype,), out_arg=jltype))
    # TODO: abs(::Union{Int32,Int64})

    # Multi-argument functions
    push!(MATH_INTRINSICS, GCNIntrinsic(:^, :pow; inp_args=(jltype,jltype), out_arg=jltype))
    push!(MATH_INTRINSICS, GCNIntrinsic(:^, :pown; inp_args=(jltype,Union{UInt32,Int32}), out_arg=jltype))
    # TODO: :sincos, :frexp, :ldexp, :copysign,
    #push!(MATH_INTRINSICS, GCNIntrinsic(:ldexp; inp_args=(jltype,), out_arg=(jltype, Int32), isinverted=true))

    push!(MATH_INTRINSICS, GCNIntrinsic(:hypot; inp_args=(jltype,jltype), out_arg=jltype))

    # Multi-output functions
    push!(MATH_INTRINSICS, GCNIntrinsic(:sincospi; inp_args=(jltype,), out_arg=jltype, isbroken=true))
end

for jltype in (Float32, Float64)
    push!(MATH_INTRINSICS, GCNIntrinsic(:isfinite; inp_args=(jltype,), out_arg=Int32, tobool=true))
    push!(MATH_INTRINSICS, GCNIntrinsic(:isinf; inp_args=(jltype,), out_arg=Int32, tobool=true))
    push!(MATH_INTRINSICS, GCNIntrinsic(:isnan; inp_args=(jltype,), out_arg=Int32, tobool=true))
    push!(MATH_INTRINSICS, GCNIntrinsic(:signbit; inp_args=(jltype,), out_arg=Int32))
end

for intr in MATH_INTRINSICS
    inp_vars = [gensym() for _ in 1:length(intr.inp_args)]
    inp_expr = [:($(inp_vars[idx])::$arg) for (idx,arg) in enumerate(intr.inp_args)]
    libname = Symbol("__$(intr.roclib)_$(intr.rocname)_$(intr.suffix)")
    mod = if isdefined(Base, intr.jlname)
        Base
    elseif isdefined(FastMath, intr.jlname)
        FastMath
    elseif isdefined(SpecialFunctions, intr.jlname)
        SpecialFunctions
    else
        nothing
    end
    ex = quote
        y = _intr($(Val(libname)), $(intr.out_arg), $(inp_expr...))
        y = $(intr.isinverted ? :(1-y) : :y)
        y = $(intr.tobool ? :(y != zero(y)) : :y)
        return y
    end
    if mod !== nothing
        @eval @device_override function $mod.$(intr.jlname)($(inp_expr...))
            $ex
        end
    else
        @eval @device_function function $(intr.jlname)($(inp_expr...))
            $ex
        end
    end
end

# ocml_sin seems broken for F16 (see #177)
@device_override Base.sin(x::Float16) = sin(Float32(x))

@device_override Base.hypot(x::T, y::T) where T <: Integer = hypot(float(x), float(y))

# Non-matching types
@device_override @inline function Base.:(^)(x::T, y::S) where {T<:Union{Float16, Float32,Float64}, S<:Union{Int64,UInt64}}
    y == -1 && return inv(x)
    y == 0 && return one(x)
    y == 1 && return x
    y == 2 && return x*x
    y == 3 && return x*x*x
    if S === Int64 && typemin(Int32) <= y <= typemax(Int32)
        return x ^ unsafe_trunc(Int32, y)
    elseif S === UInt64 && typemin(UInt32) <= y <= typemax(UInt32)
        return x ^ unsafe_trunc(UInt32, y)
    end
    x ^ T(y)
end
@device_override Base.:(^)(x::Integer, p::T) where T<:Union{Float16, Float32, Float64} = T(x) ^ p

@device_override @inline function Base.:(^)(x::Integer, p::Integer)
    p < 0 && throw("Negative integer power not supported")
    p == 0 && return one(x)
    p == 1 && return x
    p == 2 && return x*x
    p == 3 && return x*x*x

    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end


text:
using SymPy: N, Sym, diff, expand, solve, subs

"""
Quadratic Solov'ev equilibrium in (R,Z,phi) coordinates.
Based on McCarthy, Physics of Plasmas 6, 3554, 1999.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `a`, b`: free constants
"""
struct SolovevQuadratic{T <: Number} <: AbstractSolovevEquilibrium
    name::String
    R₀::T
    B₀::T
    a::T
    b::T

    function SolovevQuadratic{T}(R₀::T, B₀::T, a::T, b::T) where T <: Number
        new("QuadraticSolovevEquilibrium", R₀, B₀, a, b)
    end
end

function SolovevQuadratic(R₀::T, B₀::T, a::T, b::T) where T <: Number
    SolovevQuadratic{T}(R₀, B₀, a, b)
end

function Base.show(io::IO, equ::SolovevQuadratic)
    print(io, "Quadratic Solovev Equilibrium with\n")
    print(io, "  R₀ = ", equ.R₀, "\n")
    print(io, "  B₀ = ", equ.B₀, "\n")
    print(io, "  a  = ", equ.a, "\n")
    print(io, "  b  = ", equ.b)
end

function A₁(x::AbstractArray{T,1}, equ::SolovevQuadratic) where {T <: Number}
    zero(T)
end

function A₂(x::AbstractArray{T,1}, equ::SolovevQuadratic) where {T <: Number}
    zero(T)
end

function A₃(x::AbstractArray{T,1}, equ::SolovevQuadratic) where {T <: Number}
    - equ.a * R(x,equ)^4 / 8 - equ.b * Z(x,equ)^2 / 2
end

macro solovev_equilibrium_quadratic(R₀, B₀, a, b)
    generate_equilibrium_code(SolovevQuadratic(R₀, B₀, a, b); output=false)
end


text:
# This file was generated, do not modify it. # hide
X = select(carseats, Not([:Sales, :High]))
y = carseats.High;