text:
""" The category of finite sets and functions, and its skeleton.
"""
module FinSets
export FinSet, FinFunction, FinDomFunction, TabularSet, TabularLimit,
  force, is_indexed, preimage,
  JoinAlgorithm, SmartJoin, NestedLoopJoin, SortMergeJoin, HashJoin,
  SubFinSet, SubOpBoolean

using AutoHashEquals
using DataStructures: OrderedDict, IntDisjointSets, union!, find_root!
using Reexport
import StaticArrays
using StaticArrays: StaticVector, SVector, SizedVector, similar_type
import Tables, PrettyTables

@reexport using ..Sets
using ...GAT, ...Theories, ...CSetDataStructures, ...Graphs
using ..FinCats, ..FreeDiagrams, ..Limits, ..Subobjects
import ...Theories: Ob, meet, ∧, join, ∨, top, ⊤, bottom, ⊥
import ..Categories: ob, hom, dom, codom, compose, id, ob_map, hom_map
import ..FinCats: force, ob_generators, hom_generators, graph, is_discrete
using ..FinCats: dicttype
import ..Limits: limit, colimit, universal, pushout_complement,
  can_pushout_complement
import ..Subobjects: Subobject, SubobjectLattice
using ..Sets: IdentityFunction, SetFunctionCallable

# Finite sets
#############

""" Finite set.

A finite set has abstract type `FinSet{S,T}`. The second type parameter `T` is
the element type of the set and the first parameter `S` is the collection type,
which can be a subtype of `AbstractSet` or another Julia collection type. In
addition, the skeleton of the category **FinSet** is the important special case
`S = Int`. The set ``{1,…,n}`` is represented by the object `FinSet(n)` of type
`FinSet{Int,Int}`.
"""
abstract type FinSet{S,T} <: SetOb{T} end

FinSet(set::FinSet) = set

""" Finite set of the form ``{1,…,n}`` for some number ``n ≥ 0``.
"""
@auto_hash_equals struct FinSetInt <: FinSet{Int,Int}
  n::Int
end

FinSet{Int,Int}(n::Int) = FinSetInt(n)
FinSet(n::Int) = FinSetInt(n)

Base.iterate(set::FinSetInt, args...) = iterate(1:set.n, args...)
Base.length(set::FinSetInt) = set.n
Base.in(set::FinSetInt, elem) = in(elem, 1:set.n)

Base.show(io::IO, set::FinSetInt) = print(io, "FinSet($(set.n))")

""" Finite set given by Julia collection.

The underlying collection should be a Julia iterable of definite length. It may
be, but is not required to be, set-like (a subtype of `AbstractSet`).
"""
@auto_hash_equals struct FinSetCollection{S,T} <: FinSet{S,T}
  collection::S
end
FinSetCollection(collection::S) where S =
  FinSetCollection{S,eltype(collection)}(collection)

FinSet(collection::S) where {T, S<:Union{AbstractVector{T},AbstractSet{T}}} =
  FinSetCollection{S,T}(collection)

Base.iterate(set::FinSetCollection, args...) = iterate(set.collection, args...)
Base.length(set::FinSetCollection) = length(set.collection)
Base.in(set::FinSetCollection, elem) = in(elem, set.collection)

function Base.show(io::IO, set::FinSetCollection)
  print(io, "FinSet(")
  show(io, set.collection)
  print(io, ")")
end

""" Finite set whose elements are rows of a table.

The underlying table should be compliant with Tables.jl. For the sake of
uniformity, the rows are provided as named tuples, which assumes that the table
is not "extremely wide". This should not be a major limitation in practice but
see the Tables.jl documentation for further discussion.
"""
@auto_hash_equals struct TabularSet{Table,Row} <: FinSet{Table,Row}
  table::Table

  function TabularSet(table::Table) where Table
    schema = Tables.schema(table)
    new{Table,NamedTuple{schema.names,Tuple{schema.types...}}}(table)
  end
end

FinSet(nt::NamedTuple) = TabularSet(nt)

Base.iterate(set::TabularSet, args...) =
  iterate(Tables.namedtupleiterator(set.table), args...)
Base.length(set::TabularSet) = Tables.rowcount(set.table)
Base.collect(set::TabularSet) = Tables.rowtable(set.table)

function Base.show(io::IO, set::TabularSet)
  print(io, "TabularSet(")
  show(io, set.table)
  print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", set::TabularSet{T}) where T
  print(io, "$(length(set))-element TabularSet{$T}")
  if !get(io, :compact, false)
    println(io, ":")
    PrettyTables.pretty_table(io, set.table, nosubheader=true)
  end
end

function Base.show(io::IO, ::MIME"text/html", set::TabularSet)
  println(io, "
")
  println(io, "$(length(set))-element TabularSet")
  PrettyTables.pretty_table(io, set.table, backend=Val(:html), standalone=false,
                            nosubheader=true)
  println(io, "
")
end

# Discrete categories
#--------------------

""" Discrete category on a finite set.

The only morphisms in a discrete category are the identities, which are here
identified with the objects.
"""
@auto_hash_equals struct DiscreteCat{Ob,S<:FinSet{<:Any,Ob}} <: FinCat{Ob,Ob}
  set::S
end
DiscreteCat(n::Integer) = DiscreteCat(FinSet(n))

FinCat(s::Union{FinSet,Integer}) = DiscreteCat(s)

ob_generators(C::DiscreteCat) = C.set
hom_generators(::DiscreteCat) = ()
is_discrete(::DiscreteCat) = true
graph(C::DiscreteCat{Int,FinSetInt}) = Graph(length(C.set))

dom(C::DiscreteCat{T}, f) where T = f::T
codom(C::DiscreteCat{T}, f) where T = f::T
id(C::DiscreteCat{T}, x) where T = x::T
compose(C::DiscreteCat{T}, f, g) where T = (f::T == g::T) ? f :
  error("Nontrivial composite in discrete category: $f != $g")

hom_map(F::FinDomFunctor{<:DiscreteCat}, x) = id(codom(F), ob_map(F, x))

Base.show(io::IO, C::DiscreteCat{Int,FinSetInt}) =
  print(io, "FinCat($(length(C.set)))")

# Finite functions
##################

""" Function between finite sets.

The function can be defined implicitly by an arbitrary Julia function, in which
case it is evaluated lazily, or explictly by a vector of integers. In the vector
representation, the function (1↦1, 2↦3, 3↦2, 4↦3), for example, is represented
by the vector [1,3,2,3].

This type is mildly generalized by [`FinDomFunction`](@ref).
"""
const FinFunction{S, S′, Dom <: FinSet{S}, Codom <: FinSet{S′}} =
  SetFunction{Dom,Codom}

FinFunction(f::Function, dom, codom) =
  SetFunctionCallable(f, FinSet(dom), FinSet(codom))
FinFunction(::typeof(identity), args...) =
  IdentityFunction((FinSet(arg) for arg in args)...)
FinFunction(f::AbstractDict, args...) =
  FinFunctionDict(f, (FinSet(arg) for arg in args)...)

function FinFunction(f::AbstractVector{Int}, args...; index=false)
  if index == false
    FinDomFunctionVector(f, (FinSet(arg) for arg in args)...)
  else
    index = index == true ? nothing : index
    IndexedFinFunctionVector(f, args...; index=index)
  end
end
FinFunction(f::AbstractVector{Int}; kw...) =
  FinFunction(f, FinSet(isempty(f) ? 0 : maximum(f)); kw...)

Sets.show_type_constructor(io::IO, ::Type{<:FinFunction}) =
  print(io, "FinFunction")

""" Function out of a finite set.

This class of functions is convenient because it is exactly the class that can
be represented explicitly by a vector of values from the codomain.
"""
const FinDomFunction{S, Dom<:FinSet{S}, Codom<:SetOb} = SetFunction{Dom,Codom}

FinDomFunction(f::Function, dom, codom) =
  SetFunctionCallable(f, FinSet(dom), codom)
FinDomFunction(::typeof(identity), args...) =
  IdentityFunction((FinSet(arg) for arg in args)...)
FinDomFunction(f::AbstractDict, args...) = FinDomFunctionDict(f, args...)

function FinDomFunction(f::AbstractVector, args...; index=false)
  if index == false
    FinDomFunctionVector(f, args...)
  else
    index = index == true ? nothing : index
    IndexedFinDomFunctionVector(f, args...; index=index)
  end
end

Sets.show_type_constructor(io::IO, ::Type{<:FinDomFunction}) =
  print(io, "FinDomFunction")

# Note: Cartesian monoidal structure is implemented generically for Set but
# cocartesian only for FinSet.
@cocartesian_monoidal_instance FinSet FinFunction

Ob(C::FinCat{Int}) = FinSet(length(ob_generators(C)))
Ob(F::Functor{<:FinCat{Int}}) = FinDomFunction(collect_ob(F), Ob(codom(F)))

# Vector-based functions
#-----------------------

""" Function in **Set** represented by a vector.

The domain of this function is always of type `FinSet{Int}`, with elements of
the form ``{1,...,n}``.
"""
struct FinDomFunctionVector{T,V<:AbstractVector{T}, Codom<:SetOb{T}} <:
    FinDomFunction{Int,FinSetInt,Codom}
  func::V
  codom::Codom
end

FinDomFunctionVector(f::AbstractVector{T}) where T =
  FinDomFunctionVector(f, TypeSet{T}())

function FinDomFunctionVector(f::AbstractVector, dom::FinSet{Int}, codom)
  length(f) == length(dom) ||
    error("Length of vector $f does not match domain $dom")
  FinDomFunctionVector(f, codom)
end

dom(f::FinDomFunctionVector) = FinSet(length(f.func))

(f::FinDomFunctionVector)(x) = f.func[x]

function Base.show(io::IO, f::FinDomFunctionVector)
  print(io, "FinDomFunction($(f.func), ")
  Sets.show_domains(io, f)
  print(io, ")")
end

force(f::FinDomFunction{Int}) = FinDomFunctionVector(map(f, dom(f)), codom(f))
force(f::FinDomFunctionVector) = f

Base.collect(f::SetFunction) = force(f).func

""" Function in **FinSet** represented explicitly by a vector.
"""
const FinFunctionVector{S,T,V<:AbstractVector{T}} =
  FinDomFunctionVector{T,V,<:FinSet{S,T}}

Base.show(io::IO, f::FinFunctionVector) =
  print(io, "FinFunction($(f.func), $(length(dom(f))), $(length(codom(f))))")

Sets.do_compose(f::FinFunctionVector, g::FinDomFunctionVector) =
  FinDomFunctionVector(g.func[f.func], codom(g))

# Indexed vector-based functions
#-------------------------------

""" Indexed function out of a finite set of type `FinSet{Int}`.

Works in the same way as the special case of [`IndexedFinFunctionVector`](@ref),
except that the index is typically a dictionary, not a vector.
"""
struct IndexedFinDomFunctionVector{T,V<:AbstractVector{T},Index,Codom<:SetOb{T}} <:
    FinDomFunction{Int,FinSetInt,Codom}
  func::V
  index::Index
  codom::Codom
end

IndexedFinDomFunctionVector(f::AbstractVector{T}; kw...) where T =
  IndexedFinDomFunctionVector(f, TypeSet{T}(); kw...)

function IndexedFinDomFunctionVector(f::AbstractVector{T}, codom::SetOb{T};
                                     index=nothing) where T
  if isnothing(index)
    index = Dict{T,Vector{Int}}()
    for (i, x) in enumerate(f)
      push!(get!(index, x) do; Int[] end, i)
    end
  end
  IndexedFinDomFunctionVector(f, index, codom)
end

Base.:(==)(f::Union{FinDomFunctionVector,IndexedFinDomFunctionVector},
           g::Union{FinDomFunctionVector,IndexedFinDomFunctionVector}) =
  # Ignore index when comparing for equality.
  f.func == g.func && codom(f) == codom(g)

function Base.show(io::IO, f::IndexedFinDomFunctionVector)
  print(io, "FinDomFunction($(f.func), ")
  Sets.show_domains(io, f)
  print(io, ", index=true)")
end

dom(f::IndexedFinDomFunctionVector) = FinSet(length(f.func))
force(f::IndexedFinDomFunctionVector) = f

(f::IndexedFinDomFunctionVector)(x) = f.func[x]

""" Whether the given function is indexed, i.e., supports efficient preimages.
"""
is_indexed(f::SetFunction) = false
is_indexed(f::IdentityFunction) = true
is_indexed(f::IndexedFinDomFunctionVector) = true
is_indexed(f::FinDomFunctionVector{T,<:AbstractRange{T}}) where T = true

""" The preimage (inverse image) of the value y in the codomain.
"""
preimage(f::IdentityFunction, y) = SVector(y)
preimage(f::FinDomFunction, y) = [ x for x in dom(f) if f(x) == y ]
preimage(f::IndexedFinDomFunctionVector, y) = get_preimage_index(f.index, y)

@inline get_preimage_index(index::AbstractDict, y) = get(index, y, 1:0)
@inline get_preimage_index(index::AbstractVector, y) = index[y]

preimage(f::FinDomFunctionVector{T,<:AbstractRange{T}}, y::T) where T =
  # Both `in` and `searchsortedfirst` are specialized for AbstractRange.
  y ∈ f.func ? SVector(searchsortedfirst(f.func, y)) : SVector{0,Int}()

""" Indexed function between finite sets of type `FinSet{Int}`.

Indexed functions store both the forward map ``f: X → Y``, as a vector of
integers, and the backward map ``f: Y → X⁻¹``, as a vector of vectors of
integers, accessible through the [`preimage`](@ref) function. The backward map
is called the *index*. If it is not supplied through the keyword argument
`index`, it is computed when the object is constructed.

This type is mildly generalized by [`IndexedFinDomFunctionVector`](@ref).
"""
const IndexedFinFunctionVector{V,Index} =
  IndexedFinDomFunctionVector{Int,V,Index,FinSetInt}

function IndexedFinFunctionVector(f::AbstractVector{Int}; index=nothing)
  codom = isnothing(index) ? (isempty(f) ? 0 : maximum(f)) : length(index)
  IndexedFinFunctionVector(f, codom; index=index)
end

function IndexedFinFunctionVector(f::AbstractVector{Int}, codom; index=nothing)
  codom = FinSet(codom)
  if isnothing(index)
    index = [ Int[] for j in codom ]
    for (i, j) in enumerate(f)
      push!(index[j], i)
    end
  elseif length(index) != length(codom)
    error("Index length $(length(index)) does not match codomain $codom")
  end
  IndexedFinDomFunctionVector(f, index, codom)
end

Base.show(io::IO, f::IndexedFinFunctionVector) =
  print(io, "FinFunction($(f.func), $(length(dom(f))), $(length(codom(f))), index=true)")

# For now, we do not preserve or compose indices, only the function vectors.
Sets.do_compose(f::Union{FinFunctionVector,IndexedFinFunctionVector},
                g::Union{FinDomFunctionVector,IndexedFinDomFunctionVector}) =
  FinDomFunctionVector(g.func[f.func], codom(g))

# Dict-based functions
#---------------------

""" Function in **Set** represented by a dictionary.

The domain is a `FinSet{S}` where `S` is the type of the dictionary's `keys`
collection.
"""
@auto_hash_equals struct FinDomFunctionDict{K,D<:AbstractDict{K},Codom<:SetOb} <:
    FinDomFunction{D,FinSet{AbstractSet{K},K},Codom}
  func::D
  codom::Codom
end

FinDomFunctionDict(d::AbstractDict{K,V}) where {K,V} =
  FinDomFunctionDict(d, TypeSet{V}())

dom(f::FinDomFunctionDict) = FinSet(keys(f.func))

(f::FinDomFunctionDict)(x) = f.func[x]

function Base.show(io::IO, f::F) where F <: FinDomFunctionDict
  Sets.show_type_constructor(io, F)
  print(io, "(")
  show(io, f.func)
  print(io, ", ")
  Sets.show_domains(io, f, domain=false)
  print(io, ")")
end

force(f::FinDomFunction) =
  FinDomFunctionDict(Dict(x => f(x) for x in dom(f)), codom(f))
force(f::FinDomFunctionDict) = f

""" Function in **FinSet** represented by a dictionary.
"""
const FinFunctionDict{K,D<:AbstractDict{K},Codom<:FinSet} =
  FinDomFunctionDict{K,D,Codom}

FinFunctionDict(d::AbstractDict, codom::FinSet) = FinDomFunctionDict(d, codom)
FinFunctionDict(d::AbstractDict{K,V}) where {K,V} =
  FinDomFunctionDict(d, FinSet(Set(values(d))))

Sets.do_compose(f::FinFunctionDict{K,D}, g::FinDomFunctionDict) where {K,D} =
  FinDomFunctionDict(dicttype(D)(x => g.func[y] for (x,y) in pairs(f.func)),
                     codom(g))

# Limits
########

limit(Xs::EmptyDiagram{<:FinSet{Int}}) = Limit(Xs, SMultispan{0}(FinSet(1)))

universal(lim::Limit{<:FinSet{Int},<:EmptyDiagram}, cone::SMultispan{0}) =
  ConstantFunction(1, apex(cone), FinSet(1))

limit(Xs::SingletonDiagram{<:FinSet{Int}}) = limit(Xs, SpecializeLimit())

function limit(Xs::ObjectPair{<:FinSet{Int}})
  m, n = length.(Xs)
  indices = CartesianIndices((m, n))
  π1 = FinFunction(i -> indices[i][1], m*n, m)
  π2 = FinFunction(i -> indices[i][2], m*n, n)
  Limit(Xs, Span(π1, π2))
end

function universal(lim::Limit{<:FinSet{Int},<:ObjectPair}, cone::Span)
  f, g = cone
  m, n = length.(codom.(cone))
  indices = LinearIndices((m, n))
  FinFunction(i -> indices[f(i),g(i)], apex(cone), ob(lim))
end

function limit(Xs::DiscreteDiagram{<:FinSet{Int}})
  ns = length.(Xs)
  indices = CartesianIndices(Tuple(ns))
  n = prod(ns)
  πs = [FinFunction(i -> indices[i][j], n, ns[j]) for j in 1:length(ns)]
  Limit(Xs, Multispan(FinSet(n), πs))
end

function universal(lim::Limit{<:FinSet{Int},<:DiscreteDiagram}, cone::Multispan)
  ns = length.(codom.(cone))
  indices = LinearIndices(Tuple(ns))
  FinFunction(i -> indices[(f(i) for f in cone)...], apex(cone), ob(lim))
end

function limit(pair::ParallelPair{<:FinSet{Int}})
  f, g = pair
  m = length(dom(pair))
  eq = FinFunction(filter(i -> f(i) == g(i), 1:m), m)
  Limit(pair, SMultispan{1}(eq))
end

function limit(para::ParallelMorphisms{<:FinSet{Int}})
  @assert !isempty(para)
  f1, frest = para[1], para[2:end]
  m = length(dom(para))
  eq = FinFunction(filter(i -> all(f1(i) == f(i) for f in frest), 1:m), m)
  Limit(para, SMultispan{1}(eq))
end

function universal(lim::Limit{<:FinSet{Int},<:ParallelMorphisms},
                   cone::SMultispan{1})
  ι = collect(incl(lim))
  h = only(cone)
  FinFunction(Int[only(searchsorted(ι, h(i))) for i in dom(h)], length(ι))
end

""" Limit of finite sets with a reverse mapping or index into the limit set.

This type provides a fallback for limit algorithms that do not come with a
specialized algorithm to apply the universal property of the limit. In such
cases, you can explicitly construct a mapping from tuples of elements in the
feet of the limit cone to elements in the apex of the cone.

The index is constructed the first time it is needed. Thus there is no extra
cost to using this type if the universal property will not be invoked.
"""
mutable struct FinSetIndexedLimit{Ob<:FinSet,Diagram,Cone<:Multispan{Ob}} <:
    AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  index::Union{AbstractDict,Nothing}
end
FinSetIndexedLimit(diagram, cone::Multispan) =
  FinSetIndexedLimit(diagram, cone, nothing)

function make_limit_index(cone::Multispan{<:FinSet})
  πs = Tuple(legs(cone))
  index = Dict{Tuple{map(eltype∘codom, πs)...}, eltype(apex(cone))}()
  for x in apex(cone)
    index[map(π -> π(x), πs)] = x
  end
  return index
end

function universal(lim::FinSetIndexedLimit, cone::Multispan)
  if isnothing(lim.index)
    lim.index = make_limit_index(lim.cone)
  end
  fs = Tuple(legs(cone))
  FinFunction(Int[lim.index[map(f -> f(x), fs)] for x in apex(cone)],
              apex(cone), ob(lim))
end

""" Algorithm for limit of cospan or multicospan with feet being finite sets.

In the context of relational databases, such limits are called *joins*. The
trivial join algorithm is [`NestedLoopJoin`](@ref), which is algorithmically
equivalent to the generic algorithm `ComposeProductEqualizer`. The algorithms
[`HashJoin`](@ref) and [`SortMergeJoin`](@ref) are usually much faster. If you
are unsure what algorithm to pick, use [`SmartJoin`](@ref).
"""
abstract type JoinAlgorithm <: LimitAlgorithm end

""" Meta-algorithm for joins that attempts to pick an appropriate algorithm.
"""
struct SmartJoin <: JoinAlgorithm end

function limit(cospan::Multicospan{<:SetOb,<:FinDomFunction{Int}};
               alg::LimitAlgorithm=ComposeProductEqualizer())
  limit(cospan, alg)
end

function limit(cospan::Multicospan{<:SetOb,<:FinDomFunction{Int}}, ::SmartJoin)
  # Handle the important special case where one of the legs is a constant
  # (function out of a singleton set). In this case, we just need to take a
  # product of preimages of the constant value.
  funcs = legs(cospan)
  i = findfirst(f -> length(dom(f)) == 1, funcs)
  if !isnothing(i)
    c = funcs[i](1)
    ιs = map(deleteat(funcs, i)) do f
      FinFunction(preimage(f, c), dom(f))
    end
    x, πs = if length(ιs) == 1
      dom(only(ιs)), ιs
    else
      prod = product(map(dom, ιs))
      ob(prod), map(compose, legs(prod), ιs)
    end
    πs = insert(πs, i, ConstantFunction(1, x, FinSet(1)))
    return FinSetIndexedLimit(cospan, Multispan(πs))
  end

  # In the general case, for now we always just do a hash join, although
  # sort-merge joins can sometimes be faster.
  limit(cospan, HashJoin())
end

deleteat(vec::StaticVector, i) = StaticArrays.deleteat(vec, i)
deleteat(vec::Vector, i) = deleteat!(copy(vec), i)

insert(vec::StaticVector{N,T}, i, x::S) where {N,T,S} =
  StaticArrays.insert(similar_type(vec, typejoin(T,S))(vec), i, x)
insert(vec::Vector{T}, i, x::S) where {T,S} =
  insert!(collect(typejoin(T,S), vec), i, x)

""" [Nested-loop join](https://en.wikipedia.org/wiki/Nested_loop_join) algorithm.

This is the naive algorithm for computing joins.
"""
struct NestedLoopJoin <: JoinAlgorithm end

function limit(cospan::Multicospan{<:SetOb,<:FinDomFunction{Int}},
               ::NestedLoopJoin)
  # A nested-loop join is algorithmically the same as `ComposeProductEqualizer`,
  # but for completeness and performance we give a direct implementation here.
  funcs = legs(cospan)
  ns = map(length, feet(cospan))
  πs = map(_ -> Int[], funcs)
  for I in CartesianIndices(Tuple(ns))
    values = map((f, i) -> f(I[i]), funcs, eachindex(funcs))
    if all(==(values[1]), values)
      for i in eachindex(πs)
        push!(πs[i], I[i])
      end
    end
  end
  cone = Multispan(map((π,f) -> FinFunction(π, dom(f)), πs, funcs))
  FinSetIndexedLimit(cospan, cone)
end

""" [Sort-merge join](https://en.wikipedia.org/wiki/Sort-merge_join) algorithm.
"""
struct SortMergeJoin <: JoinAlgorithm end

function limit(cospan::Multicospan{<:SetOb,<:FinDomFunction{Int}},
               ::SortMergeJoin)
  funcs = map(collect, legs(cospan))
  sorts = map(sortperm, funcs)
  values = similar_mutable(funcs, eltype(apex(cospan)))
  ranges = similar_mutable(funcs, UnitRange{Int})

  function next_range!(i::Int)
    f, sort = funcs[i], sorts[i]
    n = length(f)
    start = last(ranges[i]) + 1
    ranges[i] = if start <= n
      val = values[i] = f[sort[start]]
      stop = start + 1
      while stop <= n && f[sort[stop]] == val; stop += 1 end
      start:(stop - 1)
    else
      start:n
    end
  end

  πs = map(_ -> Int[], funcs)
  for i in eachindex(ranges)
    ranges[i] = 0:0
    next_range!(i)
  end
  while !any(isempty, ranges)
    if all(==(values[1]), values)
      indices = CartesianIndices(Tuple(ranges))
      for i in eachindex(πs)
        append!(πs[i], (sorts[i][I[i]] for I in indices))
        next_range!(i)
      end
    else
      next_range!(argmin(values))
    end
  end
  cone = Multispan(map((π,f) -> FinFunction(π, length(f)), πs, funcs))
  FinSetIndexedLimit(cospan, cone)
end

similar_mutable(x::AbstractVector, T::Type) = similar(x, T)

function similar_mutable(x::StaticVector{N}, T::Type) where N
  # `similar` always returns an `MVector` but `setindex!(::MVector, args...)`
  # only works when the element type is a bits-type.
  isbitstype(T) ? similar(x, T) : SizedVector{N}(Vector{T}(undef, N))
end

""" [Hash join](https://en.wikipedia.org/wiki/Hash_join) algorithm.
"""
struct HashJoin <: JoinAlgorithm end

function limit(cospan::Multicospan{<:SetOb,<:FinDomFunction{Int}}, ::HashJoin)
  # We follow the standard terminology for hash joins: in a multiway hash join,
  # one function, called the *probe*, will be iterated over and need not be
  # indexed, whereas the other functions, call *build* inputs, must be indexed.
  #
  # We choose as probe the unindexed function with largest domain. If all
  # functions are already indexed, we arbitrarily choose the first one.
  i = argmax(map(legs(cospan)) do f
    is_indexed(f) ? -1 : length(dom(f))
  end)
  probe = legs(cospan)[i]
  builds = map(ensure_indexed, deleteat(legs(cospan), i))
  πs_build, π_probe = hash_join(builds, probe)
  FinSetIndexedLimit(cospan, Multispan(insert(πs_build, i, π_probe)))
end

function hash_join(builds::AbstractVector{<:FinDomFunction{Int}},
                   probe::FinDomFunction{Int})
  π_builds, πp = map(_ -> Int[], builds), Int[]
  for y in dom(probe)
    val = probe(y)
    preimages = map(build -> preimage(build, val), builds)
    n_preimages = Tuple(map(length, preimages))
    n = prod(n_preimages)
    if n > 0
      indices = CartesianIndices(n_preimages)
      for j in eachindex(π_builds)
        πb, xs = π_builds[j], preimages[j]
        append!(πb, (xs[I[j]] for I in indices))
      end
      append!(πp, (y for i in 1:n))
    end
  end
  (map(FinFunction, π_builds, map(dom, builds)), FinFunction(πp, dom(probe)))
end

function hash_join(builds::StaticVector{1,<:FinDomFunction{Int}},
                   probe::FinDomFunction{Int})
  πb, πp = hash_join(builds[1], probe)
  (SVector((πb,)), πp)
end
function hash_join(build::FinDomFunction{Int}, probe::FinDomFunction{Int})
  πb, πp = Int[], Int[]
  for y in dom(probe)
    xs = preimage(build, probe(y))
    n = length(xs)
    if n > 0
      append!(πb, xs)
      append!(πp, (y for i in 1:n))
    end
  end
  (FinFunction(πb, dom(build)), FinFunction(πp, dom(probe)))
end

ensure_indexed(f::FinFunction{Int,Int}) = is_indexed(f) ? f :
  FinFunction(collect(f), codom(f), index=true)
ensure_indexed(f::FinDomFunction{Int}) = is_indexed(f) ? f :
  FinDomFunction(collect(f), index=true)

function limit(d::BipartiteFreeDiagram{Ob,Hom}) where
    {Ob<:SetOb, Hom<:FinDomFunction{Int}}
  # As in a pullback, this method assumes that all objects in layer 2 have
  # incoming morphisms.
  @assert !any(isempty(incident(d, v, :tgt)) for v in vertices₂(d))
  d_original = d

  # For uniformity, e.g. when pairing below, ensure that all objects in layer 2
  # are type sets.
  if !all(x isa TypeSet for x in ob₂(d))
    d = map(d, ob₁=identity, ob₂=ensure_type_set, hom=ensure_type_set_codom)
  end

  # It is generally optimal to compute all equalizers (self joins) first, so as
  # to reduce the sizes of later pullbacks (joins) and products (cross joins).
  d, ιs = equalize_all(d)
  rem_vertices₂!(d, [v for v in vertices₂(d) if
                     length(incident(d, v, :tgt)) == 1])

  # Perform all pairings before computing any joins.
  d = pair_all(d)

  # Having done this preprocessing, if there are any nontrivial joins, perform
  # one of them and recurse; otherwise, we have at most a product to compute.
  #
  # In the binary case (`nv₁(d) == 2`), the preprocessing guarantees that there
  # is at most one nontrivial join, so there are no choices to make. When there
  # are multiple possible joins, do the one with smallest base cardinality
  # (product of sizes of relations to join). This is a simple greedy heuristic.
  # For more control over the order of the joins, create a UWD schedule.
  if nv₂(d) == 0
    # FIXME: Shouldn't need FinSetIndexedLimit in these special cases.
    if nv₁(d) == 1
      FinSetIndexedLimit(d_original, SMultispan{1}(ιs[1]))
    else
      πs = legs(product(SVector(ob₁(d)...)))
      FinSetIndexedLimit(d_original, Multispan(map(compose, πs, ιs)))
    end
  else
    # Select the join to perform.
    v = argmin(map(vertices₂(d)) do v
      edges = incident(d, v, :tgt)
      @assert length(edges) >= 2
      prod(e -> length(dom(hom(d, e))), edges)
    end)

    # Compute the pullback (inner join).
    join_edges = incident(d, v, :tgt)
    to_join = src(d, join_edges)
    to_keep = setdiff(vertices₁(d), to_join)
    pb = pullback(SVector(hom(d, join_edges)...), alg=SmartJoin())

    # Create a new bipartite diagram with joined vertices.
    d_joined = BipartiteFreeDiagram{Ob,Hom}()
    copy_parts!(d_joined, d, V₁=to_keep, V₂=setdiff(vertices₂(d),v), E=edges(d))
    joined = add_vertex₁!(d_joined, ob₁=apex(pb))
    for (u, π) in zip(to_join, legs(pb))
      for e in setdiff(incident(d, u, :src), join_edges)
        set_subparts!(d_joined, e, src=joined, hom=π⋅hom(d,e))
      end
    end
    rem_edges!(d_joined, join_edges)

    # Recursively compute the limit of the new diagram.
    lim = limit(d_joined)

    # Assemble limit cone from cones for pullback and reduced limit.
    πs = Vector{Hom}(undef, nv₁(d))
    for (i, u) in enumerate(to_join)
      πs[u] = compose(last(legs(lim)), legs(pb)[i], ιs[u])
    end
    for (i, u) in enumerate(to_keep)
      πs[u] = compose(legs(lim)[i], ιs[u])
    end
    FinSetIndexedLimit(d_original, Multispan(πs))
  end
end

ensure_type_set(s::FinSet) = TypeSet(eltype(s))
ensure_type_set(s::TypeSet) = s
ensure_type_set_codom(f::FinFunction) =
  SetFunctionCallable(f, dom(f), TypeSet(eltype(codom(f))))
ensure_type_set_codom(f::IndexedFinFunctionVector) =
  IndexedFinDomFunctionVector(f.func, index=f.index)
ensure_type_set_codom(f::FinDomFunction) = f

""" Compute all possible equalizers in a bipartite free diagram.

The result is a new bipartite free diagram that has the same vertices but is
*simple*, i.e., has no multiple edges. The list of inclusion morphisms into
layer 1 of the original diagram is also returned.
"""
function equalize_all(d::BipartiteFreeDiagram{Ob,Hom}) where {Ob,Hom}
  d_simple = BipartiteFreeDiagram{Ob,Hom}()
  copy_parts!(d_simple, d, V₂=vertices₂(d))
  ιs = map(vertices₁(d)) do u
    # Collect outgoing edges of u, key-ed by target vertex.
    out_edges = OrderedDict{Int,Vector{Int}}()
    for e in incident(d, u, :src)
      push!(get!(out_edges, tgt(d,e)) do; Int[] end, e)
    end

    # Equalize all sets of parallel edges out of u.
    ι = id(ob₁(d, u))
    for es in values(out_edges)
      if length(es) > 1
        fs = SVector((ι⋅f for f in hom(d, es))...)
        ι = incl(equalizer(fs)) ⋅ ι
      end
    end

    add_vertex₁!(d_simple, ob₁=dom(ι)) # == u
    for (v, es) in pairs(out_edges)
      add_edge!(d_simple, u, v, hom=ι⋅hom(d, first(es)))
    end
    ι
  end
  (d_simple, ιs)
end

""" Perform all possible pairings in a bipartite free diagram.

The resulting diagram has the same layer 1 vertices but a possibly reduced set
of layer 2 vertices. Layer 2 vertices are merged when they have exactly the same
multiset of adjacent vertices.
"""
function pair_all(d::BipartiteFreeDiagram{Ob,Hom}) where {Ob,Hom}
  d_paired = BipartiteFreeDiagram{Ob,Hom}()
  copy_parts!(d_paired, d, V₁=vertices₁(d))

  # Construct mapping to V₂ vertices from multisets of adjacent V₁ vertices.
  outmap = OrderedDict{Vector{Int},Vector{Int}}()
  for v in vertices₂(d)
    push!(get!(outmap, sort(inneighbors(d, v))) do; Int[] end, v)
  end

  for (srcs, tgts) in pairs(outmap)
    in_edges = map(tgts) do v
      sort(incident(d, v, :tgt), by=e->src(d,e))
    end
    if length(tgts) == 1
      v = add_vertex₂!(d_paired, ob₂=ob₂(d, only(tgts)))
      add_edges!(d_paired, srcs, fill(v, length(srcs)),
                 hom=hom(d, only(in_edges)))
    else
      prod = product(SVector(ob₂(d, tgts)...))
      v = add_vertex₂!(d_paired, ob₂=ob(prod))
      for (i,u) in enumerate(srcs)
        f = pair(prod, hom(d, getindex.(in_edges, i)))
        add_edge!(d_paired, u, v, hom=f)
      end
    end
  end
  d_paired
end

""" Limit of general diagram of FinSets computed by product-then-filter.

See `Limits.CompositePullback` for a very similar construction.
"""
struct FinSetCompositeLimit{Ob<:FinSet, Diagram,
                            Cone<:Multispan{Ob}, Prod<:Product{Ob},
                            Incl<:FinFunction} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  prod::Prod
  incl::Incl # Inclusion for the "multi-equalizer" in general formula.
end

limit(d::FreeDiagram{<:FinSet{Int}}) = limit(FinDomFunctor(d))

function limit(F::Functor{<:FinCat{Int},<:TypeCat{<:FinSet{Int}}})
  # Uses the general formula for limits in Set (Leinster, 2014, Basic Category
  # Theory, Example 5.1.22 / Equation 5.16). This method is simple and direct,
  # but extremely inefficient!
  J = dom(F)
  prod = product(map(x -> ob_map(F, x), ob_generators(J)))
  n, πs = length(ob(prod)), legs(prod)
  ι = FinFunction(filter(1:n) do i
    all(hom_generators(J)) do f
      s, t, h = dom(J, f), codom(J, f), hom_map(F, f)
      h(πs[s](i)) == πs[t](i)
    end
  end, n)
  cone = Multispan(dom(ι), map(x -> ι⋅πs[x], ob_generators(J)))
  FinSetCompositeLimit(F, cone, prod, ι)
end

function universal(lim::FinSetCompositeLimit, cone::Multispan{<:FinSet{Int}})
  ι = collect(lim.incl)
  h = universal(lim.prod, cone)
  FinFunction(Int[only(searchsorted(ι, h(i))) for i in dom(h)],
              apex(cone), ob(lim))
end

""" Limit of finite sets viewed as a table.

Any limit of finite sets can be canonically viewed as a table
([`TabularSet`](@ref)) whose columns are the legs of the limit cone and whose
rows correspond to elements of the limit object. To construct this table from an
already computed limit, call `TabularLimit(::AbstractLimit; ...)`. The column
names of the table are given by the optional argument `names`.

In this tabular form, applying the universal property of the limit is trivial
since it is just tupling. Thus, this representation can be useful when the
original limit algorithm does not support efficient application of the universal
property. On the other hand, this representation has the disadvantage of
generally making the element type of the limit set more complicated.
"""
const TabularLimit = Limit{<:TabularSet}

function TabularLimit(lim::AbstractLimit; names=nothing)
  πs = legs(lim)
  names = isnothing(names) ? (1:length(πs)) : names
  names = Tuple(column_name(name) for name in names)
  table = TabularSet(NamedTuple{names}(Tuple(map(collect, πs))))
  cone = Multispan(table, map(πs, eachindex(πs)) do π, i
    FinFunction(row -> Tables.getcolumn(row, i), table, codom(π))
  end)
  Limit(lim.diagram, cone)
end

function universal(lim::Limit{<:TabularSet{Table,Row}},
                   cone::Multispan) where {Table,Row}
  fs = Tuple(legs(cone))
  FinFunction(x -> Row(map(f -> f(x), fs)), apex(cone), ob(lim))
end

column_name(name) = Symbol(name)
column_name(i::Integer) = Symbol("x$i") # Same default as DataFrames.jl.

# Colimits
##########

# Colimits in Skel(FinSet)
#-------------------------

colimit(Xs::EmptyDiagram{<:FinSet{Int}}) = Colimit(Xs, SMulticospan{0}(FinSet(0)))

function universal(colim::Initial{<:FinSet{Int}}, cocone::SMulticospan{0})
  cod = apex(cocone)
  FinDomFunction(SVector{0,eltype(cod)}(), cod)
end

colimit(Xs::SingletonDiagram{<:FinSet{Int}}) = colimit(Xs, SpecializeColimit())

function colimit(Xs::ObjectPair{<:FinSet{Int}})
  m, n = length.(Xs)
  ι1 = FinFunction(1:m, m, m+n)
  ι2 = FinFunction(m+1:m+n, n, m+n)
  Colimit(Xs, Cospan(ι1, ι2))
end

function universal(colim::BinaryCoproduct{<:FinSet{Int}}, cocone::Cospan)
  f, g = cocone
  FinDomFunction(vcat(collect(f), collect(g)), ob(colim), apex(cocone))
end

function colimit(Xs::DiscreteDiagram{<:FinSet{Int}})
  ns = length.(Xs)
  n = sum(ns)
  offsets = [0,cumsum(ns)...]
  ιs = [FinFunction((1:ns[j]) .+ offsets[j],ns[j],n) for j in 1:length(ns)]
  Colimit(Xs, Multicospan(FinSet(n), ιs))
end

function universal(colim::Coproduct{<:FinSet{Int}}, cocone::Multicospan)
  cod = apex(cocone)
  FinDomFunction(mapreduce(collect, vcat, cocone, init=eltype(cod)[]),
                 ob(colim), cod)
end

function colimit(pair::ParallelPair{<:FinSet{Int}})
  f, g = pair
  m, n = length(dom(pair)), length(codom(pair))
  sets = IntDisjointSets(n)
  for i in 1:m
    union!(sets, f(i), g(i))
  end
  Colimit(pair, SMulticospan{1}(quotient_projection(sets)))
end

function colimit(para::ParallelMorphisms{<:FinSet{Int}})
  @assert !isempty(para)
  f1, frest = para[1], para[2:end]
  m, n = length(dom(para)), length(codom(para))
  sets = IntDisjointSets(n)
  for i in 1:m
    for f in frest
      union!(sets, f1(i), f(i))
    end
  end
  Colimit(para, SMulticospan{1}(quotient_projection(sets)))
end

function universal(coeq::Coequalizer{<:FinSet{Int}}, cocone::SMulticospan{1})
  pass_to_quotient(proj(coeq), only(cocone))
end

""" Create projection map π: X → X/∼ from partition of X.
"""
function quotient_projection(sets::IntDisjointSets)
  h = [ find_root!(sets, i) for i in 1:length(sets) ]
  roots = unique!(sort(h))
  FinFunction([ searchsortedfirst(roots, r) for r in h ], length(roots))
end

""" Given h: X → Y, pass to quotient q: X/~ → Y under projection π: X → X/~.
"""
function pass_to_quotient(π::FinFunction{Int,Int}, h::FinFunction{Int,Int})
  @assert dom(π) == dom(h)
  q = zeros(Int, length(codom(π)))
  for i in dom(h)
    j = π(i)
    if q[j] == 0
      q[j] = h(i)
    else
      q[j] == h(i) || error("Quotient map of colimit is ill-defined")
    end
  end
  any(==(0), q) && error("Projection map is not surjective")
  FinFunction(q, codom(h))
end

function pass_to_quotient(π::FinFunction{Int,Int}, h::FinDomFunction{Int})
  @assert dom(π) == dom(h)
  q = Vector{Union{Some{eltype(codom(h))},Nothing}}(nothing, length(codom(π)))
  for i in dom(h)
    j = π(i)
    if isnothing(q[j])
      q[j] = Some(h(i))
    else
      something(q[j]) == h(i) || error("Quotient map of colimit is ill-defined")
    end
  end
  any(isnothing, q) && error("Projection map is not surjective")
  FinDomFunction(map(something, q), codom(h))
end

function colimit(span::Multispan{<:FinSet{Int}})
  colimit(span, ComposeCoproductCoequalizer())
end

""" Colimit of general diagram of FinSets computed by coproduct-then-quotient.

See `Limits.CompositePushout` for a very similar construction.
"""
struct FinSetCompositeColimit{Ob<:FinSet, Diagram,
                              Cocone<:Multicospan{Ob}, Coprod<:Coproduct{Ob},
                              Proj<:FinFunction} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  coprod::Coprod
  proj::Proj # Projection for the "multi-coequalizer" in general formula.
end

function colimit(d::BipartiteFreeDiagram{<:FinSet{Int}})
  # As in a pushout, this method assume that all objects in layer 1 have
  # outgoing morphisms so that they can be excluded from the coproduct.
  @assert !any(isempty(incident(d, u, :src)) for u in vertices₁(d))
  coprod = coproduct(ob₂(d))
  n, ιs = length(ob(coprod)), legs(coprod)
  sets = IntDisjointSets(n)
  for u in vertices₁(d)
    out_edges = incident(d, u, :src)
    for (e1, e2) in zip(out_edges[1:end-1], out_edges[2:end])
      h1, h2 = hom(d, e1), hom(d, e2)
      ι1, ι2 = ιs[tgt(d, e1)], ιs[tgt(d, e2)]
      for i in ob₁(d, u)
        union!(sets, ι1(h1(i)), ι2(h2(i)))
      end
    end
  end
  π = quotient_projection(sets)
  cocone = Multicospan(codom(π), [ ιs[i]⋅π for i in vertices₂(d) ])
  FinSetCompositeColimit(d, cocone, coprod, π)
end

colimit(d::FreeDiagram{<:FinSet{Int}}) = colimit(FinDomFunctor(d))

function colimit(F::Functor{<:FinCat{Int},<:TypeCat{<:FinSet{Int}}})
  # Uses the general formula for colimits in Set (Leinster, 2014, Basic Category
  # Theory, Example 5.2.16).
  J = dom(F)
  coprod = coproduct(map(x -> ob_map(F, x), ob_generators(J)))
  n, ιs = length(ob(coprod)), legs(coprod)
  sets = IntDisjointSets(n)
  for f in hom_generators(J)
    s, t, h = dom(J, f), codom(J, f), hom_map(F, f)
    for i in dom(h)
      union!(sets, ιs[s](i), ιs[t](h(i)))
    end
  end
  π = quotient_projection(sets)
  cocone = Multicospan(codom(π), map(x -> ιs[x]⋅π, ob_generators(J)))
  FinSetCompositeColimit(F, cocone, coprod, π)
end

function universal(colim::FinSetCompositeColimit, cocone::Multicospan)
  h = universal(colim.coprod, cocone)
  pass_to_quotient(colim.proj, h)
end

# Colimits with names
#--------------------

""" Compute colimit of finite sets whose elements are meaningfully named.

This situation seems to be mathematically uninteresting but is practically
important. The colimit is computed by reduction to the skeleton of **FinSet**
(`FinSet{Int}`) and the names are assigned afterwards, following some reasonable
conventions and add tags where necessary to avoid name clashes.
"""
struct NamedColimit <: ColimitAlgorithm end

function colimit(::Type{<:Tuple{<:FinSet{<:Any,T},<:FinFunction}}, d) where
    {T <: Union{Symbol,AbstractString}}
  colimit(d, NamedColimit())
end

function colimit(d::FixedShapeFreeDiagram{<:FinSet{<:Any,T},Hom},
                 alg::NamedColimit) where {T,Hom}
  # Reducing to the case of bipartite free diagrams is a bit lazy, but at least
  # using `SpecializeColimit` below should avoid some gross inefficiencies.
  colimit(BipartiteFreeDiagram{FinSet{<:Any,T},Hom}(d), alg)
end
function colimit(d::BipartiteFreeDiagram{<:FinSet{<:Any,T}}, ::NamedColimit) where T
  # Compute colimit of diagram in the skeleton of FinSet (`FinSet{Int}`).
  # Note: no performance would be gained by using `DisjointSets{T}` from
  # DataStructures.jl because it is just a wrapper around `IntDisjointSets` that
  # internally builds the very same indices that we use below.
  sets₁_skel = map(set -> skeletize(set, index=false), ob₁(d))
  sets₂_skel = map(set -> skeletize(set, index=true), ob₂(d))
  funcs = map(edges(d)) do e
    skeletize(hom(d,e), sets₁_skel[src(d,e)], sets₂_skel[tgt(d,e)])
  end
  d_skel = BipartiteFreeDiagram{FinSetInt,eltype(funcs)}()
  add_vertices₁!(d_skel, nv₁(d), ob₁=dom.(sets₁_skel))
  add_vertices₂!(d_skel, nv₂(d), ob₂=dom.(sets₂_skel))
  add_edges!(d_skel, src(d), tgt(d), hom=funcs)
  colim_skel = colimit(d_skel, SpecializeColimit())

  # Assign elements/names to the colimit set.
  elems = Vector{T}(undef, length(apex(colim_skel)))
  for (ι, Y) in zip(colim_skel, sets₂_skel)
    for i in dom(Y)
      elems[ι(i)] = Y(i)
    end
  end
  # The vector should already be filled, but to reduce arbitrariness we prefer
  # names from the layer 1 sets whenever possible. For example, when computing a
  # pushout, we prefer names from the apex of cospan to names from the feet.
  for (u, X) in zip(vertices₁(d_skel), sets₁_skel)
    e = first(incident(d_skel, u, :src))
    f, ι = hom(d_skel, e), legs(colim_skel)[tgt(d_skel, e)]
    for i in dom(X)
      elems[ι(f(i))] = X(i)
    end
  end
  # Eliminate clashes in provisional list of names.
  unique_by_tagging!(elems)

  ιs = map(colim_skel, sets₂_skel) do ι, Y
    FinFunction(Dict(Y(i) => elems[ι(i)] for i in dom(Y)), FinSet(elems))
  end
  Colimit(d, Multicospan(FinSet(elems), ιs))
end

function skeletize(set::FinSet; index::Bool=false)
  # FIXME: We should support `unique_index` and it should be used here.
  FinDomFunction(collect(set), set, index=index)
end
function skeletize(f::FinFunction, X, Y)
  FinFunction(i -> only(preimage(Y, f(X(i)))), dom(X), dom(Y))
end

""" Make list of elements unique by adding tags if necessary.

If the elements are already unique, they will not be mutated.
"""
function unique_by_tagging!(elems::AbstractVector{T}; tag=default_tag) where T
  tag_count = Dict{T,Int}()
  for x in elems
    tag_count[x] = haskey(tag_count, x) ? 1 : 0
  end
  for (i, x) in enumerate(elems)
    (j = tag_count[x]) > 0 || continue
    tagged = tag(x, j)
    @assert !haskey(tag_count, tagged) # Don't conflict with original elems!
    elems[i] = tagged
    tag_count[x] += 1
  end
  elems
end

default_tag(x::Symbol, t) = Symbol(x, "#", t)
default_tag(x::AbstractString, t) = string(x, "#", t)

# Pushout complements
#--------------------

""" Compute a pushout complement of finite sets, if possible.

Given functions ``l: I → L`` and ``m: L → G`` to form a pushout square

    l
  L ← I
m ↓   ↓k
  G ← K
    g

define the set ``K := G / m(L / l(I))`` and take ``g: K ↪ G`` to be the
inclusion. Then the map ``k: I → K`` is determined by the map ``l⋅m: I → G``
from the requirement that the square commutes.

Pushout complements exist only if the identification condition is satisfied. An
error will be raised if the pushout complement cannot be constructed. To check
this in advance, use [`can_pushout_complement`](@ref).
"""
function pushout_complement(pair::ComposablePair{<:FinSet{Int}})
  l, m = pair
  I, L, G = dom(l), codom(l), codom(m)

  # Construct inclusion g: K ↪ G.
  l_image = Set(collect(l))
  m_image = Set([ m(x) for x in L if x ∉ l_image ])
  g = FinFunction([x for x in G if x ∉ m_image], G)
  K = dom(g)

  # Construct morphism k: I → K using partial inverse of g.
  g_inv = Dict{Int,Int}(zip(collect(g), K))
  k = FinFunction(map(I) do x
    y = m(l(x))
    get(g_inv, y) do; error("Identification failed for domain element $x") end
  end, I, K)

  return ComposablePair(k, g)
end

can_pushout_complement(pair::ComposablePair{<:FinSet{Int}}) =
  all(isempty, id_condition(pair))

""" Check identification condition for pushout complement of finite sets.

The identification condition says that the functions do not map (1) both a
deleted item and a preserved item in L to the same item in G or (2) two distinct
deleted items to the same item. It is trivially satisfied for injective functions.

Returns pair of iterators of

  (1) a nondeleted item that maps to a deleted item in G
  (2) a pair of distinct items in L that are deleted yet mapped to the same
      item in G.
"""
function id_condition(pair::ComposablePair{<:FinSet{Int}})
  l, m = pair
  l_image = Set(collect(l))
  l_imageᶜ = [ x for x in codom(l) if x ∉ l_image ]
  m_image = Set(map(m, l_imageᶜ))
  ((i for i in l_image if m(i) ∈ m_image),
   ((i, j) for i in eachindex(l_imageᶜ)
           for j in i+1:length(l_imageᶜ)
           if m(l_imageᶜ[i]) == m(l_imageᶜ[j])))
end

# Subsets
#########

""" Subset of a finite set.
"""
const SubFinSet{S,T} = Subobject{<:FinSet{S,T}}

Subobject(X::FinSet, f) = Subobject(FinFunction(f, X))
SubFinSet(X, f) = Subobject(FinFunction(f, X))

force(A::SubFinSet{Int}) = Subobject(force(hom(A)))
Base.collect(A::SubFinSet) = collect(hom(A))
Base.sort(A::SubFinSet) = SubFinSet(ob(A), sort(collect(A)))

const AbstractBoolVector = Union{AbstractVector{Bool},BitVector}

""" Subset of a finite set represented as a boolean vector.

This is the subobject classifier representation since `Bool` is the subobject
classifier for `Set`.
"""
@auto_hash_equals struct SubFinSetVector{S<:FinSet} <: Subobject{S}
  set::S
  predicate::AbstractBoolVector

  function SubFinSetVector(X::S, pred::AbstractBoolVector) where S<:FinSet
    length(pred) == length(X) ||
      error("Size of predicate $pred does not equal size of object $X")
    new{S}(X, pred)
  end
end

Subobject(X::FinSet, pred::AbstractBoolVector) = SubFinSetVector(X, pred)
SubFinSet(pred::AbstractBoolVector) = Subobject(FinSet(length(pred)), pred)

ob(A::SubFinSetVector) = A.set
hom(A::SubFinSetVector) = FinFunction(findall(A.predicate), A.set)
predicate(A::SubFinSetVector) = A.predicate

function predicate(A::SubFinSet)
  f = hom(A)
  pred = falses(length(codom(f)))
  for x in dom(f)
    pred[f(x)] = true
  end
  pred
end

@instance SubobjectLattice{FinSet,SubFinSet} begin
  @import ob
  meet(A::SubFinSet, B::SubFinSet) = meet(A, B, SubOpBoolean())
  join(A::SubFinSet, B::SubFinSet) = join(A, B, SubOpBoolean())
  top(X::FinSet) = top(X, SubOpWithLimits())
  bottom(X::FinSet) = bottom(X, SubOpWithLimits())
end

""" Algorithm to compute subobject operations using elementwise boolean logic.
"""
struct SubOpBoolean <: SubOpAlgorithm end

meet(A::SubFinSet{Int}, B::SubFinSet{Int}, ::SubOpBoolean) =
  SubFinSet(predicate(A) .& predicate(B))
join(A::SubFinSet{Int}, B::SubFinSet{Int}, ::SubOpBoolean) =
  SubFinSet(predicate(A) .| predicate(B))
top(X::FinSet{Int}, ::SubOpBoolean) = SubFinSet(trues(length(X)))
bottom(X::FinSet{Int}, ::SubOpBoolean) = SubFinSet(falses(length(X)))

end


text:
module Move_FourBar2

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)
massProperties = "Aluminium"
vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly Bar(;Lx = 0.1, Ly=Lx/5, Lz=Ly) begin
   frame0 = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vmat1))
   frame1 = Modia3D.Object3D(frame0; r=[-Lx/2, 0.0, 0.0])
   frame2 = Modia3D.Object3D(frame0; r=[ Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2)
   cyl1   = Modia3D.Object3D(frame1, cyl)
   cyl2   = Modia3D.Object3D(frame2, cyl)
end

@assembly Fourbar2(;Lx = 0.1, Ly=Lx/5, Lz=Ly, groundWidth=Lx, groundHeight=0.1*Lx) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))
   frame1 = Modia3D.Object3D(world ; r=[Lx/2, 0.0, groundWidth/2])
   frame2 = Modia3D.Object3D(frame1; r=[Lx/2, 0.0, 0.0])
   ground = Modia3D.Object3D(world , Modia3D.Box(3*Lx, groundWidth, groundHeight; material=groundMaterial);
                                              r = [1.5*Lx, -groundHeight/2, groundWidth/2],
                                              R = ModiaMath.rot_nxy([-1,0,0], [0,0,1]))
   bar1   = Bar(Lx=Lx, Ly=Ly)
   bar2   = Bar(Lx=Lx, Ly=Ly)
   rev1   = Modia3D.Revolute(frame1     , bar1.frame1; phi_start =  pi/2)
   rev2   = Modia3D.Revolute(bar1.frame2, bar2.frame1; phi_start = -pi/2)

   Modia3D.updatePosition!(world)    # need to be improved (made automatic)
   L3         = Modia3D.distance(frame2,bar2.frame2)
   phi3_start = Modia3D.planarRotationAngle(frame2, bar2.frame2)
   bar3       = Bar(Lx=L3, Ly=Ly)

   rev3   = Modia3D.Revolute(frame2     , bar3.frame1; phi_start = phi3_start)
   rev4   = Modia3D.Revolute(bar3.frame2, bar2.frame2; phi_start = NaN)
end

@signal Sine(;yStart = 0.0, A = 1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Sine, sim::ModiaMath.SimulationState)
    signal.y.value = signal.yStart + signal.A*sin(2*pi*signal.freqHz*sim.time)
end

Lx = 1.0
@assembly Move2 begin
   fourbar = Fourbar2(Lx=Lx)

   sine = Sine(yStart=fourbar.rev1.phi.start, A=-pi/4, freqHz=0.5)
   sig  = Modia3D.SignalToFlangeAngle(sine.y)
   Modia3D.connect(sig, fourbar.rev1)
end

model  = Modia3D.SimulationModel( Move2(), analysis=ModiaMath.KinematicAnalysis )
result = ModiaMath.simulate!(model, stopTime=3.0, log=true)

println("... success of Move_FourBar2.jl!")
end

ext: jl


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# -----------------------------------------------------------------------------------------
# Constants for the tests.
const TOTAL_SIZE  = 1024
const END_SIZE    = 256
const STEP_SIZE   = 32

# -----------------------------------------------------------------------------------------
# Define the iterator to iterate over all the hash chunks.
struct StressChunks
    max::Int
end

function first(
            iter::StressChunks
        )::Vector{Int}
    return [ iter.max ]
end

function next(
            iter::StressChunks,
            prev::Union{Vector{Int}, Nothing} = nothing
        )::Bool

    while length(prev) > 0 && prev[end] <= STEP_SIZE
        pop!(prev)
    end

    if length(prev) == 0
        return false
    end

    prev[end] -= STEP_SIZE
    push!(prev, iter.max - sum(prev))

    return true
end

function increased_size(
            iter::StressChunks
        )::Union{StressChunks, Nothing}
    if iter.max == END_SIZE
        return nothing
    end
    return StressChunks(iter.max + STEP_SIZE)
end

mutable struct Stress
    front::Union{StressChunks, Nothing}
    back::Union{StressChunks, Nothing}
    Stress() = new(nothing, nothing)
end

struct StressState
    front::Union{Vector{Int}, Nothing}
    back::Union{Vector{Int}, Nothing}
end

function Base.iterate(
            iter::Stress,
            state::Union{StressState, Nothing} = nothing
        )::Union{Tuple{Vector{Int}, StressState}, Nothing}

    function result(st::StressState)
        if st.front === nothing
            return vcat(TOTAL_SIZE - sum(st.back), st.back)
        elseif st.back === nothing
            return vcat(st.front, TOTAL_SIZE - sum(st.front))
        end
        return vcat(st.front, TOTAL_SIZE - sum(st.front) - sum(st.back), st.back)
    end

    # -------------------------------------------------------------------------------------
    # Not state, let's get this show on the road.
    if state === nothing
        iter.back = StressChunks(STEP_SIZE)
        state     = StressState(nothing, first(iter.back))
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # If the back iterator is nothing, create the smallest one.
    if iter.back === nothing
        iter.back = StressChunks(STEP_SIZE)
        state     = StressState(state.front, first(iter.back))
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # Increment the back iterator.
    if next(iter.back, state.back) == true
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # Back iterator completed, increase it's size and grab the first item.
    iter.back = increased_size(iter.back)

    if iter.back !== nothing
        state = StressState(state.front, first(iter.back))
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # Back iterator has finished with END_SIZE, clear it out before incrementing the front
    # iterator.
    iter.back = nothing

    # -------------------------------------------------------------------------------------
    # Back iterator has finished with END_SIZE, create the front iterator if needed.
    if iter.front === nothing
        iter.front = StressChunks(STEP_SIZE)
        state      = StressState(first(iter.front), nothing)
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # Back iterator has finished with END_SIZE, increment the front iterator.
    if next(iter.front, state.front) == true
        state = StressState(state.front, nothing)
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # Front iterator completed, increase it's size and grab the first item.
    iter.front = increased_size(iter.front)

    if iter.front !== nothing
        state = StressState(first(iter.front), nothing)
        return (result(state), state)
    end

    # -------------------------------------------------------------------------------------
    # Front iterator has finsihed with END_SIZE, all done.
    return nothing
end

struct StressIter
    data::Vector{UInt8}
    stress::Stress
    StressIter(data) = new(data, Stress())
end

function Base.iterate(iter::StressIter, state=nothing)
    next = iterate(iter.stress, state)

    if next === nothing
        return nothing
    end

    data = Vector()

    pos = 1
    for size in next[1]
        push!(data, view(iter.data, pos:pos+size-1))
        pos += size
    end

    return (data, next[2])
end


text:
#
# Sets inicial velocities for particles in p.
#

function initial_velocity(T, m, N)
    R = 8.314_462_618_153_24
    vmod = sqrt.( -2R*T * log.(-rand(N) .+ 1)./m )
    angle = 2*π*rand(N)
    return [ [vmod[i] * cos(angle[i]), vmod[i] * sin(angle[i])] for i=1:N ]
end


text:
## Preloads
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra
using Measures

#internal AREPO units in cgs
arepoLength = 3.0856e20
arepoMass = 1.991e33
arepoVel = 1.0e5

#AREPO conversion factors
arepoTime = arepoLength/arepoVel
arepoDensity = arepoMass/arepoLength/arepoLength/arepoLength
arepoEnergy= arepoMass*arepoVel*arepoVel
arepoColumnDensity = arepoMass/arepoLength/arepoLength

#Helium abundance
ABHE = 0.1

#Mass of proton
mp=1.67e-24

#Grid size
nx,ny,nz=500,500,500

#Trial Read in
y =  Array{Float32}(undef, 125000000)

test = read("/Users/aksaydjari/Dropbox/GradSchool_AKS/Doug/ExtData/CloudFactory/density_grid_x1480y1390z1395",y)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
include("../src/util.jl")
include("../src/LUMND.jl")

using .LUMND
using MatrixMarket

function write_array(io, v)
    if length(v) == 0
        return
    end
    @printf(io, "%d", v[1])
    for v_ in v[2:end]
        @printf(io, ",%d", v_)
    end
end

# Inputs:
#   - A, a sparse matrix in CSC format of size N x N
#   - lvl, the number of level in the ND ordering
#   - prefix, something describing the problem. The 2 files (see below) will be saved to
#       prefix_ord_lvl.txt
#       prefix_clust_lvl.txt
#     where lvl are the # of level (lvl)
# Outputs: none. Write the ordering & clustering to file
#
# Separators are numbered from leaves to top, i.e., for lvl=3, we have 2^3-1=7 separators, numbered 0 to 6
#          6
#        4   5
#       0 1 2 3
#
# Description of the format:
#   prefix_ord_lvl.txt
#       Line 1: lvl #sep
#           where lvl is lvl, and #sep the total number of separators (including leaves), i.e., 2^(lvl)-1
#       Line 2...#sep+1
#           line i corresponds to separator i-2
#           Each line is
#               sepid;i0,i1,i2,...,in
#           where sepid is the separator id [0...2^lvl-2]
#           and i0, i1, ... are the degrees of freedom in the original matrix (i.e., rows or columns of A), from 0 to N-1
#
#   prefix_clust_lvl.txt
#       Line 1: lvl #sep
#           where lvl is lvl, and #sep the total number of separators (including leaves), i.e., 2^(lvl)-1
#       Line 2...#sep+1
#           line i corresponds to separator i-2
#           A separator at lvl l as a hierarchy with l+1 levels
#           The hierarchy can be seen as a sequence of nested interval, where the top (root) one is 'everything'
#           Each level as at least one interval, but potentially empty
#           Each line contains
#               sepid;interval 0;interval 1;interval 2;...;interval k
#           Interval 0 indicates how we 'group' the dofs (described in prefix_ord_lvl)
#           For instance, if the dofs are
#                       4,7,2,6,11,15,17,25,41,12
#           And interval 0 is
#                       0,3,4,8,10
#           This means we have 4 clusters at the leaf-level:
#                       4,7,2 - 6 - 11,15,17,25 - 41,12
#           Then if inteval 1 is
#                       0,2,4
#           It means at the next level, the clusters become:
#                       4,7,2,6 - 11,15,17,25,41,12
#           Then finally if interval 2 is
#                       0,2
#           It means at the next level, the clusters become:
#                       4,7,2,6,11,15,17,25,41,12
#
# Everything is stored in base-0 by convention

function order_and_write_file(A, lvl, prefix)

    # Compute ordering
    (subseps, hrch, dofs) = LUMND.mnd(A+A', lvl) # Works on the sparsity pattern of A+A'

    # Sort dofs
    for l = 1:lvl
        for s = 1:length(dofs[l])
            d_ = dofs[l][s]
            s_ = subseps[d_]
            p_ = sortperm(s_)
            dofs[l][s] = dofs[l][s][p_]
        end
    end

    # Write
    ordering_name = @sprintf("%s_ord_%d.txt", prefix, lvl)
    clustering_name = @sprintf("%s_clust_%d.txt", prefix, lvl)
    @printf("Data saved to %s and %s\n", ordering_name, clustering_name)

    # Write ordering
    file = open(ordering_name, "w")
    @printf(file, "%d %d\n", lvl, 2^lvl-1) # Lvl, # seps
    k = 0 # Base-0
    for l = 1:lvl
        for s = 1:length(dofs[l])
            d_ = dofs[l][s]
            print(file, @sprintf("%d;", k))
            write_array(file, dofs[l][s] .- 1) # Base-0
            println(file)
            k += 1
        end
    end
    close(file)

    # Write partitionning
    file = open(clustering_name, "w")
    @printf(file, "%d %d\n", lvl, 2^lvl-1) # Lvl, # seps
    k = 0 # Base-0
    for l = 1:lvl
        for s = 1:length(dofs[l])
            d_ = dofs[l][s]
            s_ = subseps[d_]
            h_ = hrch[l][s]
            # The leaf-level pointers
            old_ptr = sorted_to_ptr(s_)
            print(file, @sprintf("%d;", k));
            write_array(file, old_ptr .- 1) # Base-0
            print(file, ";")
            # The hierarchy of intervals
            for d = length(h_.nit)-1:-1:1
                ptr = [n.c[1].i for n in h_.nit[d]]
                push!(ptr, length(h_.nit[d+1]) + 1)
                write_array(file, ptr .- 1) # Base-0
                print(file,";")
                old_ptr = ptr
            end
            println(file)
            k += 1
        end
    end
    close(file)

end

n = 20
d = 2
A = Ad(n,d) # n^d laplacian
lvl = 5
prefix = @sprintf("lapl_%d_%d", n, d)
order_and_write_file(A, lvl, prefix)

# We can also finally save the matrix
matfile = @sprintf("%s.mtx", prefix)
mmwrite(matfile, A)
@printf("Matrix saved in MatrixMarket format to %s\n", matfile)


text:
const TEST_NAME = "Philox2x64"

include("common.jl")

using RandomNumbers.Random123

r = Philox2x(UInt64, 123)

test_all(r, 100_000_000)


text:
using Inflector
using Test

@testset "Inflector.jl" begin
    # Write your tests here.
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2018, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    Crimisini(tilesize)

Examplar-based inpainting based on confidence
and isophote maps. `tilesize` is the patch size
as a tuple of integers.

## References

Crimisini, A., Pérez, P., Toyama, K., 2004. Region Filling
and Object Removal by Examplar-based Image Inpainting.
"""
struct Crimisini{N} <: InpaintAlgo
  tilesize::Dims{N} # patch size
end

Crimisini(tilesize::Vararg{Int,N}) where {N} = Crimisini{N}(tilesize)

# implementation follows the notation in the paper
function inpaint_impl(img::AbstractArray{T,N}, mask::BitArray{N}, algo::Crimisini{N}) where {T,N}
  # use all CPU cores in FFT
  set_num_threads(cpucores())

  # patch (or tile) size
  tilesize = algo.tilesize

  # already filled region
  ϕ = .!mask

  # initialize confidence map
  C = Float64.(ϕ)

  # pad arrays
  prepad  = [(tilesize[i]-1) ÷ 2 for i=1:N]
  postpad = [(tilesize[i])   ÷ 2 for i=1:N]
  padimg = parent(padarray(img, Pad(:symmetric, prepad, postpad)))
  ϕ = parent(padarray(ϕ, Fill(true, prepad, postpad)))
  C = parent(padarray(C, Fill(0.0,  prepad, postpad)))

  # fix any invalid pixel value (e.g. NaN) inside of the mask
  padimg[isnan.(padimg)] .= zero(T)

  # inpainting frontier
  δΩ = findall(dilate(ϕ) .& .!ϕ)

  while !isempty(δΩ)
    # update confidence values in frontier
    for p in δΩ
      c, b = selectpatch((C, ϕ), tilesize, p)
      C[p] = sum(c[b]) / prod(tilesize)
    end

    # isophote map
    grads = pointgradients(padimg, δΩ)
    direc = pointgradients(ϕ, δΩ)
    D = vec(abs.(sum(grads.*direc, dims=2)))
    D /= maximum(D)

    # select patch in frontier
    idx = argmax(C[δΩ].*D)
    p = δΩ[idx]
    ψₚ, bₚ = selectpatch((padimg, ϕ), tilesize, p)

    # compute distance to all other patches
    Δ = convdist(padimg, ψₚ, weights=bₚ)

    # only consider patches in filled region
    Δ[mask] .= Inf

    # find index in padded arrays
    sub = argmin(Δ)
    q = ntuple(i -> sub[i] + (tilesize[i]-1)÷2, N)

    # select best candidate
    ψᵦ, bᵦ = selectpatch((padimg, ϕ), tilesize, q)

    # paste patch and mark pixels as painted
    b = bᵦ .& .!bₚ
    ψₚ[b] .= ψᵦ[b]
    bₚ[b] .= true

    # update frontier
    δΩ = findall(dilate(ϕ) .& .!ϕ)
  end

  view(padimg, [1+prepad[i]:size(padimg,i)-postpad[i] for i=1:N]...)
end


text:
@system C4Base(CBase) begin
    Cm(Ci): mesophyll_co2 ~ track(u"μbar")

    gbs: bundle_sheath_conductance => 0.003 ~ preserve(u"mol/m^2/s/bar" #= CO2 =#, parameter) # bundle sheath conductance to CO2, mol m-2 s-1
    # gi => 1.0 ~ preserve(parameter) # conductance to CO2 from intercelluar to mesophyle, mol m-2 s-1, assumed
end

@system C4c(C4Base) begin
    # Kp25: Michaelis constant for PEP caboxylase for CO2
    Kp25: pep_carboxylase_constant_for_co2_at_25 => 80 ~ preserve(u"μbar", parameter)
    Kp(Kp25, kTQ10): pep_carboxylase_constant_for_co2 => begin
        Kp25 * kTQ10
    end ~ track(u"μbar")

    Vpm25: maximum_pep_carboxylation_rate_for_co2_at_25 => 70 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    EaVp: activation_energy_for_pep_carboxylation => 75.1 ~ preserve(u"kJ/mol", parameter)
    Vpmax(Vpm25, kT, EaVp, kN): maximum_pep_carboxylation_rate => begin
        Vpm25 * kT(EaVp) * kN
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    # PEP regeneration limited Vp, value adopted from vC book
    Vpr25: pep_regeneration_rate_for_co2_at_25 => 80 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    Vpr(Vpr25, kTQ10): pep_regeneration_rate => begin
        Vpr25 * kTQ10
    end ~ track(u"μmol/m^2/s" #= CO2 =#)
    Vp(Vpmax, Vpr, Cm, Kp): pep_carboxylation_rate => begin
        # PEP carboxylation rate, that is the rate of C4 acid generation
        (Cm * Vpmax) / (Cm + Kp)
    end ~ track(u"μmol/m^2/s" #= CO2 =#, max=Vpr)

    Vcm25: maximum_carboxylation_rate_at_25 => 50 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    # EaVc: Sage (2002) JXB
    EaVc: activation_energy_for_carboxylation => 55.9 ~ preserve(u"kJ/mol", parameter)
    Vcmax(Vcm25, kT, EaVc, kN): maximum_carboxylation_rate => begin
        Vcm25 * kT(EaVc) * kN
    end ~ track(u"μmol/m^2/s" #= CO2 =#)
end

@system C4j(C4Base) begin
    Jm25: maximum_electron_transport_rate_at_25 => 300 ~ preserve(u"μmol/m^2/s" #= Electron =#, parameter)
    Haj: activation_enthalpy_for_electron_transport => 32.8 ~ preserve(u"kJ/mol", parameter)
    Hdj: deactivation_enthalpy_for_electron_transport => 220 ~ preserve(u"kJ/mol", parameter)
    Toj: optimum_temperature_for_electron_transport => 33.6 ~ preserve(u"°C", parameter)
    Jmax(Jm25, kTpeak, Haj, Hdj, Toj, kN): maximum_electron_transport_rate => begin
        Jm25 * kTpeak(Haj, Hdj, Toj) * kN
    end ~ track(u"μmol/m^2/s" #= Electron =#)

    # θ: sharpness of transition from light limitation to light saturation
    θ: light_transition_sharpness => 0.5 ~ preserve(parameter)
    J(I2, Jmax, θ): electron_transport_rate => begin
        a = θ
        b = -(I2+Jmax)
        c = I2*Jmax
        a*J^2 + b*J + c ⩵ 0
    end ~ solve(lower=0, upper=Jmax, u"μmol/m^2/s")
end

@system C4r(C4Base) begin
    # Kc25: Michaelis constant of rubisco for CO2 of C4 plants (2.5 times that of tobacco), ubar, Von Caemmerer 2000
    Kc25: rubisco_constant_for_co2_at_25 => 650 ~ preserve(u"μbar", parameter)
    Eac: activation_energy_for_co2 => 59.4 ~ preserve(u"kJ/mol", parameter)
    Kc(kT, Kc25, Eac): rubisco_constant_for_co2 => begin
        Kc25 * kT(Eac)
    end ~ track(u"μbar")

    # Ko25: Michaelis constant of rubisco for O2 (2.5 times C3), mbar
    Ko25: rubisco_constant_for_o2_at_25 => 450 ~ preserve(u"mbar", parameter)
    # Activation energy for Ko, Bernacchi (2001)
    Eao: activation_energy_for_o2 => 36 ~ preserve(u"kJ/mol", parameter)
    Ko(Ko25, kT, Eao): rubisco_constant_for_o2 => begin
        Ko25 * kT(Eao)
    end ~ track(u"mbar")

    # mesophyll O2 partial pressure
    Om: mesophyll_o2_partial_pressure => 210 ~ preserve(u"mbar", parameter)
    Km(Kc, Om, Ko): rubisco_constant_for_co2_with_o2 => begin
        # effective M-M constant for Kc in the presence of O2
        Kc * (1 + Om / Ko)
    end ~ track(u"μbar")

    # Rd25: Values in Kim (2006) are for 31C, and the values here are normalized for 25C. SK
    Rd25: dark_respiration_at_25 => 2 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    Ear: activation_energy_for_respiration => 39.8 ~ preserve(u"kJ/mol", parameter)
    Rd(Rd25, kT, Ear): dark_respiration => begin
        Rd25 * kT(Ear)
    end ~ track(u"μmol/m^2/s")
    Rm(Rd) => 0.5Rd ~ track(u"μmol/m^2/s")
end

@system C4Rate(C4c, C4j, C4r) begin
    # Enzyme limited A (Rubisco or PEP carboxylation)
    Ac1(Vp, gbs, Cm, Rm) => (Vp + gbs*Cm - Rm) ~ track(u"μmol/m^2/s" #= CO2 =#)
    Ac2(Vcmax, Rd) => (Vcmax - Rd) ~ track(u"μmol/m^2/s" #= CO2 =#)
    Ac(Ac1, Ac2): enzyme_limited_photosynthesis_rate => begin
        #Ac1 = max(0, Ac1) # prevent Ac1 from being negative Yang 9/26/06
        min(Ac1, Ac2)
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    # x: Partitioning factor of J, yield maximal J at this value
    x: electron_transport_partitioning_factor => 0.4 ~ preserve(parameter)
    # Light and electron transport limited A mediated by J
    Aj1(x, J, Rm, gbs, Cm) => (x * J/2 - Rm + gbs*Cm) ~ track(u"μmol/m^2/s" #= CO2 =#)
    Aj2(x, J, Rd) => (1-x) * J/3 - Rd ~ track(u"μmol/m^2/s" #= CO2 =#)
    Aj(Aj1, Aj2): transport_limited_photosynthesis_rate => begin
        min(Aj1, Aj2)
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    # smoothing the transition between Ac and Aj
    β: photosynthesis_transition_factor => 0.99 ~ preserve(parameter)
    A_net(Ac, Aj, β): net_photosynthesis => begin
        x = A_net
        a = β
        b = -(Ac+Aj)
        c = Ac*Aj
        a*x^2 + b*x + c ⩵ 0
    end ~ solve(pick=:minimum, u"μmol/m^2/s")

    A_gross(A_net, Rd): gross_photosynthesis => begin
        A_gross = A_net + Rd
        # gets negative when PFD = 0, Rd needs to be examined, 10/25/04, SK
        #max(A_gross, zero(A_gross))
    end ~ track(u"μmol/m^2/s" #= CO2 =#)
end

@system C4(C4Rate) begin
    #FIXME: currently not used variables

    # FIXME are they even used?
    # beta_ABA => 1.48e2 # Tardieu-Davies beta, Dewar (2002) Need the references !?
    # delta => -1.0
    # alpha_ABA => 1.0e-4
    # lambda_r => 4.0e-12 # Dewar's email
    # lambda_l => 1.0e-12
    # K_max => 6.67e-3 # max. xylem conductance (mol m-2 s-1 MPa-1) from root to leaf, Dewar (2002)

    # alpha: fraction of PSII activity in the bundle sheath cell, very low for NADP-ME types
    α: bundle_sheath_PSII_activity_fraction => 0.0001 ~ preserve(parameter)
    # Bundle sheath O2 partial pressure, mbar
    Os(A_net, gbs, Om, α): bundle_sheath_o2 => begin
        α * A_net / (0.047gbs) + Om
    end ~ track(u"mbar")

    Cbs(A_net, Vp, Cm, Rm, gbs): bundle_sheath_co2 => begin
        Cm + (Vp - A_net - Rm) / gbs # Bundle sheath CO2 partial pressure, ubar
    end ~ track(u"μbar")

    # half the reciprocal of rubisco specificity, to account for O2 dependence of CO2 comp point,
    # note that this become the same as that in C3 model when multiplied by [O2]
    Γ1 => 0.193 ~ preserve(parameter)
    Γ★(Γ1, Os) => Γ1 * Os ~ track(u"μbar")
    Γ(Rd, Km, Vcmax, Γ★): co2_compensation_point => begin
        (Rd*Km + Vcmax*Γ★) / (Vcmax - Rd)
    end ~ track(u"μbar")
end


text:
# Check for expected base covariance
@test length(getcovariance(m)) == 0
@test length(getcovariance(m2)) == 0

# Test the get/set functions
setcovariance(m, Matrix(I, flex_data.numRVs, flex_data.numRVs))
@test all(getcovariance(m) .== Matrix(I, flex_data.numRVs, flex_data.numRVs))

# Check the safeguards and reset the covariance matrix
@test_throws ErrorException setcovariance(m, [1 1])
@test_throws ErrorException setcovariance(m, [1 2; 3 4])
@test_throws ErrorException setcovariance(m, [-2 0; 0 -3])
flex_data.covariance = Matrix(undef, 0, 0)


text:
"""
Inherit this class for muzero to play
"""
struct Game
    board::Array <: Any
    player::Int
end

"""
Apply action to the game.

Args:
    action : action of the action_space to take.

Returns:
    The new observation, the reward and a boolean if the game has ended.
"""
function execute_step(action)
    nothing
end

"""
Return the current player.

Returns:
    The current player, it should be an element of the players list in the config.
"""
function to_play()
    return 0
end

"""
Should return the legal actions at each turn, if it is not available, it can return
the whole action space. At each turn, the game have to be able to handle one of returned actions.

For complex game where calculating legal moves is too long, the idea is to functionize the legal actions
equal to the action space but to return a negative reward if the action is illegal.

Returns:
    An array of integers, subset of the action space.
"""
function legal_actions()
    nothing
end

"""
Reset the game for a new game.

Returns:
    Initial observation of the game.
"""
function reset_game()
    nothing
end

"""
Properly close the game.
"""
function close_game()
    nothing
end

"""
For multiplayer games, ask the user for a legal action
and return the corresponding action number.

Returns:
An integer from the action space.
"""
function human_to_action()
    choice = input("Enter the action to play for the player {to_play()}: ")
    while int(choice) ∉ legal_actions()
        choice = input("Ilegal action. Enter another action : ")
    end
    return int(choice)
end

"""
Hard coded agent that MuZero faces to assess his progress in multiplayer games.
It doesn't influence training

Returns:
Action as an integer to take in the current game state
"""
function expert_agent()
    error("unimplemented")
end

"""
Convert an action number to a string representing the action.

Args:
action_number: an integer from the action space.

Returns:
String representing the action.
"""
function action_to_string(action_number)
    return "action_number"
end

# Extras

function get_observation(game::Game)
    return observation
end

function have_winner()::Bool end

function expert_action()

end

"""
Display the game observation.
"""
function render_game(game::Game)
    println(game.board)
end


text:
#!/usr/bin/env julia

#load path to qjulia home directory
push!(LOAD_PATH, joinpath(@__DIR__, "..", "core"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "libs/quda-routines"))
push!(LOAD_PATH, joinpath(@__DIR__, "..", "main/fields"))

import QJuliaGrid
import QJuliaFields
import QJuliaFieldUtils
import QJuliaBlas
import QJuliaReduce
import QJuliaUtils
import QJuliaEnums
import QJuliaInterface
import QJuliaGaugeUtils
import QJuliaComms
import QJuliaSolvers
import QUDARoutines

using Random
using LinearAlgebra
using MPI

##############################################################################################

[QJuliaUtils.gridsize_from_cmdline[i] = 1 for i = 1:length(QJuliaUtils.gridsize_from_cmdline)]
QJuliaUtils.get_rank_order("col")

#initialize MPI
MPI.Init()

QUDARoutines.initCommsGridQuda_qj(length(QJuliaUtils.gridsize_from_cmdline), QJuliaUtils.gridsize_from_cmdline, QJuliaUtils.lex_rank_from_coords_t_c, C_NULL)

QUDARoutines.initQuda_qj(0)

Random.seed!(2019)

const lx = 32
const ly = 32
const lz = 32
const lt = 32
const ls = 1

data_type = ComplexF32
data_prec_type = ComplexF32
data_sloppy_type = ComplexF32

grid_desc = QJuliaGrid.QJuliaGridDescr_qj{data_type}(QJuliaEnums.QJULIA_CPU_FIELD_LOCATION, 0, (lx,ly,lz,lt,ls))

cs_in = QJuliaFields.CreateColorSpinor(grid_desc)
cs_ou = QJuliaFields.CreateColorSpinor(grid_desc)

QJuliaFieldUtils.gen_random_spinor!(cs_in)

gauge_param = QJuliaInterface.QJuliaGaugeParam_qj()
gauge_param.X = (lx, ly, lz, lt)
gauge_param.cpu_prec   = data_type == ComplexF64 ? QJuliaEnums.QJULIA_DOUBLE_PRECISION : QJuliaEnums.QJULIA_SINGLE_PRECISION
gauge_param.t_boundary = QJuliaEnums.QJULIA_ANTI_PERIODIC_T
gauge_param.gtype      = QJuliaEnums.QJULIA_WILSON_LINKS
gauge_param.anisotropy = 2.38

gauge_param.cuda_prec                     = data_type == ComplexF64 ? QJuliaEnums.QJULIA_DOUBLE_PRECISION : QJuliaEnums.QJULIA_SINGLE_PRECISION
gauge_param.reconstruct                   = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.cuda_prec_sloppy              = data_sloppy_type == ComplexF64 ? QJuliaEnums.QJULIA_DOUBLE_PRECISION : QJuliaEnums.QJULIA_SINGLE_PRECISION
gauge_param.reconstruct_sloppy            = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.cuda_prec_precondition        = data_prec_type == ComplexF64 ? QJuliaEnums.QJULIA_DOUBLE_PRECISION : QJuliaEnums.QJULIA_SINGLE_PRECISION
gauge_param.reconstruct_precondition      = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.reconstruct_refinement_sloppy = QJuliaEnums.QJULIA_RECONSTRUCT_12
gauge_param.cuda_prec_refinement_sloppy   = QJuliaEnums.QJULIA_HALF_PRECISION

println("======= Gauge parameters =======")
QUDARoutines.printQudaGaugeParam_qj(gauge_param)

gauge_field = QJuliaFields.CreateGaugeField(grid_desc)

QJuliaFieldUtils.construct_gauge_field!(gauge_field, 1, gauge_param)

inv_param = QJuliaInterface.QJuliaInvertParam_qj()
inv_param.residual_type = QJuliaEnums.QJULIA_L2_RELATIVE_RESIDUAL

println("======= Invert parameters =======")
QUDARoutines.printQudaInvertParam_qj(inv_param)

gauge_param.gtype       = QJuliaEnums.QJULIA_SU3_LINKS 		#currently cannot set QJULIA_WILSON_LINKS (=QJULIA_SU3_LINKS)  for QUDA

x_face_size = gauge_param.X[2]*gauge_param.X[3]*Int(gauge_param.X[4]/2);
y_face_size = gauge_param.X[1]*gauge_param.X[3]*Int(gauge_param.X[4]/2);
z_face_size = gauge_param.X[1]*gauge_param.X[2]*Int(gauge_param.X[4]/2);
t_face_size = gauge_param.X[1]*gauge_param.X[2]*Int(gauge_param.X[3]/2);

gauge_param.ga_pad = max(x_face_size, y_face_size, z_face_size, t_face_size);

QUDARoutines.loadGaugeQuda_qj(gauge_field.v, gauge_param)

#Check plaquette
plaq = Array{Cdouble, 1}(undef, 3)
QUDARoutines.plaqQuda_qj(plaq)
println("Computed plaquette is ", plaq[1], ", (spatial = ",  plaq[2], ", temporal = ", plaq[3], ")")

mass = -0.9

inv_param.mass = mass
inv_param.kappa = 1.0 / (2.0 * (1 + 3/gauge_param.anisotropy + mass))
inv_param.maxiter = 500
inv_param.tol = 5e-8

inv_param.cpu_prec = gauge_param.cpu_prec
inv_param.cuda_prec = gauge_param.cuda_prec
inv_param.cuda_prec_sloppy = gauge_param.cuda_prec_sloppy
inv_param.cuda_prec_precondition = QJuliaEnums.QJULIA_HALF_PRECISION
inv_param.solution_type = QJuliaEnums.QJULIA_MATPC_SOLUTION
inv_param.inv_type = QJuliaEnums.QJULIA_PIPEPCG_INVERTER
#inv_param.inv_type = QJuliaEnums.QJULIA_PCG_INVERTER

println("Kappa = ",  inv_param.kappa)

mdagm(out, inp) = QUDARoutines.MatDagMatQuda_qj(out, inp, inv_param)
Doe(out, inp)   = QUDARoutines.dslashQuda_qj(out, inp, inv_param, QJuliaEnums.QJULIA_EVEN_PARITY)
Deo(out, inp)   = QUDARoutines.dslashQuda_qj(out, inp, inv_param, QJuliaEnums.QJULIA_ODD_PARITY )

# Setup preconditioner
precond_param = QJuliaInterface.QJuliaInvertParam_qj()

precond_param.residual_type            = QJuliaEnums.QJULIA_L2_RELATIVE_RESIDUAL
precond_param.inv_type                 = QJuliaEnums.QJULIA_PCG_INVERTER
#precond_param.inv_type                 = QJuliaEnums.QJULIA_LANMR_INVERTER #wroks for naive and fails for pipelined
precond_param.inv_type                 = QJuliaEnums.QJULIA_INVALID_INVERTER
precond_param.dslash_type_precondition = QJuliaEnums.QJULIA_WILSON_DSLASH
precond_param.kappa                    = 1.0 / (2.0 * (1 + 3/gauge_param.anisotropy + mass))
precond_param.cuda_prec                = data_prec_type == ComplexF64 ? QJuliaEnums.QJULIA_DOUBLE_PRECISION : QJuliaEnums.QJULIA_SINGLE_PRECISION
precond_param.cuda_prec_sloppy         = precond_param.cuda_prec
precond_param.cuda_prec_precondition   = precond_param.cuda_prec
precond_param.solution_type            = QJuliaEnums.QJULIA_MATPC_SOLUTION
precond_param.maxiter                  = precond_param.inv_type == QJuliaEnums.QJULIA_PCG_INVERTER ? 30 : 6
precond_param.Nsteps    	       = 1

mdagmPre(out, inp)  = QUDARoutines.MatDagMatQuda_qj(out, inp, precond_param)

pre_solv_param = QJuliaSolvers.QJuliaSolverParam_qj()

pre_solv_param.dtype = data_prec_type == ComplexF64 ? Float64 : Float32

pre_solv_param.inv_type  = precond_param.inv_type
pre_solv_param.tol       = 1e-2
#
pre_solv_param.maxiter   = precond_param.maxiter
pre_solv_param.Nsteps    = 1
pre_solv_param.global_reduction = false

K(out, inp) = QJuliaSolvers.solve(out, inp, mdagmPre, mdagmPre, pre_solv_param)

x_even = QJuliaFields.Even(cs_ou)
x_odd  = QJuliaFields.Odd(cs_ou)

b_even = QJuliaFields.Even(cs_in)
b_odd  = QJuliaFields.Odd(cs_in)

#Auxiliary field
tmp = QJuliaFields.CreateColorSpinor(spinor_field_desc; NSpin=4)

t_even = QJuliaFields.Even(tmp)
t_odd  = QJuliaFields.Odd(tmp)

#prepare source/solution:
if inv_param.matpc_type == QJuliaEnums.QJULIA_MATPC_EVEN_EVEN
# src = b_e + k D_eo b_o
Deo(t_even.v, b_odd.v)
x_odd.v .=@. b_even.v + inv_param.kappa*t_even.v
end

#
x_odd2 = dot(x_odd.v, x_odd.v)
#
println("Initial source norm: ", sqrt(real(x_odd2)), " ,source norm2 ", x_odd2, " , requested tolerance: ", inv_param.tol)

solv_param = QJuliaSolvers.QJuliaSolverParam_qj()

solv_param.dtype = data_type == ComplexF64 ? Float64 : Float32
solv_param.dtype_sloppy = data_sloppy_type == ComplexF64 ? Float64 : Float32
solv_param.dtype_precondition = data_prec_type == ComplexF64 ? Float64 : Float32
# Set up parameters
solv_param.inv_type               = inv_param.inv_type
solv_param.inv_type_precondition  = precond_param.inv_type
solv_param.tol                    = inv_param.tol
#
solv_param.maxiter                = inv_param.maxiter
solv_param.Nsteps                 = 1

#matrix form
#sol = view(reinterpret(cs_ou.field_desc.prec, x_even.v), :, :)
#src = view(reinterpret(cs_in.field_desc.prec, x_odd.v), :, :)

#vector form
sol = view(reinterpret(cs_ou.field_desc.prec, x_even.v), :)
src = view(reinterpret(cs_in.field_desc.prec, x_odd.v), :)

if precond_param.inv_type != QJuliaEnums.QJULIA_INVALID_INVERTER
  QJuliaSolvers.solve(sol, src, mdagm, mdagm, solv_param, K)
else
  QJuliaSolvers.solve(sol, src, mdagm, mdagm, solv_param)
end

#compute true residual:
r = t_odd.v
mdagm(r, x_even.v)
r  .=@. x_odd.v - r
r2 = dot(r, r)
println("True residual: ", sqrt(real(r2)))

#reconstruct source/solution:
if inv_param.matpc_type == QJuliaEnums.QJULIA_MATPC_EVEN_EVEN
# x_o = b_o + k D_oe x_e
Doe(t_odd.v, x_even.v)
x_odd.v .=@. b_odd.v + inv_param.kappa*t_odd.v
end

QUDARoutines.endQuda_qj()
MPI.Finalize()

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
"""
    add_reactions!(model::StandardModel, rxns::Vector{Reaction})

Add `rxns` to `model` based on reaction `id`.
"""
function add_reactions!(model::StandardModel, rxns::Vector{Reaction})
    for rxn in rxns
        model.reactions[rxn.id] = rxn
    end
    nothing
end

"""
    add_reaction!(model::StandardModel, rxn::Reaction)

Add `rxn` to `model` based on reaction `id`.
"""
add_reaction!(model::StandardModel, rxn::Reaction) = add_reactions!(model, [rxn])

"""
    add_metabolites!(model::StandardModel, mets::Vector{Metabolite})

Add `mets` to `model` based on metabolite `id`.
"""
function add_metabolites!(model::StandardModel, mets::Vector{Metabolite})
    for met in mets
        model.metabolites[met.id] = met
    end
    nothing
end

"""
    add_metabolite!(model::StandardModel, met::Metabolite)

Add `met` to `model` based on metabolite `id`.
"""
add_metabolite!(model::StandardModel, met::Metabolite) = add_metabolites!(model, [met])

"""
    add_genes!(model::StandardModel, genes::Vector{Gene})

Add `genes` to `model` based on gene `id`.
"""
function add_genes!(model::StandardModel, genes::Vector{Gene})
    for gene in genes
        model.genes[gene.id] = gene
    end
    nothing
end

"""
    add_gene!(model::StandardModel, genes::Gene)

Add `gene` to `model` based on gene `id`.
"""
add_gene!(model::StandardModel, gene::Gene) = add_genes!(model, [gene])

"""
    @add_reactions!(model::Symbol, ex::Expr)

Shortcut to add multiple reactions and their lower and upper bounds

Call variants
-------------
```
@add_reactions! model begin
    reaction_name, reaction
end

@add_reactions! model begin
    reaction_name, reaction, lower_bound
end

@add_reactions! model begin
    reaction_name, reaction, lower_bound, upper_bound
end
```

Examples
--------
```
@add_reactions! model begin
    "v1", nothing → A, 0, 500
    "v2", A ↔ B + C, -500
    "v3", B + C → nothing
end
```
"""
macro add_reactions!(model::Symbol, ex::Expr)
    model = esc(model)
    all_reactions = Expr(:block)
    for line in MacroTools.striplines(ex).args
        args = line.args
        id = esc(args[1])
        reaction = esc(args[2])
        push!(all_reactions.args, :(r = $reaction))
        push!(all_reactions.args, :(r.id = $id))
        if length(args) == 3
            lb = args[3]
            push!(all_reactions.args, :(r.lb = $lb))
        elseif length(args) == 4
            lb = args[3]
            ub = args[4]
            push!(all_reactions.args, :(r.lb = $lb))
            push!(all_reactions.args, :(r.ub = $ub))
        end
        push!(all_reactions.args, :(add_reaction!($model, r)))
    end
    return all_reactions
end

"""
    remove_genes!(
        model::StandardModel,
        ids::Vector{String};
        knockout_reactions::Bool = false,
    )

Remove all genes with `ids` from `model`. If `knockout_reactions` is true, then also
constrain reactions that require the genes to function to carry zero flux.

# Example
```
remove_genes!(model, ["g1", "g2"])
```
"""
function remove_genes!(
    model::StandardModel,
    gids::Vector{String};
    knockout_reactions::Bool = false,
)
    if knockout_reactions
        rm_reactions = String[]
        for (rid, r) in model.reactions
            if !isnothing(r.grr) &&
               all(any(in.(gids, Ref(conjunction))) for conjunction in r.grr)
                push!(rm_reactions, rid)
            end
        end
        pop!.(Ref(model.reactions), rm_reactions)
    end
    pop!.(Ref(model.genes), gids)
    nothing
end

"""
    remove_gene!(
        model::StandardModel,
        id::Vector{String};
        knockout_reactions::Bool = false,
    )

Remove gene with `id` from `model`. If `knockout_reactions` is true, then also
constrain reactions that require the genes to function to carry zero flux.

# Example
```
remove_gene!(model, "g1")
```
"""
remove_gene!(model::StandardModel, gid::String; knockout_reactions::Bool = false) =
    remove_genes!(model, [gid]; knockout_reactions = knockout_reactions)

@_change_bounds_fn StandardModel String inplace begin
    isnothing(lower) || (model.reactions[rxn_id].lb = lower)
    isnothing(upper) || (model.reactions[rxn_id].ub = upper)
    nothing
end

@_change_bounds_fn StandardModel String inplace plural begin
    for (i, l, u) in zip(rxn_ids, lower, upper)
        change_bound!(model, i, lower = l, upper = u)
    end
end

@_change_bounds_fn StandardModel String begin
    change_bounds(model, [rxn_id], lower = [lower], upper = [upper])
end

@_change_bounds_fn StandardModel String plural begin
    n = copy(model)
    n.reactions = copy(model.reactions)
    for i in rxn_ids
        n.reactions[i] = copy(n.reactions[i])
    end
    change_bounds!(n, rxn_ids, lower = lower, upper = upper)
    return n
end

@_remove_fn reaction StandardModel String inplace begin
    if !(reaction_id in reactions(model))
        @_models_log @info "Reaction $reaction_id not found in model."
    else
        delete!(model.reactions, reaction_id)
    end
    nothing
end

@_remove_fn reaction StandardModel String inplace plural begin
    remove_reaction!.(Ref(model), reaction_ids)
    nothing
end

@_remove_fn reaction StandardModel String begin
    remove_reactions(model, [reaction_id])
end

@_remove_fn reaction StandardModel String plural begin
    n = copy(model)
    n.reactions = copy(model.reactions)
    remove_reactions!(n, reaction_ids)
    return n
end

@_remove_fn metabolite StandardModel String inplace begin
    remove_metabolites!(model, [metabolite_id])
end

@_remove_fn metabolite StandardModel String inplace plural begin
    !all(in.(metabolite_ids, Ref(metabolites(model)))) &&
        @_models_log @info "Some metabolites not found in model."
    remove_reactions!(
        model,
        [
            rid for (rid, rn) in model.reactions if
            any(haskey.(Ref(rn.metabolites), metabolite_ids))
        ],
    )
    delete!.(Ref(model.metabolites), metabolite_ids)
    nothing
end

@_remove_fn metabolite StandardModel String begin
    remove_metabolites(model, [metabolite_id])
end

@_remove_fn metabolite StandardModel String plural begin
    n = copy(model)
    n.reactions = copy(model.reactions)
    n.metabolites = copy(model.metabolites)
    remove_metabolites!(n, metabolite_ids)
    return n
end

"""
    change_objective!(
        model::StandardModel,
        rxn_ids::Vector{String};
        weights = ones(length(rxn_ids)),
    )

Change the objective for `model` to reaction(s) with `rxn_ids`, optionally specifying their `weights`. By default,
assume equal weights. If no objective exists in model, sets objective.
"""
function change_objective!(
    model::StandardModel,
    rxn_ids::Vector{String};
    weights = ones(length(rxn_ids)),
)
    all(!haskey(model.reactions, rid) for rid in rxn_ids) &&
        throw(DomainError(rxn_ids, "Some reaction ids were not found in model."))

    for (_, rxn) in model.reactions # reset to zero
        rxn.objective_coefficient = 0.0
    end
    for (k, rid) in enumerate(rxn_ids)
        model.reactions[rid].objective_coefficient = weights[k]
    end
end

change_objective!(model::StandardModel, rxn_id::String) = change_objective!(model, [rxn_id])


text:
# Test fluent functionality
path = joinpath(dirname(pathof(PDDL)), "..", "test", "fluents")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == Symbol("zeno-travel")
@test domain.functions[:distance] == pddl"(distance ?c1 ?c2)"
@test domain.functypes[:fuel] == [:aircraft]

# Test for static functions
@test pddl"(capacity ?a)" in get_static_functions(domain)
@test length(get_static_functions(domain)) == 5

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.metric == (-1, pddl"(+ (* 4 (total-time)) (* 5 (total-fuel-used)))")

state = init_state(problem)

# Person 1 boards plane 1
state = execute(pddl"(board person1 plane1 city0)", state, domain)
@test state[domain, pddl"(onboard plane1)"] == 1

# Plane 1 flies from city 0 to city 2
state = execute(pddl"(fly plane1 city0 city2)", state, domain)
@test state[domain, pddl"(total-fuel-used)"] == 3100

# Person 1 debarks at city 2
state = execute(pddl"(debark person1 plane1 city2)", state, domain)
@test state[domain, pddl"(at person1 city2)"] == true

# Plane 1 refuels at city 2
state = execute(pddl"(refuel plane1 city2)", state, domain)
@test state[domain, pddl"(fuel plane1)"] == 10232

# Person 2 boards plane 1 at city2
state = execute(pddl"(board person2 plane1 city2)", state, domain)
@test state[domain, pddl"(in person2 plane1)"] == true

# Plane 1 flies from city 2 to city 0
state = execute(pddl"(fly plane1 city2 city0)", state, domain)
@test state[domain, pddl"(total-fuel-used)"] == 6200

# Person 2 debarks at city 0
state = execute(pddl"(debark person2 plane1 city0)", state, domain)
@test state[domain, pddl"(at person2 city0)"] == true

# Plane 1 refuels at city 0
state = execute(pddl"(refuel plane1 city0)", state, domain)
@test state[domain, pddl"(fuel plane1)"] == 10232

# Plane 1 zooms from city 0 to city 1
state = execute(pddl"(zoom plane1 city0 city1)", state, domain)
@test state[domain, pddl"(total-fuel-used)"] == 16370

# The whole plan took 9 steps
@test state[domain, pddl"(total-time)"] == 9

# Check if goal is satisfied
@test satisfy(problem.goal, state, domain)[1] == true

# Test execution of entire plan
state = init_state(problem)
plan = @pddl(
    "(board person1 plane1 city0)",
    "(fly plane1 city0 city2)",
    "(debark person1 plane1 city2)",
    "(refuel plane1 city2)",
    "(board person2 plane1 city2)",
    "(fly plane1 city2 city0)",
    "(debark person2 plane1 city0)",
    "(refuel plane1 city0)",
    "(zoom plane1 city0 city1)"
)
state = execute(plan, state, domain)
@test satisfy(problem.goal, state, domain)[1] == true


text:
function qp1_test()
    #=
        min x^2 + xy + y^2 + yz + z^2
    s.t.
        x + 2y + 3z >= 4 (c1)
        x +  y      >= 1 (c2)
        x,y \in R
    =#
    model = TestModel{Float64}()

    v = MOI.add_variables(model, 3)

    cf1 =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], v), 0.0)

    c1 = MOI.add_constraint(model, cf1, MOI.GreaterThan(4.0))
    c2 = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [v[1], v[2]]),
            0.0,
        ),
        MOI.GreaterThan(1.0),
    )

    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(
            [2.0, 1.0, 2.0, 1.0, 2.0],
            v[[1, 1, 2, 2, 3]],
            v[[1, 2, 2, 3, 3]],
        ),
        MOI.ScalarAffineTerm{Float64}[],
        0.0,
    )
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end

function qp2_test()
    #=
        min 2 x^2 + y^2 + xy + x + y + 1
    s.t.
        x, y >= 0
        x + y = 1
    =#

    model = TestModel{Float64}()

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)

    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
            0.0,
        ),
        MOI.EqualTo(1.0),
    )

    vc1 = MOI.add_constraint(model, x, MOI.GreaterThan(0.0))
    vc2 = MOI.add_constraint(model, y, MOI.GreaterThan(0.0))

    obj = MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.([4.0, 2.0, 1.0], [x, y, x], [x, y, y]),
        MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]),
        1.0,
    )

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
        obj,
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return model
end


text:
#=
Author: Stephan Scholz
Year: 2021

This file contains a machine learning approach (SciML) for the van der Pol oscillator
=#

# ODE function of van der Pol oscillator
function vdp(dx, x, p, t)

    μ = p[1]

    dx[1] = μ * (x[1] - 1/3 * x[1]^3 - x[2] )
    dx[2] = 1/μ * x[1]
end

tspan = (0.0, 100.0) # time span
x0    = [0.5, 0.5]   # initial values
p     = [1.0]        # parameter: μ
tsteps= 0.1          # saved time steps

using DifferentialEquations
prob = ODEProblem( vdp, x0, tspan, p )      # Build ODEProblem
sol = solve(prob, Tsit5(), saveat = tsteps) # Solve ODEProblem

# using Plots
# plot(sol, title="Solution of the van der Pol oscillator")

# MACHINE LEARNING
function loss(p)
    sol_pred = solve(prob, Tsit5(), p=p, saveat = tsteps)
    loss = sum(abs2, sol - sol_pred)
    return loss, sol_pred
end

pred_param = [] # Stores all found parameters

callback = function (p, l, pred)
    display(l)
    append!(pred_param, p)

    # Tell sciml_train to not halt the optimization. If return true, then
    # optimization stops.
    return false
end

using DiffEqFlux

p_ml = [0.1] # Initial parameter
loss(p_ml)   # Test whether loss function works

# Run ML training
result_ode = DiffEqFlux.sciml_train(loss, p_ml, ADAM(0.1), cb = callback, maxiters=100)

# Build and plot Cost function
x_min = 0.1
x_max = 2
x_step = 0.01
x_range = range(x_min,x_max,step=x_step)
Nlength = length(x_range)

loss_results = zeros(Nlength)
idx = 1;
for i in x_range
   loss_results[idx],_ = loss([i])
   idx += 1;
end

# Cost function
plot(x_range, loss_results, label="loss", title="Cost function", xaxis=("Parameter μ"))

# Save Machine Learning results as gif
@gif for pp in pred_param
    sol_remade = solve(remake(prob, p=[pp]), Tsit5(), saveat=tsteps)

    # Phase portrait
    pp_text = string( "μ = ", string(round(pp; digits=3)) )
    plot(sol_remade[1,:], sol_remade[2,:], annotations=(0.1, 0, Plots.text(pp_text, :left)), title="Phase portrait of van der Pol oscillator", label=false)

    # Oscillation
    # plot(sol_remade, title="Machine Learning of van der Pol oscillator", label=["x1" "x2"], legend=:bottomright)

end


text:
# CodeMatch: Pattern -> matching code
module CodeMatch
using Graph, Toivo
export code_match, code_pat

# ---- sequence: construct an evaluation order --------------------------------

type Seq
    p::Pattern
    evaluated::Set{Node}
    seq::Vector{Node}

    Seq(p::Pattern) = new(p, Set{Node}(), Node[])
end

function sequence(p::Pattern)
    c = Seq(p)
    for g in guardsof(p);  sequence!(c, g);  end
    c.seq
end

function sequence!(c::Seq, node::Node)
    if has(c.evaluated, node); return; end

    for dep in depsof(c.p, node);  sequence!(c, dep);  end
    add(c.evaluated, node)
    push(c.seq, node)
end

# ---- code_pred: code the matching predicate of a pattern --------------------

function code_pat(p::Pattern)
    pred, results = code_pred(p)
    code = {}
    for b in p.bindings
        push(code, :($(b.name) = $(evaluate!(results, b.arg))))
    end
    pred, expr(:block, code)
end

function code_pred(p::Pattern)
    seq = sequence(p)
    results = Dict{Node,Any}()
    code = {}
    preds = {}
    for node in seq
        ex = evaluate!(results, node)
        if isa(node, Guard)
            push(preds, isempty(code) ? ex : quote; $(code...); $ex; end)
            code = {}
        else
            if isa(ex, Symbol)
                results[node] = ex
            else
                val = gensym("v")
                push(code, :( $val = $ex ))
                results[node] = val
            end
        end
    end
    if isempty(preds); return quot(true); end
    pred = preds[1]
    for factor in preds[2:end]; pred = :($pred && $factor); end
    pred, results
end

function evaluate!(results::Dict{Node,Any}, node::Node)
    if has(results, node) results[node]
    else                  encode(results, node)
    end
end

encode(c, v::Arg)      = argsym
encode(c, v::TupleRef) = :( $(evaluate!(c,v.arg))[$(v.index)] )
encode(c, g::Egal) = :(is( $(evaluate!(c, g.arg)), $(quot(g.value))))
encode(c, g::Isa)  = :(isa($(evaluate!(c, g.arg)), $(quot(g.typ  ))))

# ---- Old stuff -------------------------------------------------------------

type Ctx
    p::Pattern
    code::Vector
    results::Dict{Node,Any}
    Ctx(p::Pattern) = new(p, {}, Dict{Node,Any}())
end

emit(c::Ctx, ex) = (push(c.code, ex); nothing)
emit_guard(c::Ctx, ex) = emit(c, :( if !$ex; return (false,nothing); end ))

function code_match(p::Pattern)
    c = Ctx(p)
    for g in values(p.guards); evaluate(c, g); end

    bound = Set{Symbol}()
    for b in p.bindings
        if has(bound, b.name)
            emit_guard(c, :( is($(evaluate(c,b.arg)), $(b.name)) ))
        else
            emit(c, :( $(b.name) = $(evaluate(c,b.arg)) ))
            add(bound, b.name)
        end
    end
    quote; $(c.code...); end
end

function evaluate(c::Ctx, node::Node)
    if has(c.results, node)
        c.results[node]
    else
        ex = code_match(c, node)
        if (isa(ex, Symbol) || ex === nothing)
            c.results[node] = ex
        else
            val = gensym("v")
            emit(c, :( $val = $ex ))
            c.results[node] = val
        end
    end
end

evalguardof(c::Ctx, node::Value) = evaluate(c, c.p.guards[node])

code_match(c::Ctx, v::Arg) = argsym
function code_match(c::Ctx, v::TupleRef)
    evalguardof(c, v.arg)
    :( $(evaluate(c,v.arg))[$(v.index)] )
end

code_match(c::Ctx, g::Guard) = emit_guard(c, code_pred(c, g))

code_pred(c::Ctx,g::Egal) = :(is( $(evaluate(c, g.arg)), $(quot(g.value))))
code_pred(c::Ctx,g::Isa)  = :(isa($(evaluate(c, g.arg)), $(quot(g.typ  ))))

end # module