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

using Reexport

using AutoHashEquals
using DataStructures: OrderedDict, IntDisjointSets, union!, find_root!
using FunctionWrappers: FunctionWrapper
import StaticArrays
using StaticArrays: StaticVector, SVector, SizedVector, similar_type

@reexport using ..Sets
using ...GAT, ...Theories
using ...CSetDataStructures, ...Graphs, ..FreeDiagrams, ..Limits, ..Subobjects
import ...Theories: dom, codom, meet, ∧, join, ∨, top, ⊤, bottom, ⊥
import ..Limits: limit, colimit, universal
import ..Subobjects: Subobject, SubobjectLattice
using ..Sets: SetFunctionCallable, SetFunctionIdentity

# Data types
############

# Finite sets
#------------

""" Finite set.

This generic type encompasses the category **FinSet** of finite sets and
functions, through types `FinSet{S} where S <: AbstractSet`, as well as the
skeleton of this category, through the type `FinSet{Int}`. In the latter case,
the object `FinSet(n)` represents the set ``{1,...,n}``.
"""
@auto_hash_equals struct FinSet{S,T} <: SetOb{T}
  set::S
end

FinSet(n::Int) = FinSet{Int,Int}(n)
FinSet(set::S) where {T,S<:AbstractSet{T}} = FinSet{S,T}(set)
FinSet(s::FinSet) = s

Base.iterate(s::FinSet, args...) = iterate(iterable(s), args...)
Base.length(s::FinSet) = length(iterable(s))
Base.in(s::FinSet, elem) = in(s, iterable(s))
iterable(s::FinSet{Int}) = 1:s.set
iterable(s::FinSet{<:AbstractSet}) = s.set

Base.show(io::IO, s::FinSet) = print(io, "FinSet($(s.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...) =
  SetFunctionIdentity((FinSet(arg) for arg in args)...)
FinFunction(f::AbstractVector{Int}; kw...) =
  FinDomFunctionVector(f, FinSet(isempty(f) ? 0 : maximum(f)); kw...)

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

Sets.show_type(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} = SetFunction{Dom,Codom}

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

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

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

""" 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}``.
"""
@auto_hash_equals struct FinDomFunctionVector{T,V<:AbstractVector{T},
    Codom<:SetOb{T}} <: FinDomFunction{Int,FinSet{Int,Int},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]

Base.show(io::IO, f::FinDomFunctionVector) =
  print(io, "FinDomFunction($(f.func), $(dom(f)), $(codom(f)))")

""" Force evaluation of lazy function or relation.
"""
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.compose_impl(f::FinFunctionVector, g::FinDomFunctionVector) =
  FinDomFunctionVector(g.func[f.func], codom(g))

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

# Indexed functions
#------------------

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

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

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

function IndexedFinDomFunction(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
  IndexedFinDomFunction(f, index, codom)
end

Base.show(io::IO, f::IndexedFinDomFunction) =
  print(io, "FinDomFunction($(f.func), $(dom(f)), $(codom(f)), index=true)")

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

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

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

""" The preimage (inverse image) of the value y in the codomain.
"""
preimage(f::SetFunctionIdentity, y) = SVector(y)
preimage(f::FinDomFunction, y) = [ x for x in dom(f) if f(x) == y ]
preimage(f::IndexedFinDomFunction, 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 [`IndexedFinDomFunction`](@ref).
"""
const IndexedFinFunction{V,Index} =
  IndexedFinDomFunction{Int,V,Index,FinSet{Int,Int}}

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

function IndexedFinFunction(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
  IndexedFinDomFunction(f, index, codom)
end

Base.show(io::IO, f::IndexedFinFunction) =
  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.compose_impl(f::Union{FinFunctionVector,IndexedFinFunction},
                  g::Union{FinDomFunctionVector,IndexedFinDomFunction}) =
  FinDomFunctionVector(g.func[f.func], codom(g))

# Limits
########

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

function universal(lim::Terminal{<:FinSet{Int}},
                   cone::SMultispan{0,<:FinSet{Int}})
  ConstantFunction(1, apex(cone), FinSet(1))
end

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::BinaryProduct{<:FinSet{Int}}, cone::Span{<:FinSet{Int}})
  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::Product{<:FinSet{Int}}, cone::Multispan{<:FinSet})
  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::Equalizer{<:FinSet{Int}},
                   cone::SMultispan{1,<:FinSet{Int}})
  ι = collect(incl(lim))
  h = only(cone)
  FinFunction(Int[only(searchsorted(ι, h(i))) for i in dom(h)], length(ι))
end

""" Algorithm for limit of spans or multispans out of finite sets.

In the context of relational databases, such limits are joins.
"""
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 Limit(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
  Limit(cospan, Multispan(map((π,f) -> FinFunction(π, dom(f)), πs, funcs)))
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
  Limit(cospan, Multispan(map((π,f) -> FinFunction(π, length(f)), πs, funcs)))
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)
  Limit(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)

""" Limit of free diagram of FinSets.

See `CompositePullback` for a very similar construction.
"""
struct FinSetFreeDiagramLimit{Ob<:FinSet, Diagram<:AbstractFreeDiagram{Ob},
                              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

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

  # 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
    if nv₁(d) == 1
      Limit(d_original, SMultispan{1}(ιs[1]))
    else
      πs = legs(product(SVector(ob₁(d)...)))
      Limit(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
    Limit(d_original, Multispan(πs))
  end
end

""" 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 ``V₁`` set but a possibly reduced ``V₂``.
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

function limit(d::FreeDiagram{<: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!
  prod = product(ob(d))
  n, πs = length(ob(prod)), legs(prod)
  ι = FinFunction(filter(1:n) do i
    all(begin
          s, t, h = src(d,e), tgt(d,e), hom(d,e)
          h(πs[s](i)) == πs[t](i)
        end for e in edges(d))
    end, n)
  cone = Multispan(dom(ι), [ι⋅πs[i] for i in vertices(d)])
  FinSetFreeDiagramLimit(d, cone, prod, ι)
end

function universal(lim::FinSetFreeDiagramLimit, 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

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

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

function universal(colim::Initial{<:FinSet{Int}},
                   cocone::SMulticospan{0,<:FinSet{Int}})
  FinFunction(Int[], apex(cocone))
end

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{<:FinSet{Int}})
  f, g = cocone
  FinFunction(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{<:FinSet{Int}})
  FinFunction(reduce(vcat, (collect(f) for f in cocone), init=Int[]),
              ob(colim), apex(cocone))
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,<:FinSet{Int}})
  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
      @assert q[j] == h(i) "Quotient map out of coequalizer is ill-defined"
    end
  end
  @assert all(i > 0 for i in q) "Projection map is not surjective"
  FinFunction(q, codom(h))
end

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

""" Colimit of free diagram of FinSets.

See `CompositePushout` for a very similar construction.
"""
struct FinSetFreeDiagramColimit{Ob<:FinSet, Diagram<:AbstractFreeDiagram{Ob},
                                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) ])
  FinSetFreeDiagramColimit(d, cocone, coprod, π)
end

function colimit(d::FreeDiagram{<:FinSet{Int}})
  # Uses the general formula for colimits in Set (Leinster, 2014, Basic Category
  # Theory, Example 5.2.16).
  coprod = coproduct(ob(d))
  n, ιs = length(ob(coprod)), legs(coprod)
  sets = IntDisjointSets(n)
  for e in edges(d)
    s, t, h = src(d,e), tgt(d,e), hom(d,e)
    for i in dom(h)
      union!(sets, ιs[s](i), ιs[t](h(i)))
    end
  end
  π = quotient_projection(sets)
  cocone = Multicospan(codom(π), [ ιs[i]⋅π for i in vertices(d) ])
  FinSetFreeDiagramColimit(d, cocone, coprod, π)
end

function universal(colim::FinSetFreeDiagramColimit,
                   cocone::Multicospan{<:FinSet{Int}})
  h = universal(colim.coprod, cocone)
  pass_to_quotient(colim.proj, h)
end

# Subsets
#########

""" Subset of a finite set represented as an inclusion map.
"""
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}

function Subobject(X::FinSet, pred::AbstractBoolVector)
  length(pred) == length(X) ||
    error("Size of predicate $pred does not equal size of object $X")
  Subobject(X, findall(pred))
end
function SubFinSet(pred::AbstractBoolVector)
  Subobject(FinSet(length(pred)), findall(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(as_predicate(A) .& as_predicate(B))
join(A::SubFinSet{Int}, B::SubFinSet{Int}, ::SubOpBoolean) =
  SubFinSet(as_predicate(A) .| as_predicate(B))
top(X::FinSet{Int}, ::SubOpBoolean) = SubFinSet(trues(length(X)))
bottom(X::FinSet{Int}, ::SubOpBoolean) = SubFinSet(falses(length(X)))

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

end


text:
"""
    is_leap_year(year)

Return `true` if `year` is a leap year in the gregorian calendar.

on every year that is evenly divisible by 4
  except every year that is evenly divisible by 100
    unless the year is also evenly divisible by 400

;; from a Clojure program I wrote - first implementation
(defn leap-year? [year]
 (cond
   (zero? (mod year 400))
   true
   (and (zero? (mod year 4)) (not (zero? (mod year 100))))
   true
   :else
   false))

;; 2nd look - this one I translated into Julia
(defn leap-year? [year]
  (or (zero? (mod year 400))

      (and (zero? (mod year 4))
        (not (zero? (mod year 100))))))

"""
function is_leap_year(year::Int64)
  ((year % 400) == 0) || (((year % 4) == 0) && ((year % 100) != 0))
end

ext: jl

max_stars_repo_name: jawuku/julia

max_issues_repo_name: jawuku/julia

max_forks_repo_name: jawuku/julia


text:
include("lib_for_robot.jl")
function task3(r)
    x,y=tostart(r)
    while isborder(r,Ost)==0
        putmarks(r,Nord)
        putmarker!(r)
        move!(r,Ost)
        putmarks(r,Sud)
        putmarker!(r)
    end
    tostart(r)
    for i in 1:x
        move!(r,Ost)
    end
    for i in 1:y
        move!(r,Nord)
    end
end


text:
#=
    H2 tank storage modelling
 =#

mutable struct H2Tank  <: AbstractStorage
     # Paramètres
     α_p_ch::Float64
     α_p_dch::Float64
     η_ch::Float64
     η_dch::Float64
     η_self::Float64
     α_soc_min::Float64
     α_soc_max::Float64
     lifetime::Int64
     bounds::NamedTuple{(:lb, :ub), Tuple{Float64, Float64}}
     # Initial conditions
     Erated_ini::Float64
     soc_ini::Float64
     soh_ini::Float64
     # Variable
     Erated::AbstractArray{Float64,2}
     carrier::Hydrogen
     soc::AbstractArray{Float64,3}
     # Eco
     cost::AbstractArray{Float64,2}
     # Inner constructor
     H2Tank(; α_p_ch = 1.5,
        α_p_dch = 1.5,
        η_ch = 1.,
        η_dch = 1.,
        η_self = 0.,
        α_soc_min = 0.,
        α_soc_max = 1.,
        lifetime = 25,
        bounds = (lb = 0., ub = 10000.),
        Erated_ini = 1e-6,
        soc_ini = 0.5,
        soh_ini = 1.) =
        new(α_p_ch, α_p_dch, η_ch, η_dch, η_self, α_soc_min, α_soc_max, lifetime, bounds, Erated_ini, soc_ini, soh_ini)
end

### Preallocation
function preallocate!(h2tank::H2Tank, nh::Int64, ny::Int64, ns::Int64)
    h2tank.Erated = convert(SharedArray,zeros(ny+1, ns)) ; h2tank.Erated[1,:] .= h2tank.Erated_ini
    h2tank.carrier = Hydrogen()
    h2tank.carrier.power = convert(SharedArray,zeros(nh, ny, ns))
    h2tank.soc = convert(SharedArray,zeros(nh+1, ny+1, ns)) ; h2tank.soc[1,1,:] .= h2tank.soc_ini
    h2tank.cost = convert(SharedArray,zeros(ny, ns))
    return h2tank
end

### Operation dynamic
function compute_operation_dynamics!(h::Int64, y::Int64, s::Int64, h2tank::H2Tank, decision::Float64, Δh::Int64)
     h2tank.soc[h+1,y,s], h2tank.carrier.power[h,y,s] = compute_operation_dynamics(h2tank, (Erated = h2tank.Erated[y,s], soc = h2tank.soc[h,y,s]), decision, Δh)
end

function compute_operation_dynamics(h2tank::H2Tank, state::NamedTuple{(:Erated, :soc), Tuple{Float64, Float64}}, decision::Float64, Δh::Int64)
     # Control power constraint and correction
     power_dch = max(min(decision, h2tank.α_p_dch * state.Erated, h2tank.η_dch * (state.soc * (1. - h2tank.η_self * Δh) - h2tank.α_soc_min) * state.Erated / Δh), 0.)
     power_ch = min(max(decision, -h2tank.α_p_ch * state.Erated, (state.soc * (1. - h2tank.η_self * Δh) - h2tank.α_soc_max) * state.Erated / Δh / h2tank.η_ch), 0.)
     # SoC dynamic
     soc_next = state.soc * (1. - h2tank.η_self * Δh) - (power_ch * h2tank.η_ch + power_dch / h2tank.η_dch) * Δh / state.Erated
     return soc_next, power_dch + power_ch
end

### Investment dynamic
function compute_investment_dynamics!(y::Int64, s::Int64, h2tank::H2Tank, decision::Union{Float64, Int64})
    h2tank.Erated[y+1,s], h2tank.soc[1,y+1,s] = compute_investment_dynamics(h2tank, (Erated = h2tank.Erated[y,s], soc = h2tank.soc[end,y,s]), decision)
end

function compute_investment_dynamics(h2tank::H2Tank, state::NamedTuple{(:Erated, :soc), Tuple{Float64, Float64}}, decision::Union{Float64, Int64})
    if decision > 1e-2
        Erated_next = decision
        soc_next = h2tank.soc_ini
    else
        Erated_next = state.Erated
        soc_next = state.soc
    end
    return Erated_next, soc_next
end


text:
#Uitlity functions

" Spectral heat flux "
function transfer_w(T1,T2,w,τ)
    return (bose_einstein(w,T1)-bose_einstein(w,T2))*τ
end

"Total heat flux obtained from simulation"
function total_transfer(T1,T2,w,τ)
  t_w   = transfer_w.(T1,T2,w,τ)
  q_tot = trapz(w,t_w)
  return q_tot
end

"Generates frequency file that serves as imput to scuff-em"
function frequency_file(freqnum,T,typespace :: String ; f1=0.01 ,f2=10.0)
    freq = zeros(Float64,freqnum)
    if typespace == "linspace"
        freq =collect(linspace(f1*wien(T)/w0,f2*wien(T)/w0,freqnum))
        dffreq = DataFrame(Freq=freq)
        writetable("frequency_N="*string(freqnum)*"_"*typespace*".dat", dffreq,header=false)
    elseif typespace == "logspace"
        wi = log10(f1*wien(T)/w0)
        wf = log10(f2*wien(T)/w0)
        freq =collect(logspace(wi,wf,freqnum))
        dffreq = DataFrame(Freq=freq)
        writetable("frequency_N="*string(freqnum)*"_"*typespace*".dat", dffreq,header=false)
    end
end


text:
@testset "enumeration" begin
    @jyro function model()
        a ~ Bernoulli(0.5)
        b ~ Bernoulli(0.5)
        c ~ Bernoulli(0.5)
        # Use this as a deterministic function.
        s ~ DiscreteNonParametric([a + b + c], [1.0])
        return s
    end

    dist = discrete_enumeration(model, :_return)
    analytic_dist = DiscreteNonParametric([0, 1, 2, 3], [0.125, 0.375, 0.375, 0.125])
    @test dist == analytic_dist
    dist_s = discrete_enumeration(model, :s)
    @test dist_s == analytic_dist
end

@testset "helper functions" begin
    D = 2
    @jyro function model()
        x ~ Normal()
        y ~ MvNormal(D, 1)
        z ~ Normal()
    end

    condition!(model, Dict(:z => 0))
    param_info, num_params = Minijyro.get_param_info(model, ())

    @test num_params == (D + 1)
    @test !haskey(param_info, :z)
    @test param_info[:x][:is_scalar]
    @test !param_info[:y][:is_scalar]
    @test param_info[:y][:num_elems] == D
    @test param_info[:y][:shape] == (D,)

    params_dict = Dict()
    new_params = [0.0, 1.0, 1.0]
    Minijyro.load_params!(params_dict, param_info, new_params)

    params_dict2 = Dict()
    Minijyro.load_params!(params_dict2, param_info, new_params)

    # Because we loop through a dict we do not know the order of keys beforehand
    # but it should always be consistent.
    @test params_dict[:x] == params_dict2[:x]
    @test params_dict[:y] == params_dict2[:y]
end


text:
module ScanDir

export scandir, scandirtree, DirEntry

import Base.Filesystem: uv_dirent_t
import Base: uv_error, _sizeof_uv_fs

module PosixFileTypes
    @enum PosixFileType Unknown File Directory Link FIFO Socket CharDev BlockDev
end
using .PosixFileTypes: PosixFileType

struct DirEntry
    name::String
    path::String
    type::PosixFileType
end

function Base.show(io::IO, ::MIME"text/plain", e::DirEntry)
    print(io, "<", lowercase(string(e.type)), " ", repr(e.path), ">")
end

_islink(e::DirEntry) = e.type == PosixFileTypes.Link
for (i,s) in enumerate((:isfile, :isdir, :islink, :isfifo, :issocket, :ischardev, :isblockdev))
    @eval Base.Filesystem.$s(e::DirEntry) = e.type == PosixFileType($i) || _islink(e) && $s(e.path)
end

filename(e::DirEntry) = e.name

# Implementation copied from Base.readdir and modified to return DirEntry's
function scandir(dir::AbstractString="."; sort=true)
    # Allocate space for uv_fs_t struct
    req = Libc.malloc(_sizeof_uv_fs)
    try
        # defined in sys.c, to call uv_fs_readdir, which sets errno on error.
        err = ccall(:uv_fs_scandir, Int32, (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cint, Ptr{Cvoid}),
                    C_NULL, req, dir, 0, C_NULL)
        err < 0 && uv_error("readdir($(repr(dir)))", err)

        # iterate the listing into entries
        entries = DirEntry[]
        ent = Ref{uv_dirent_t}()
        while Base.UV_EOF != ccall(:uv_fs_scandir_next, Cint, (Ptr{Cvoid}, Ptr{uv_dirent_t}), req, ent)
            ent_name = unsafe_string(ent[].name)
            ent_path = joinpath(dir, ent_name)
            push!(entries, DirEntry(ent_name, ent_path, PosixFileType(ent[].typ)))
        end

        # Clean up the request string
        # on newer julia versions this can be: Base.Filesystem.uv_fs_req_cleanup(req)
        ccall(:uv_fs_req_cleanup, Cvoid, (Ptr{Cvoid},), req)

        # sort entries unless opted out
        sort && sort!(entries; by=filename)

        return entries
    finally
        Libc.free(req)
    end
end

_channel_try_io(f, channel, onerror) = try
    f()
catch err
    isa(err, Base.IOError) || rethrow()
    try
        onerror(err)
    catch err2
        close(channel, err2)
    end
    return
end

# Implementation copied from Base.walkdir and modified to use scandir,
# to avoid unnecessary stat()ing
function scandirtree(root="."; topdown=true, follow_symlinks=false, onerror=throw, prune=_->false)
    function _scandirtree(chnl, root)
        isfilelike(e) = (!follow_symlinks && islink(e)) || !isdir(e)
        tryf(f, p) = _channel_try_io(()->f(p), chnl, onerror)

        content = tryf(scandir, root)
        content === nothing && return
        dirs = DirEntry[]
        files = DirEntry[]
        for entry in content
            prune(entry) && continue
            filelike = tryf(isfilelike, entry)
            filelike === nothing && return
            push!(filelike ? files : dirs, entry)
        end

        if topdown
            push!(chnl, (; root, dirs, files))
        end
        for dir in dirs
            _scandirtree(chnl, joinpath(root, dir.name))
        end
        if !topdown
            push!(chnl, (; root, dirs, files))
        end
        nothing
    end
    TreeEntry = NamedTuple{(:root, :dirs, :files), Tuple{String, Vector{DirEntry}, Vector{DirEntry}}}
    return Channel{TreeEntry}(chnl -> _scandirtree(chnl, root))
end

function walkdir(root="."; topdown=true, follow_symlinks=false, onerror=throw, prune=_->false)
    scan_channel = scandirtree(root; topdown, follow_symlinks, onerror, prune)
    WalkdirEntry = NamedTuple{(:root, :dirs, :files), Tuple{String, Vector{String}, Vector{String}}}
    return Channel{WalkdirEntry}() do channel
        _channel_try_io(channel, onerror) do
            for (root, dirs, files) in scan_channel
                push!(channel, (; root, dirs = [e.name for e in dirs], files = [e.name for e in files]))
            end
        end
    end
end

end # module


text:
using NCDatasets
using Formatting

macro pf(exps...)
    return :(print(format($(eval(exps)...))))
end

macro pfln(exps...)
    return :(println(format($(eval(exps)...))))
end

filename = ARGS[1]
output_filename = ARGS[2]

@pfln("Input file: {}", filename)
@pfln("Output file: {}", output_filename)

first_file = true

ds_o = Dataset(output_filename, "c")

cnt = 1
for line in eachline(filename)
    global first_file, cnt, missing_value

    println("Reading data: ", line)

    ds_i = Dataset(line, "r")

    if first_file

        first_file = false

        missing_value = ds_i["T_ML"].attrib["_FillValue"]
        Nx = ds_i.dim["Nx"]
        Ny = ds_i.dim["Ny"]

        defDim(ds_o, "Nx", Nx)
        defDim(ds_o, "Ny", Ny)
        defDim(ds_o, "time", Inf)

        defVar(ds_o, "T_ML", Float64, ("Nx", "Ny", "time"), fillvalue=missing_value)

        global buffer = zeros(Float64, Nx, Ny)
    end

    buffer[:, :] = nomissing( ds_i["T_ML"][:, :], missing_value)

    ds_o["T_ML"][:, :, cnt] = buffer
    close(ds_i)

    cnt += 1
end

close(ds_o)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
Sys.iswindows() && (ENV["MPLBACKEND"]="qt4agg")
using PyPlot, Test, Random, Revise, Statistics, LinearAlgebra
srcdir = dirname(pwd())*"/src"
any(srcdir .== LOAD_PATH) || push!(LOAD_PATH, srcdir)
push!(LOAD_PATH, pwd())
import GP, TransD_GP
## make options for the multichannel lengthscale GP
nminlog10λ, nmaxlog10λ = 2, 200
pnorm = 2.
Klog10λ = GP.Mat32()
λx,λy,λz = 1, 1, 1
x = 0:(0.01λx):λx-0.01λx
y = 0:(0.01λy):λy-0.01λy
z = 0:(0.01λy):λy/2-0.01λz
λlog10λ = [0.05maximum(x), 0.05maximum(y), 0.05maximum(z)]
demean = false
sdev_poslog10λ = [0.01maximum(y), 0.01maximum(x), 0.01maximum(z)]
log10bounds = [-1 -0.69; -1 -0.69; -1 0.69]
δlog10λ = 0.1
sdev_proplog10λ = [0.1, 0.1, 0.1]
xall = zeros(3,length(x)*length(y)*length(z))
for i in 1:size(xall,2)
    xid, yid, zid = Tuple(CartesianIndices((length(x),length(y), length(z)))[i])
    xall[:,i] = [x[xid]; y[yid]; z[zid]]
end
xbounds = zeros(Float64,size(xall, 1), 2)
for dim in 1:size(xall, 1)
    xbounds[dim,:] = [minimum(xall[dim,:]), maximum(xall[dim,:])]
end
## Initialize a model using these options
Random.seed!(12)
optlog10λ = TransD_GP.OptionsStat(nmin = nminlog10λ,
                        nmax = nmaxlog10λ,
                        xbounds = xbounds,
                        fbounds = log10bounds,
                        xall = xall,
                        λ = λlog10λ,
                        δ = δlog10λ,
                        demean = demean,
                        sdev_prop = sdev_proplog10λ,
                        sdev_pos = sdev_poslog10λ,
                        pnorm = pnorm,
                        quasimultid = false,
                        K = Klog10λ,
                        timesλ = 3.6
                        )
@time  log10λ = TransD_GP.init(optlog10λ)
## make options for the nonstationary GP
nmin, nmax = 2, 200
fbounds = [-2. 2]
δ = 0.1
sdev_prop = [0.1]
sdev_pos = [0.05, 0.05, 0.05]
K = GP.Mat32()
demean_ns = true
## Initialize model for the nonstationary GP
Random.seed!(13)
opt = TransD_GP.OptionsNonstat(optlog10λ,
                        nmin = nmin,
                        nmax = nmax,
                        fbounds = fbounds,
                        δ = δ,
                        demean = demean_ns,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        K = K
                        )
@time m = TransD_GP.init(opt, log10λ)
## timing for s birth in s model for ns mode
log10λ = TransD_GP.init(optlog10λ)
m = TransD_GP.init(opt, log10λ)
for i = 1:98
    TransD_GP.birth!(log10λ, optlog10λ, m, opt)
end
for i = 1:98
    TransD_GP.birth!(m, opt, log10λ)
end
NTIMES = 1
ntimes = 150
T = zeros(NTIMES)
for I = 1:NTIMES
    T[I] = time()
    for i = 1:ntimes
        TransD_GP.death!(log10λ, optlog10λ, m, opt)
        TransD_GP.birth!(log10λ, optlog10λ, m, opt)
    end
    T[I] = (time() - T[I])/2ntimes
end
@info "time for $ntimes birth/death is $(mean(T)) +- $(std(T)/sqrt(NTIMES))"


text:
using Test, QBase

@testset "./src/types/probabilities.jl" begin

@testset "is_probability_distribution()" begin
    @testset "valid cases" begin
        @test is_probability_distribution([1])
        @test is_probability_distribution([0,0,0,1,0])
        @test is_probability_distribution(1/100*fill(1,100))
        @test is_probability_distribution([0.2,0,0,0.1,0.1,0.1,0.3,0.2,0])
    end

    @testset "invalid cases" begin
        @test !is_probability_distribution([0])
        @test !is_probability_distribution([0.2,0.1,0.11,0.4])
        @test !is_probability_distribution([-0.1,1,0.1])
        @test !is_probability_distribution([-0.1,-0.2])
    end

    @testset "atol" begin
        @test !is_probability_distribution([0.3,0.3,0.39])
        @test is_probability_distribution([0.3,0.3,0.39], atol=2e-2)
        @test !is_probability_distribution([0.3,0.3,0.41,-0.01])
        @test is_probability_distribution([0.3,0.3,0.41,-0.01], atol=2e-2)
    end

    @testset "distribution types" begin
        @test is_probability_distribution(Probabilities([1,0,0]))
        @test is_probability_distribution(JointProbabilities([1 0 0;0 0 0;0 0 0]))
    end
end

@testset "is_conditional_distribution()" begin
    @testset "valid cases" begin
        @test is_conditional_distribution([1 0 0;0 1 0;0 0 1])
        @test is_conditional_distribution([1/3 1/2 1/8;1/3 1/4 0;1/3 1/4 7/8])
    end

    @testset "invalid cases" begin
        @test !is_conditional_distribution([1 1 1;1 0 0;1 0 0])
        @test !is_conditional_distribution([1 0 0;0 1 0;0 0 0])
    end

    @testset "atol" begin
        @test !is_conditional_distribution([1 0 0;0 1 0;0 0 1.1])
        @test is_conditional_distribution([1 0 0;0 1 0;0 0 1.1], atol=0.11)
    end

    @test is_conditional_distribution(Conditionals([1 0 0;0 1 0;0 0 1]))
end

@testset "Probabilities()" begin
    @testset "valid cases" begin
        probs = Probabilities([0.5,0.5])
        @test probs isa Probabilities{Float64}
        @test probs isa ProbabilityDistribution
        @test probs == [0.5,0.5]

        probs_int = Probabilities([1,0])
        @test probs_int isa Probabilities{Int64}
        @test probs_int isa ProbabilityDistribution
        @test probs_int == [1,0]

        probs_atol = Probabilities([1,0.1],atol=.11)
        @test probs_atol isa Probabilities{Float64}
        @test probs_atol isa ProbabilityDistribution
        @test probs_atol == [1,0.1]
    end

    @testset "invalid cases" begin
        @test_throws DomainError Probabilities([0.5,0.4,0.2])
        @test_throws DomainError Probabilities([0.5,-1,1.5])
    end
end

@testset "Conditionals()" begin
    @testset "valid cases" begin
        conditionals = Conditionals([1/3 1/2 1/8;1/3 1/4 0;1/3 1/4 7/8])
        @test conditionals == [1/3 1/2 1/8;1/3 1/4 0;1/3 1/4 7/8]
        @test conditionals isa Conditionals{Float64}
        @test conditionals isa ConditionalDistribution
    end

    @testset "inexact edge case" begin
        @test_throws DomainError Conditionals([
            0.666667 0.166667 0.166667; 0.166667 0.666667 0.166667; 0.166667 0.166667 0.666667
        ])

        conditionals = Conditionals([
            0.666667 0.166667 0.166667; 0.166667 0.666667 0.166667; 0.166667 0.166667 0.666667
        ], atol=1e-5)
        @test conditionals == [0.666667 0.166667 0.166667; 0.166667 0.666667 0.166667; 0.166667 0.166667 0.666667]
    end

    @testset "invalid cases" begin
        @test_throws DomainError Conditionals([0.5 -0.5;0.5 1.5])
    end
end

@testset "JointProbabilities()" begin
    @testset "matrix constructor" begin
        joints = JointProbabilities([1 0 0;0 0 0;0 0 0])
        @test joints isa JointProbabilityDistribution
        @test joints isa JointProbabilities{Int64}
        @test joints == [1 0 0;0 0 0;0 0 0]
    end

    @testset "errors" begin
        @test_throws DomainError JointProbabilities([1 1 0;0 0 0;0 0 0])
    end

    @testset "atol" begin
        @test_throws DomainError JointProbabilities([1.1 0 0;0 0 0;0 0 -0.1])
        @test JointProbabilities([1.1 0 0;0 0 0;0 0 -0.1], atol=0.11) isa JointProbabilityDistribution

        @test_throws DomainError JointProbabilities([1,0,0],[1,0,0.1])
        @test JointProbabilities([1,0,0],[1,0,0.1], atol=0.11) isa JointProbabilityDistribution
    end

    @testset "priors and conditionals" begin
        priors_vec = [0.3,0.6,0.1]
        priors = Probabilities(priors_vec)

        conditionals_mat = [0.3 0.5 0.4;0.7 0.5 0.6]
        conditionals = Conditionals(conditionals_mat)

        joint_probs = JointProbabilities(priors,conditionals)
        @test joint_probs ≈ [0.09 0.30 0.04;0.21 0.30 0.06]
        @test joint_probs isa JointProbabilityDistribution
        @test joint_probs isa JointProbabilities{Float64}

        joint_probs = JointProbabilities(priors_vec,conditionals_mat)
        @test joint_probs ≈ [0.09 0.30 0.04;0.21 0.30 0.06]
        @test joint_probs isa JointProbabilityDistribution
        @test joint_probs isa JointProbabilities{Float64}

        joint_probs = JointProbabilities(conditionals_mat, priors_vec)
        @test joint_probs ≈ [0.09 0.30 0.04;0.21 0.30 0.06]
        @test joint_probs isa JointProbabilityDistribution
        @test joint_probs isa JointProbabilities{Float64}

        joint_probs = JointProbabilities(priors_vec,conditionals)
        @test joint_probs ≈ [0.09 0.30 0.04;0.21 0.30 0.06]
        @test joint_probs isa JointProbabilityDistribution
        @test joint_probs isa JointProbabilities{Float64}
    end

    @testset "priors and priors" begin
        priors1_vec = [0.9,0,0.1]
        priors1 = Probabilities(priors1_vec)

        priors2_vec = [0,0,0.5,0.5]
        priors2 = Probabilities(priors2_vec)

        joint_probs = JointProbabilities(priors1, priors2)
        @test joint_probs == [0 0 0;0 0 0;0.45 0 0.05;0.45 0 0.05]
        @test joint_probs isa JointProbabilities{Float64}

        joint_probs = JointProbabilities(priors2, priors1)
        @test joint_probs == [0 0 0.45 0.45;0 0 0 0;0 0 0.05 0.05]

        joint_probs = JointProbabilities(priors1_vec, priors2_vec)
        @test joint_probs == [0 0 0;0 0 0;0.45 0 0.05;0.45 0 0.05]
    end
end

@testset "outcome_probabilities()" begin
    @testset "multiplication" begin
        priors = Probabilities([0.3,0.3,0.4])
        conditionals = Conditionals([1 0 0;0 1 0;0 0 1])

        probs = conditionals*priors
        @test probs isa Probabilities{Float64}
        @test probs isa ProbabilityDistribution
        @test probs == [0.3,0.3,0.4]
    end

    @testset "ordering of inputs" begin
        priors_vec = [0.3,0.6,0.1]
        priors = Probabilities(priors_vec)

        conditionals_mat = [0.3 0.5 0.4;0.7 0.5 0.6]
        conditionals = Conditionals(conditionals_mat)

        out_probs = outcome_probabilities(priors,conditionals)
        @test out_probs ≈ [0.43,0.57]
        @test out_probs isa Probabilities{Float64}

        out_probs = outcome_probabilities(conditionals,priors)
        @test out_probs ≈ [0.43,0.57]
        @test out_probs isa Probabilities{Float64}

        out_probs = outcome_probabilities(conditionals_mat,priors_vec)
        @test out_probs ≈ [0.43,0.57]
        @test out_probs isa Probabilities{Float64}

        out_probs = outcome_probabilities(priors_vec,conditionals_mat)
        @test out_probs ≈ [0.43,0.57]
        @test out_probs isa Probabilities{Float64}
    end
end

end


text:
#!/usr/bin/env julia
# Comprehension and Generators
# https://teaching.healthtech.dtu.dk/22110/index.php/Comprehension_and_Generators

# Ex 1
matrix1 = readdlm("data/matrix1", '\t', Int, '\n')
matrix2 = readdlm("data/matrix2", '\t', Int, '\n')

dot = matrix1*matrix2
matrix_product(a,b) = [sum([a[i,r]*b[r,j] for r = 1:size(a,2)])
    for i = 1:size(a,1), j = 1:size(b,2)]

matrix_product(matrix1,matrix2)

# Ex 2
function _read_experiment(f, acc_match)
    status = Float64[]
    for line = eachline(f)
        line_split = split(line, "\t")
        acc = string(line_split[2])
        if line_split[1] == "COL_CLASSES"
            status = [parse(Int64, x) for x in line_split[4:end]]
        elseif acc == acc_match
            nums = [parse(Float64, x) for x in line_split[4:end]]
            return nums[status .== 0], nums[status .== 1]
        end
    end
    return Float64[], Float64[]
end

function read_experiment(acc_search, outfilename)
    out = open(outfilename, "w")
    cancer, control = open(f->_read_experiment(f, acc_search), "data/dna-array.dat")
    if length(cancer) == 0 && length(control) == 0
        println("Accession number not found.")
        close(out)
        return nothing
    end
    min_len = min(length(cancer), length(control))
    for i = 1:min_len
        println(out, cancer[i], "\t", control[i])
    end
    if length(cancer) > length(control)
        for i = (min_len+1):length(cancer)
            println(out, cancer[i])
        end
    else
        for i = (min_len+1):length(control)
            println(out, "\t", cancer[i])
        end
    end
    close(out)
end

println("Enter accession number:")
read_experiment(chomp(readline()), "data/dna-array-subset.dat")


text:
@testset "NfAbs" begin
  include("NfAbs/Conjugates.jl")
  include("NfAbs/NonSimple.jl")
end


text:
export
    KnotPoint,
    Traj,
    state,
    control

abstract type AbstractKnotPoint{T,N,M,NM} end

""" $(TYPEDEF)
Stores critical information corresponding to each knot point in the trajectory optimization
problem, including the state and control values, as well as the time and time step length.

# Getters
Use the following methods to access values from a `KnotPoint`:
```julia
x  = state(z::KnotPoint)    # returns the n-dimensional state as a SVector
u  = control(z::KnotPoint)  # returns the m-dimensional control vector as a SVector
t  = z.t                    # current time
dt = z.dt                   # time step length
```

# Setters
Use the following methods to set values in a `KnotPoint`:
```julia
set_state!(z::KnotPoint, x)
set_control!(z::KnotPoint, u)
z.t = t
z.dt = dt
```

# Constructors
```julia
KnotPoint(x, u, dt, t=0.0)
KnotPoint(x, m, t=0.0)  # for terminal knot point
```

Use `is_terminal(z::KnotPoint)` to determine if a `KnotPoint` is a terminal knot point (e.g.
has no time step length and z.t == tf).
"""
mutable struct KnotPoint{T,N,M,NM} <: AbstractKnotPoint{T,N,M,NM}
    z::SVector{NM,T}
    _x::SVector{N,Int}
    _u::SVector{M,Int}
    dt::T # time step
    t::T  # total time
end

function KnotPoint(x::AbstractVector, u::AbstractVector, dt::Float64, t=0.0)
    n = length(x)
    m = length(u)
    xinds = ones(Bool, n+m)
    xinds[n+1:end] .= 0
    _x = SVector{n}(1:n)
    _u = SVector{m}(n .+ (1:m))
    z = SVector{n+m}([x;u])
    KnotPoint(z, _x, _u, dt, t)
end

# Constructor for terminal time step
function KnotPoint(x::AbstractVector, m::Int, t=0.0)
    u = zeros(m)
    KnotPoint(x, u, 0., t)
end

@inline state(z::AbstractKnotPoint) = z.z[z._x]
@inline control(z::AbstractKnotPoint) = z.z[z._u]
@inline is_terminal(z::AbstractKnotPoint) = z.dt == 0

const Traj = AbstractVector{<:KnotPoint}
traj_size(Z::Vector{<:KnotPoint{T,N,M}}) where {T,N,M} = N,M,length(Z)

function Base.copy(Z::Vector{KnotPoint{T,N,M,NM}}) where {T,N,M,NM}
    [KnotPoint((@SVector ones(NM)) .* z.z, z._x, z._u, z.dt, z.t) for z in Z]
end

function Traj(n::Int, m::Int, dt::AbstractFloat, N::Int, equal=false)
    x = NaN*@SVector ones(n)
    u = @SVector zeros(m)
    Traj(x,u,dt,N,equal)
end

function Traj(x::SVector, u::SVector, dt::AbstractFloat, N::Int, equal=false)
    equal ? uN = N : uN = N-1
    Z = [KnotPoint(x,u,dt,(k-1)*dt) for k = 1:uN]
    if !equal
        m = length(u)
        push!(Z, KnotPoint(x,m,(N-1)*dt))
    end
    return Z
end

function Traj(X::Vector, U::Vector, dt::Vector, t=cumsum(dt) .- dt[1])
    Z = [KnotPoint(X[k], U[k], dt[k], t[k]) for k = 1:length(U)]
    if length(U) == length(X)-1
        push!(Z, KnotPoint(X[end],length(U[1]),t[end]))
    end
    return Z
end

@inline states(Z::Traj) = state.(Z)
@inline controls(Z::Traj) = control.(Z[1:end-1])

states(Z::Traj, i::Int) = [state(z)[i] for z in Z]

set_state!(z::KnotPoint, x) = z.z = [x; control(z)]
set_control!(z::KnotPoint, u) = z.z = [state(z); u]

function set_states!(Z::Traj, X)
    for k in eachindex(Z)
        Z[k].z = [X[k]; control(Z[k])]
    end
end

function set_controls!(Z::Traj, U)
    for k in 1:length(Z)-1
        Z[k].z = [state(Z[k]); U[k]]
    end
end

function set_controls!(Z::Traj, u::SVector)
    for k in 1:length(Z)-1
        Z[k].z = [state(Z[k]); u]
    end
end

function set_times!(Z::Traj, ts)
    for k in eachindex(ts)
        Z[k].t = ts[k]
    end
end

function get_times(Z::Traj)
    [z.t for z in Z]
end

function shift_fill!(Z::Traj)
    N = length(Z)
    for k in eachindex(Z)
        Z[k].t += Z[k].dt
        if k < N
            Z[k].z = Z[k+1].z
        else
            Z[k].t += Z[k-1].dt
        end
    end
end

struct StaticKnotPoint{T,N,M,NM} <: AbstractKnotPoint{T,N,M,NM}
    z::SVector{NM,T}
    _x::SVector{N,Int}
    _u::SVector{M,Int}
    dt::T # time step
    t::T  # total time
end


text:
# -----------------------------------------------------------------------
# Radiative forcing from carbon dioxide (Etminan et al. 2016 equations).
# -----------------------------------------------------------------------

@defcomp rf_co2_etminan begin

    a₁           = Parameter()             # Carbon dioxide forcing constant.
    b₁           = Parameter()             # Carbon dioxide forcing constant.
    c₁           = Parameter()             # Carbon dioxide forcing constant.
    CO₂_0        = Parameter()             # Initial (pre-industrial) carbon dioxide concentration (ppm).
    N₂O_0        = Parameter()             # Initial (pre-industrial) nitrous oxide concentration (ppb).
    rf_scale_CO₂ = Parameter()             # Radiative forcing scaling term (to keep forcing from 2x CO₂ consistent).
    N₂O          = Parameter(index=[time]) # Atmospheric nitrous oxide concentration (ppb).
    CO₂          = Parameter(index=[time]) # Atmospheric carbon dioxide concentration (ppm).

    rf_CO₂       = Variable(index=[time])  # Forcing from atmospheric carbon dioxide concentrations (Wm⁻²).

    function run_timestep(p, v, d, t)

        if is_first(t)
            # Set initial radiative forcing to 0.
            v.rf_CO₂[t] = 0.0

        else
            # Create temporary variables.
            CO₂_diff = p.CO₂[t]-p.CO₂_0
            N̄ = 0.5 * (p.N₂O[t] + p.N₂O_0)

            # Calculate carbon dioxide radiative forcing.
            v.rf_CO₂[t] = ((p.a₁*CO₂_diff^2 + p.b₁*abs(CO₂_diff) + p.c₁*N̄ + 5.36) * log(p.CO₂[t] / p.CO₂_0)) * p.rf_scale_CO₂
        end
    end
end


text:
using Test

include("framework.jl")

# TODO test integrate with ts[1] both zero and nonzero

"""
Consistency checks for operators and matrices:

julia> TT ≈ op2mat(T)

julia> VV ≈ diagm(0=>V[:])

julia> UU ≈ diagm(0=>U(ψ)[:])

julia> JJ ≈ op2mat(J)

function op2mat(f)
    M = similar(z,length(z),length(z))
    u = similar(z)
    for j = eachindex(u)
        u .= 0
        u[j] = 1
        M[:,j] = f(u)[:]
    end
    M
end

Check zero mode for BdGmat
"""

# Check that relaxed orbiting vortices and lattices have the right
# winding number around the trap edge.

# Check consistency when hbm, C and grid are changed invariantly
# Given one set, fuzerate identical ones

@testset "Derivatives of cloud" begin
    # Cloud is cos, so the second derivative flips the sign
    for d = fuzerate(Discretisation)
        l = d.h*(d.n+1)
        u = cloud(d)
#         @test D(u, 1, d, 2) ≈ -(π/l)^2*u
#         @test D(D(u, 1, d), 1, d) ≈ D(u, 1, d, 2)
#         @test D(u, 2, d, 2) ≈ -(π/l)^2*u
#         @test D(D(u, 2, d), 2, d) ≈ D(u, 2, d, 2)
    end
end

# TODO extend fuzerate with a "this is no less converged than the last time the test ran" macro and an order "this test data should converge better than that does"

@testset "Compare operators to explicit matrices" begin
    d = FDDiscretisation{2}(10, 10/9)
    s = Superfluid{2}(500, (x,y) -> (x^2+y^2)/2, hbm=rand())
    Ω = rand()
    # TODO fuzerate wave functions
    ψ = cloud(d)
    es = [reshape(e, d.n, d.n)
        for e = eachcol(Matrix{Complex{Float64}}(I, d.n^2, d.n^2))]

    function compare(a, f)
        A = Superfluids.operators(s,d,a) |> only
        M = Superfluids.matrices(s,d,Ω,ψ,a) |> only
        @test Superfluids.expand(q->f(A,q),es) ≈ M
    end

    compare(:U, (U,q)->U(q,abs2.(ψ)))
    compare(:L, (L,q)->L(q,abs2.(ψ);Ω))
    compare(:H, (H,q)->H(q,abs2.(ψ);Ω))
    compare(:T, (T,q)->T(q))
    compare(:V, (V,q)->V(q))
    compare(:J, (J,q)->J(q))
end

@testset "Compare Hartree modes to matrix eigenvectors" begin
    hbm = rand()
    Ω = 0.01*hbm	# break degeneracy
    s = Superfluid{2}(500, (x,y)->(x^2+y^2)/2; hbm)
    d = FDDiscretisation{2}(40, 15/39)
    L, J = operators(s,d,:L,:J)
    ψ = steady_state(s,d;Ω)
    μ = dot(L(ψ;Ω), ψ) |> real
    M = matrices(s,d,Ω,ψ,:L) |> only
    ew, ev = eigen(M)
    ws, us = hartree_modes(s,d,ψ,10,Ω)
    @test norm(imag(ws), Inf) < 1e-10
    ws = real(ws)
    @test ws[1] ≈ μ
    @test ws ≈ ew[eachindex(ws)]
    for (j, u) = pairs(us)
        @test abs(dot(u, ev[:,j])) ≈ 1
    end
    jj = [dot(u,J(u)) for u in us]
    @test norm(imag(jj), Inf) < 1e-10
    @test norm(round.(jj) - jj, Inf) < 0.1
end

@testset "Compare BdG modes with matrix" begin
    # leave hbm = 1 to keep core smaller than trap
    Ω = 0.01
    s = Superfluid{2}(500, (x,y)->(x^2+y^2)/2, hbm=0.5)
    d = FDDiscretisation{2}(40, 15/39)
    L, H, J = Superfluids.operators(s,d,:L,:H,:J)
    ψ = steady_state(s,d;Ω)
    @test norm(ψ) ≈ 1

    # Check that this problem is solvable
    M = Superfluids.BdGmatrix(s,d,Ω,ψ)
    ew, ev = eigs(M, nev=20, which=:SM)
    @test norm(imag(ew), Inf) < 1e-3
    ew = real(ew)
    @test ew[1] ≈ 0 atol=1e-6
    for j = 2:2:length(ew)
        @test ew[j-1] ≈ -ew[j] rtol=1e-3
    end
    uu = [reshape(u[1:d.n^2],d.n,d.n) for u in eachcol(ev)]
    vv = [reshape(v[d.n^2+1:end],d.n,d.n) for v in eachcol(ev)]
    @test norm(uu[1]) ≈ 1/√2
    @test norm(vv[1]) ≈ 1/√2
    @test abs(dot(uu[1],ψ)) ≈ 1/√2 rtol=1e-3
    @test abs(dot(vv[1],ψ)) ≈ 1/√2 rtol=1e-3
    ix = ew .≥ 0
    ew = ew[ix]
    uu = uu[ix]
    vv = vv[ix]

    # Check expansions
    hw, hv = hartree_modes(s,d,ψ,100,Ω)
    @test norm(imag(hw), Inf) < 1e-6
    hw = real(hw)

    uns = [norm(u)^2 for u in uu]
    ups = [sum(abs2, dot(h,u) for h in hv) for u in uu]
    @test uns ≈ ups rtol=0.01
    vns = [norm(v)^2 for v in vv]
    vps = [sum(abs2, dot(conj(h),v) for h in hv) for v in vv]
    @test vns ≈ vps rtol=0.01

    ws, us, vs = bdg_modes(s, d, ψ, Ω, 100)
    @test norm(imag(ws), Inf) < 1e-3
    ws = real(ws)
    @test ew ≈ ws[eachindex(ew)] rtol=0.05
    for j = eachindex(ew)
        u, v = us[j], vs[j]
        u0, v0 = uu[j], vv[j]
        @test abs(dot([u0[:]; v0[:]], [u[:]; v[:]])) ≈ 1 rtol=0.03
    end
end

@testset "SHO ground state expectation values" begin
    s = Superfluid{2}(0, (x,y) -> (x^2+y^2)/2)
    for d = fuzerate(Discretisation)
        T, V, U, J = operators(s,d,:T,:V,:U,:J)
        z = argand(d)
        w = similar(z)
        w .= normalize(@. exp(-abs2(z)/2))
        @test dot(w,T(w)) ≈ 0.5 atol=0.02
        @test dot(w,V(w)) ≈ 0.5 atol=0.02
        @test dot(w,J(w)) ≈ 0.0 atol=1e-10
        @test dot(w,U(w)) ≈ 0.0
    end
end

@testset "SHO ground state is an eigenstate" begin
    s = Superfluid{2}(0, (x,y) -> (x^2+y^2)/2)
    for d = fuzerate(Discretisation)
        L = operators(s,d,:L) |> only
        z = argand(d)
        w = similar(z)
        w .= normalize(@. exp(-abs2(z)/2))
        @test dot(w,L(w)) ≈ 1.0 atol=0.02
        @test norm(L(w)-w) < 0.05
    end
end

@testset "chemical potential" begin
    # TODO check this analytically
    s = Superfluid{2}(500, (x,y) -> (x^2+y^2)/2)
    for d = fuzerate(Discretisation)
        L = operators(s,d,:L) |> only
        q = steady_state(s,d)
        @test dot(q, L(q)) ≈ 12.678 atol=0.02
    end
end

@testset "Hartree BdG matrix compared to explicit form"

end

@testset "Order parameter is zero mode" begin

end

@testset "BdG matrix annihilates zero modes"  begin
    # ew, us = hartree_modes(s, d, ψ, nev, Ω);
    # Umat = hcat([u[:] for u in us]...);
    # cs = Umat'*q[:];
    # ds = conj(cs)
    # M = Superfluids.BdGmatrix(s, d, Ω, ψ, us)
    # A = M[1:10, 1:10];
    # B = M[1:10, 11:20];
    # C = M[11:20, 1:10];
    # D = M[11:20, 11:20];
    # Mu = Umat*(A*cs+B*ds);
    # Mu = reshape(Mu, d.n, d.n);
    # Mv = conj(Umat)*(C*cs+D*ds);
    # Mv = reshape(Mv, d.n, d.n);
end


text:
# Lewis, A. L. (2001)
# A simple option formula for general jump-diffusion and other exponential Lévy processes
# Available at SSRN 282110

struct Lewis <: AbstractPricer
  x::Array{Float64,1}
  w::Array{Float64,1}
  function Lewis(n=512)
    x, w = laguerre(n, 0.0)
    idx = w .!= 0.0
    new(x[idx], w[idx])
  end
end

_lewis_integrand(m, x0, x, k, r, q, T, u) =
  real(exp(u + logmgf(m, x0, T, 0.5 + u * im) + im * u * (x - k + (r - q) * T))) / (u^2 + 0.25)

function integral(m, p::Lewis, x0, x, k, r, q, T)
  z = 0.0
  @simd for i in 1:length(p.x)
    @inbounds z += p.w[i] * _lewis_integrand(m, x0, x, k, r, q, T, p.x[i])
  end
  z
end

function price(model, pricer::Lewis, x0, strike, rate, time; yield=0.0)
  stock = x0[1]
  stock * exp(-yield * time) - 1 / pi * sqrt(stock * strike) * exp(-0.5 * (rate + yield) * time) * integral(model, pricer, x0, log(stock), log(strike), rate, yield, time)
end


text:
using Logging
using VectorLogging
using NahaJuliaLib

@testset "Recognize empty Sequence" begin
    @test recognize(Sequence(), "abcd") == (true, [], 1)
    @test recognize(CharacterLiteral('a'), "ab") == (true, 'a', 2)
    @test recognize(Alternatives(CharacterLiteral('a')), "bb") == (false, nothing, 1)
    @test recognize(Alternatives(CharacterLiteral('a')), "ab") == (true, 'a', 2)
    @test recognize(Alternatives(Empty(), CharacterLiteral('a')),
                    "bb") == (true, nothing, 1)
    @test recognize(Alternatives(Empty(),
                                 CharacterLiteral('a')), "ab") == (true, 'a', 2)
end

@testset "Test hand coded BNF grammar" begin
    @test recognize(BNFRef(:BootstrapBNFGrammar, "<character1>"), "abcd") ==
        (true, 'a', 2)
    let
        matched, v, i = recognize(BNFRef(:BootstrapBNFGrammar, "<literal>"),
                                  "'abcd'")
        @test matched == true
        @test i == 7
        @test v isa StringLiteral
        @test v.str == "abcd"
    end
    #=
    logger = VectorLogger()
    @eval(AnotherParser, trace_recognize = true)
    @test AnotherParser.trace_recognize == true
    with_logger(logger) do
        try
            r = recognize(BNFRef(:BootstrapBNFGrammar, "<literal>"),
                          "'abcd'")
            @test r == (true, "abcd", 7)
        finally
            @test length(logger.log) > 0
            show_trace(analyze_traces(logger)[1])
        end
    end
    =#
end


text:
__precompile__()

"""
This module implements computation of kinship and other identity coefficients.
The kinship matrix is computed recursively by the standard formula.
The delta7 matrix is derived from the kinship matrix for non-inbred pedigrees.
The remaining identity coefficients of Jacquard are computed by simulation.
"""
module MendelKinship
#
# Required OpenMendel packages and modules.
#
using MendelBase
using SnpArrays
#
# Required external modules.
#
using CSV
using DataFrames
using PlotlyJS
using StatsBase

export Kinship

"""
  Kinship(control_file = ""; args...)

This is the wrapper function for the Kinship analysis option. Keyword options
should be specified in the control file.
"""
function Kinship(control_file = ""; args...)
  #
  # Print the logo. Store the initial directory.
  #
  print(" \n \n")
  println("     Welcome to OpenMendel's")
  println("     Kinship analysis option")
  print(" \n \n")
  println("Reading the data.\n")
  initial_directory = pwd()
  #
  # The user specifies the analysis to perform via a set of keywords.
  # Start the keywords at their default values.
  #
  keyword = set_keyword_defaults!(Dict{AbstractString, Any}())
  #
  # Keywords unique to this analysis should be first defined here
  # by setting their default values using the format:
  # keyword["some_keyword_name"] = default_value
  #
  keyword["kinship_output_file"] = ""
  keyword["repetitions"] = 1
  keyword["xlinked_analysis"] = false
  keyword["compare_kinships"] = false
  keyword["maf_threshold"] = 0.01
  keyword["grm_method"] = :MoM # more robust & less sensitive to rare SNPs;
                               # alternatives are :GRM and :Robust
  keyword["deviant_pairs"] = 0
  keyword["kinship_plot"] = ""
  keyword["z_score_plot"] = ""
  keyword["full_pedigree_file"] = ""
  #
  # Process the run-time user-specified keywords that will control the analysis.
  # This will also initialize the random number generator.
  #
  process_keywords!(keyword, control_file, args)
  if keyword["kinship_output_file"] == ""
    keyword["kinship_output_file"] = keyword["output_file"]
  end
  #
  # Check that the correct analysis option was specified.
  #
  lc_analysis_option = lowercase(keyword["analysis_option"])
  if (lc_analysis_option != "" &&
      lc_analysis_option != "kinship")
     throw(ArgumentError(
       "An incorrect analysis option was specified.\n \n"))
  end
  keyword["analysis_option"] = "Kinship"
  #
  # Read the genetic data from the external files named in the keywords.
  #
  (pedigree, person, nuclear_family, locus, snpdata,
    locus_frame, phenotype_frame, person_frame, snp_definition_frame) =
    read_external_data_files(keyword)
  #
  # Check if SNP data were read and names are unique.
  #
  if keyword["compare_kinships"]
    if snpdata.snps == 0
      throw(ArgumentError(
        "\n\nERROR: This analysis requires data from SNP files!\n"))
    elseif snpdata.snps != length(unique(snpdata.snpid))
      throw(ArgumentError(
        " \n\nERROR: There are non-unique names of SNPs.\n"))
    elseif person.people != length(unique(person.name))
      throw(ArgumentError(
        " \n\nERROR: There are non-unique names of individuals.\n"))
    end
  elseif snpdata.snps != 0
      throw(ArgumentError(
        " \n\nERROR: This analysis does not use data from SNP files!\n"))
  end
  #
  # Execute the specified analysis.
  #
  println(" \nAnalyzing the data.\n")
  execution_error = false
  if keyword["compare_kinships"]
    kinship_frame = compare_kinships(pedigree, person, snpdata,
                                     person_frame, keyword)
  else
    kinship_frame = theoretical_coefficients(pedigree, person, keyword)
  end
  # display(kinship_frame)

  if execution_error
    println(" \n \nERROR: Mendel terminated prematurely!\n")
  else
    println(" \n \nMendel's analysis is finished.\n")
  end
  #
  # Finish up by closing, and thus flushing, any output files.
  # Return to the initial directory.
  #
  close(keyword["output_unit"])
  cd(initial_directory)
  return nothing
end # function Kinship

"""
    compare_kinships(genotyped_pedigree, genotyped_person, snpdata,person_frame, keyword)

This function orchestrates the comparison of theoretical
and empirical kinships including Fisher z-scores.
The results are placed in a dataframe.
"""
function compare_kinships(genotyped_pedigree::Pedigree,
  genotyped_person::Person, snpdata::SnpDataStruct,
  person_frame::DataFrame, keyword::Dict{AbstractString, Any})
  #
  # Initialize constants
  #
  genotyped_people = genotyped_person.people
  deviant_pairs = keyword["deviant_pairs"]
  if deviant_pairs == 0
    deviant_pairs = genotyped_people^2
  end
  #
  # Record the order of the individuals in the user data upon entry,
  # before MendelBase permutes them.
  #
  entryorder = convert(Vector{Int64}, person_frame[:, :EntryOrder])
  #
  # If the full pedigree contains people without SNP information,
  # initialize the full pedigree and person including people without SNPs.
  # The boolean vector `contains_snp` (which is not used)
  # records if there are SNPs: contains_snp[i] means person i has SNP data.
  #
  if keyword["full_pedigree_file"] != ""
    full_pedigree, full_person, contains_snp =
      compute_full_person_and_pedigree(keyword, genotyped_person)
  else
    full_pedigree = genotyped_pedigree
    full_person = genotyped_person
  end
  #
  # Initialize more constants.
  #
  full_pedigrees = full_pedigree.pedigrees
  xlinked = keyword["xlinked_analysis"]
  maf_threshold = keyword["maf_threshold"]
  grm_method = Symbol(keyword["grm_method"])
  #
  # People are reordered, and some may not have SNP info.
  # Compute identification maps to match positions of IDs within two vectors.
  # Recall that snpdata.personid has been permuted; thus, we inverse permute.
  #
  personid_in_snpdata = copy(snpdata.personid)
  invpermute!(personid_in_snpdata, entryorder)
  (name_to_id, id_to_name) = correspond(full_person.name, personid_in_snpdata)
  #
  # Compute the genetic relationship matrix and apply Fisher's transform.
  #
  GRM = grm(snpdata.snpmatrix, method=grm_method, minmaf=maf_threshold)
  clamp!(GRM, -0.99999, 0.99999)
  GRM .= atanh.(GRM)
  #
  # For each pedigrees, compute empiric kinship coefficients
  # and subtract theoretical kinship coefficients.
  #
  kinship = Vector{Matrix{Float64}}(undef, full_pedigrees)
  for ped = 1:full_pedigrees
    kinship[ped] = kinship_matrix(full_pedigree, full_person, ped, xlinked)
    q = full_pedigree.pedstart[ped] - 1
    for j = 1:full_pedigree.individuals[ped]
      jj = name_to_id[j + q]
      jj == 0 && continue
      for i = j:full_pedigree.individuals[ped]
        ii = name_to_id[i + q]
        ii == 0 && continue
        GRM[ii, jj] = GRM[ii, jj] - atanh(kinship[ped][i, j])
        GRM[jj, ii] = GRM[ii, jj]
      end
    end
  end
  #
  # Find the indices of the most deviant pairs.
  #
  p = partialsortperm(vec(GRM), 1:deviant_pairs, by = abs, rev = true)
  genotyped_indices = CartesianIndices((genotyped_people, genotyped_people))[p]
  #
  # Enter the most deviant pairs in a dataframe.
  #
  kinship_frame = DataFrame(Pedigree1 = String[], Pedigree2 = String[],
    Person1 = String[], Person2 = String[], theoretical_kinship= Float64[],
    empiric_kinship=Float64[])
  r = 0.0
  for k = 1:deviant_pairs
    #
    # Find the index of the kth largest deviation in the GRM matrix.
    #
    (ii, jj) = (genotyped_indices[k][1], genotyped_indices[k][2])
    #
    # Map this snpid index back to the individual's name.
    #
    (i, j) = (id_to_name[ii], id_to_name[jj])
    j > i && continue
    (pedi, pedj) = (full_person.pedigree[i], full_person.pedigree[j])
    if pedi == pedj
      q = full_pedigree.pedstart[pedi] - 1
      r = GRM[ii, jj] + atanh(kinship[pedi][i - q, j - q])
      push!(kinship_frame,
        [full_pedigree.name[pedi], full_pedigree.name[pedj],
         full_person.name[i], full_person.name[j],
         kinship[pedi][i - q, j - q], tanh(r)])
    else
      push!(kinship_frame,
        [full_pedigree.name[pedi], full_pedigree.name[pedj],
         full_person.name[i], full_person.name[j],
         0.0, tanh(GRM[ii, jj])])
    end
  end
  #
  # Compute the Fisher z-scores.
  #
  kinship_frame = compute_fishers_z(kinship_frame)
  #
  # Create and save plots based on user's requests.
  #
  plot_pair = make_plot_pairs(kinship_frame)
  if keyword["kinship_plot"] != ""
    my_compare_plot = make_compare_plot(kinship_frame, plot_pair)
    # PlotlyBase.savefig(my_compare_plot, keyword["kinship_plot"] * ".pdf")
    PlotlyJS.savefig(my_compare_plot.plot, keyword["kinship_plot"] * ".html")
    println("Kinship plot saved.")
  end
  if keyword["z_score_plot"] != ""
    my_fisher_plot = plot_fisher_z(kinship_frame, plot_pair)
    # PlotlyBase.savefig(my_fisher_plot, keyword["z_score_plot"] * ".pdf")
    PlotlyJS.savefig(my_fisher_plot.plot, keyword["z_score_plot"] * ".html")
    println("Fisher's plot saved.")
  end
  #
  # Write the combined theoretical and empirical kinship frame to a file.
  #
  kinship_table_file = string(keyword["kinship_output_file"])
  CSV.write(kinship_table_file, kinship_frame;
    header = true, delim = keyword["output_field_separator"],
    missingstring = keyword["output_missing_value"])

  return kinship_frame
end # function compare_kinships

"""
    theoretical_coefficients(pedigree::Pedigree, person::Person, keyword::Dict{AbstractString, Any}

This function orchestrates the computation via gene dropping
of the Kinship and Delta7 coefficients deterministically
and Jacquard's 9 identity coefficients stochastically.
The results are placed in a dataframe.
"""
function theoretical_coefficients(pedigree::Pedigree, person::Person,
  keyword::Dict{AbstractString, Any})

  pedigrees = pedigree.pedigrees
  repetitions = keyword["repetitions"]
  xlinked = keyword["xlinked_analysis"]
  #
  # Define dataframes for theoretical Kinship, Delta7, and Jacquard
  # coefficients.
  #
  kinship_dataframe = DataFrame(ped1 = Int[], Pedigree = String[],
    Person1 = String[], Person2 = String[], Kinship = Float64[])
  if !xlinked
    delta_dataframe = DataFrame(ped2 = Int[], Pedigree = String[],
      Person1 = String[], Person2 = String[], Delta7 = Float64[])
  end
  coefficient_dataframe = DataFrame(ped3 = Int[], Pedigree = String[],
    Person1 = String[], Person2 = String[],
    delta1 = Float64[], delta2 = Float64[], delta3 = Float64[],
    delta4 = Float64[], delta5 = Float64[], delta6 = Float64[],
    delta7 = Float64[], delta8 = Float64[], delta9 = Float64[])
  #
  # Loop over all pedigrees.
  #
  for ped = 1:pedigrees
    #
    # Retrieve the various coefficients.
    #
    kinship = kinship_matrix(pedigree, person, ped, xlinked)
    if !xlinked
      delta7 = delta7_matrix(pedigree, person, kinship, ped)
    end
    coefficient = jacquard_coefficients(pedigree, person, ped,
      repetitions, xlinked)
    #
    # Insert the computed coefficients into the appropriate dataframes.
    #
    q = pedigree.pedstart[ped] - 1
    for i = 1:pedigree.individuals[ped]
      for j = i:pedigree.individuals[ped]
        push!(kinship_dataframe, [ped, pedigree.name[ped],
          person.name[i + q], person.name[j + q], kinship[i, j]])
        if !xlinked
          push!(delta_dataframe, [ped, pedigree.name[ped],
            person.name[i + q], person.name[j + q], delta7[i, j]])
        end
        push!(coefficient_dataframe, [ped, pedigree.name[ped],
          person.name[i + q], person.name[j + q],
          coefficient[1, i, j], coefficient[2, i, j], coefficient[3, i, j],
          coefficient[4, i, j], coefficient[5, i, j], coefficient[6, i, j],
          coefficient[7, i, j], coefficient[8, i, j], coefficient[9, i, j]])
      end
    end
  end
  #
  # Join the separate dataframes, and sort the combined dataframe so
  # that relative pairs are lumped by pedigree..
  #
  if xlinked
    combined_dataframe = join(kinship_dataframe, coefficient_dataframe,
      on = [:Pedigree, :Person1, :Person2])
    sort!(combined_dataframe, [:ped1, :Person1, :Person2])
##    deletecols!(combined_dataframe, [:ped1, :ped3])
    names_list = names(combined_dataframe)
    deleteat!(names_list, findall((in)([:ped1, :ped3]), names_list))
    select!(combined_dataframe, names_list)
  else
    temp_dataframe = join(kinship_dataframe, delta_dataframe,
      on = [:Pedigree, :Person1, :Person2])
    combined_dataframe = join(temp_dataframe, coefficient_dataframe,
      on = [:Pedigree, :Person1, :Person2])
    sort!(combined_dataframe, [:ped1, :Person1, :Person2])
##    deletecols!(combined_dataframe, [:ped1, :ped2, :ped3])
    names_list = names(combined_dataframe)
    deleteat!(names_list, findall((in)([:ped1, :ped2, :ped3]), names_list))
    select!(combined_dataframe, names_list)
  end
  #
  # Write the combined coefficient frame to a file and return.
  #
  kinship_table_file = string(keyword["kinship_output_file"])
  CSV.write(kinship_table_file, combined_dataframe;
    header = true, delim = keyword["output_field_separator"],
    missingstring = keyword["output_missing_value"])
  return combined_dataframe
end # function theoretical_coefficients

"""
    kinship_matrix(pedigree::Pedigree, person::Person, ped::Int, xlinked::Bool)

This function recursively computes a kinship matrix for each pedigree.
Children must appear after parents. Each pedigree occurs in a block
demarked by its start and finish and is defined implicitly by
the provided fathers and mothers. Person i is a pedigree founder
if and only if mother[i] = 0. X-linked inheritance is allowed.
"""
function kinship_matrix(pedigree::Pedigree, person::Person,
  ped::Int, xlinked::Bool)
  #
  # Allocate a kinship matrix and an offset.
  #
  pedstart = pedigree.pedstart[ped]
  pedfinish = pedigree.pedfinish[ped]
  kinship = zeros(Float64, pedigree.individuals[ped], pedigree.individuals[ped])
  q = pedstart - 1
  for i = pedstart:pedfinish
    #
    # Initialize the kinship coefficients for founders.
    #
    j = person.mother[i]
    if j == 0
      if xlinked && person.male[i]
        kinship[i - q, i - q] = 1.0
      else
        kinship[i - q, i - q] = 0.5
      end
    else
      #
      # Compute the kinship coefficients of descendants by averaging.
      #
      if xlinked && person.male[i]
        kinship[i - q, i - q] = 1.0
        for m = pedstart:i - 1
          kinship[i - q, m - q] = kinship[j - q, m - q]
          kinship[m - q, i - q] = kinship[i - q, m - q]
        end
      else
        k = person.father[i]
        kinship[i - q, i - q] = 0.5 + 0.5 * kinship[j - q, k - q]
        for m = pedstart:i - 1
          kinship[i - q, m - q] = 0.5 * (kinship[j - q, m - q]
             + kinship[k - q, m - q])
          kinship[m - q, i - q] = kinship[i - q, m - q]
        end
      end
    end
  end
  #
  # Extend the kinship matrix to co-twins.
  #
  cotwin_extension!(kinship, pedigree, person, ped)
  return kinship
end # function kinship_matrix

"""
    delta7_matrix(pedigree::Pedigree, person::Person, kinship::Matrix{Float64}, ped::Int)

This function computes the Delta7 matrix for non-inbred pedigrees.
Children must appear after parents. Each pedigree occurs in a block
and is defined implicitly by the provided fathers and mothers.
Person i is a pedigree founder if and only if mother[i] = 0.
"""
function delta7_matrix(pedigree::Pedigree, person::Person,
  kinship::Matrix{Float64}, ped::Int)

  pedstart = pedigree.pedstart[ped]
  pedfinish = pedigree.pedfinish[ped]
  delta7 = zeros(Float64, pedigree.individuals[ped], pedigree.individuals[ped])
  q = pedstart - 1
  for i = pedstart:pedfinish
    for j = pedstart:i
      if j == i
        delta7[i - q, i - q] = 1.0
      elseif person.mother[i] > 0 && person.mother[j] > 0
        k = person.father[i]
        l = person.mother[i]
        m = person.father[j]
        n = person.mother[j]
        delta7[i - q, j - q] = kinship[k - q, m - q] * kinship[l - q, n - q] +
          kinship[k - q, n - q] * kinship[l - q, m - q]
        delta7[j - q, i - q] = delta7[i - q, j - q]
      end
    end
  end
  #
  # Extend the delta7 matrix to co-twins.
  #
  cotwin_extension!(delta7, pedigree, person, ped)
  return delta7
end # function delta7_matrix

"""
    jacquard_coefficients(pedigree, person, ped::Int, replicates::Int, xlinked::Bool)

This function computes Jacquard's nine condensed identity
coefficients for each pedigree by simulation. Children
must appear after parents. A pedigree is defined implicitly
by giving fathers and mothers. Person i is a pedigree founder
if and only if mother[i] = 0. The more replicates requested,
the more precision attained. X-linked inheritance is possible.
"""
function jacquard_coefficients(pedigree::Pedigree, person::Person,
  ped::Int, replicates::Int, xlinked::Bool)

  pedstart = pedigree.pedstart[ped]
  pedfinish = pedigree.pedfinish[ped]
  ped_size = pedigree.individuals[ped]
  coefficient = zeros(Float64, 9, ped_size, ped_size)
  x = zeros(Float64, ped_size, ped_size)
  source = zeros(Int, 2, ped_size)
  #
  # Loop over all replicates.
  #
  q = pedstart - 1
  for rep = 1:replicates
    f = 0
    for j = pedstart:pedfinish
      if person.mother[j] == 0
        f = f + 1
        source[1, j - q] = f
        if xlinked && person.male[j]
        else
          f = f + 1
        end
        source[2, j - q] = f
      else
        if xlinked && person.male[j]
          u = rand(1:2, 1)
          source[1, j - q] = source[u[1], person.mother[j] - q]
          source[2, j - q] = source[1, j - q]
        else
          u = rand(1:2, 2)
          source[1, j - q] = source[u[1], person.father[j] - q]
          source[2, j - q] = source[u[2], person.mother[j] - q]
        end
      end
    end
    #
    # Loop over all pairs.
    #
    for j = pedstart:pedfinish
      for k = pedstart:pedfinish
        l = identity_state(source[:, j - q], source[:, k - q])
        coefficient[l, j - q, k - q] = coefficient[l, j - q, k - q] + 1.0
      end
    end
  end
  #
  # Extend the coefficients to co-twins.
  #
  for i = 1:9
    x = reshape(coefficient[i, :, :], ped_size, ped_size)
    cotwin_extension!(x, pedigree, person, ped)
    coefficient[i, :, :] = reshape(x, 1, ped_size, ped_size)
  end
  return coefficient / replicates
end # function jacquard_coefficients

"""
    cotwin_extension!(x::Matrix{Float64}, pedigree::Pedigree, person::Person, ped::Int)

This function extends identity coefficients to co-twins.
"""
function cotwin_extension!(x::Matrix{Float64}, pedigree::Pedigree,
  person::Person, ped::Int)

  q = pedigree.pedstart[ped] - 1
  for i = pedigree.pedfinish[ped] + 1:pedigree.twin_finish[ped]
    m = person.primary_twin[i]
    for j = pedigree.pedstart[ped]:pedigree.pedfinish[ped]
      x[i - q, j - q] = x[m - q, j - q]
      x[j - q, i - q] = x[i - q, j - q]
    end
    for j = pedigree.pedfinish[ped] + 1:pedigree.twin_finish[ped]
      n = person.primary_twin[j]
      x[i - q, j - q] = x[m - q, n - q]
    end
  end
end # function cotwin_extension!

"""
    identity_state(source1::Vector{Int}, source2::Vector{Int})

This function assigns one of the 9 condensed identity states
to two sourced genotypes. The integer n is the number of
unique genes among the 4 sampled genes.
"""
function identity_state(source1::Vector{Int}, source2::Vector{Int})

  n = length(union(BitSet(source1), BitSet(source2)))
  if n == 4
    return 9
  elseif n == 3
    if source1[1] == source1[2]
      return 4
    elseif source2[1] == source2[2]
      return 6
    else
      return 8
    end
  elseif n == 2
    if source1[1] == source1[2]
      if source2[1] == source2[2]
        return 2
      else
        return 3
      end
    else
      if source2[1] == source2[2]
        return 5
      else
        return 7
      end
    end
  else
    return 1
  end
end # function identity_state

"""
This function assembles pairs of names and values for printing.
"""
function make_plot_pairs(x::DataFrame)

  plot_pair = Array{String}(undef, size(x, 1))
  for i in 1:length(plot_pair)
    plot_pair[i] = "Person1=" * x[i, 3] * ", " * "Person2=" * x[i, 4]
  end
  return plot_pair
end # function make_plot_pairs

"""
    make_compare_plot(x::DataFrame, plot_pair::Vector{String})

This function composes the comparison plot using PlotlyJS.
"""
function make_compare_plot(x::DataFrame, plot_pair::Vector{String})

  trace1 = scatter(;x=x[!, :theoretical_kinship],
    y=x[!, :empiric_kinship], mode="markers",
    name="empiric kinship", text=plot_pair)
  trace2 = scatter(;x=[1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 0.0],
    y=[1/2, 1/4, 1/8, 1/16, 1/32, 1/64, 1/128, 0.0],
    mode="markers", name="marker for midpoint")
  layout = Layout(;title="Compare Empiric and Theoretical Kinships",
    hovermode="closest",
    xaxis=attr(title="Theoretical Kinship",
      showgrid=false, zeroline=false),
    yaxis=attr(title="Empiric Kinship", zeroline=false))
  data = [trace1, trace2]
  plot(data, layout)
end # function make_compare_plot

"""
    compute_fishers_z(x::DataFrame)

This function computes Fisher's atanh transform and standardizes.
"""
function compute_fishers_z(x::DataFrame)

  theoretical_transformed = map(atanh, x[!, :theoretical_kinship])
  empiric_transformed = map(atanh, x[!, :empiric_kinship])
  difference = empiric_transformed - theoretical_transformed
  new_zscore = zscore(difference) #in StatsBase package
  return [x DataFrame(fishers_zscore = new_zscore)]
end # function compute_fishers_z

"""
    plot_fisher_z(x::DataFrame, plot_pair::Vector{String})

This function plots the transformed data using PlotlyJS.
"""
function plot_fisher_z(x::DataFrame, plot_pair::Vector{String})

    trace1 = histogram(x=x[!, :fishers_zscore], text=plot_pair)
    data = [trace1]
    layout = Layout(barmode="overlay",
        title="Z-score plot for Fisher's statistic",
        xaxis=attr(title="Standard Deviations"),
        yaxis=attr(title="Count"))
    plot(data, layout)
end # function plot_fisher_z

"""

"""
function correspond(x::AbstractVector{AbstractString},
  y::AbstractVector{AbstractString})

  (m, n) = (length(x), length(y))
  xperm = sortperm(x)
  yperm = sortperm(y)
  x_to_y = zeros(Int, m)
  y_to_x = zeros(Int, n)
  (i, j) = (1, 1)
  done = false
  while !done
    ii = xperm[i]
    jj = yperm[j]
    if x[ii] == y[jj]
      x_to_y[ii] = jj
      y_to_x[jj] = ii
      (i, j) = (i + 1, j + 1)
    elseif x[ii] < y[jj]
      i = i +1
    elseif y[jj] < x[ii]
      j = j + 1
    end
    done = i > m || j > n
  end
  return (x_to_y, y_to_x)
end # function correspond

# x = ["the", "and", "a", "be", "an"]
# y = ["but", "be", "a", "an", "so"]
# (x_to_y, y_to_x) = correspond(x, y)

"""

"""
function compute_full_person_and_pedigree(keyword :: Dict{AbstractString, Any},
  subset_person :: Person)
#
# Use MendelBase's function to construct the "full" pedigree and person data.
# First, store original key values.
#
  original_pedigree_file = keyword["pedigree_file"]
  snpdefinition_file = keyword["snpdefinition_file"]
  snpdata_file = keyword["snpdata_file"]

  keyword["pedigree_file"] = keyword["full_pedigree_file"]
  keyword["snpdefinition_file"] = "" # Avoid reading in snpdata again.
  keyword["snpdata_file"] = ""
  full_pedigree, full_person, = read_external_data_files(keyword)
  keyword["pedigree_file"] = original_pedigree_file
  keyword["snpdefinition_file"] = snpdefinition_file
  keyword["snpdata_file"] = snpdata_file
#
# Now compute indicator of person i having SNP information.
#
  total_people = full_person.people
  contains_snp = falses(total_people)
  for i in 1:total_people
    person_i = full_person.name[i]
    (person_i in subset_person.name) && (contains_snp[i] = true)
  end

  return full_pedigree, full_person, contains_snp
end # function compute_full_person_and_pedigree

#
# Method to obtain path to this package's data files
# so they can be used in the documentation and testing routines.
# For example, datadir("Control file.txt") will return
# "/path/to/package/data/Control file.txt"
#
datadir(parts...) = joinpath(@__DIR__, "..", "data", parts...)

end # module MendelKinship


text:
#!/usr/bin/julia
#
# File containing inner product operations on the grid space.
#
export ip, norm, norm2, normalise!, orthog!, orthog;

"""
    do_ip(v::AbstractVector, w::AbstractVector, Δ::Number)

Implementation of inner product calculation.
"""
function do_ip(v::AbstractVector, w::AbstractVector, Δ::Number)
  # Approximate integral over the discrete space.
  return sum(v' * w) * Δ;
end

"""
    check_ip_compatable(ψ::WaveObject{T,N,S,M,G},
                        ϕ::WaveObject{T,K,S,M,G}) where {T,N,S,M,K,G<:UniformGrid}

Check the lengths and uniform grid spacings of the `WaveObjects` are compatible
to allow the calculation of an inner product.
"""
function check_ip_compatable(ψ::WaveObject{T,N,S,M,G},
                             ϕ::WaveObject{P,K,S,M,G}) where {T,P,N,S,M,K,G<:UniformGrid}

  if (N === K) check_same_size(ψ, ϕ);
  else check_ip_compatable_dim(ψ,ϕ); end

  if (ψ.Δ != ϕ.Δ) error("grid spacings of WaveObjects are different"); end
end

"""
    check_same_size(ψ::WaveObject{T,N}, ϕ::WaveObject{P,N}) where {T,P,N}

Check `WaveObjects` of the same dimension have the same size.
"""
function check_same_size(ψ::WaveObject{T,N}, ϕ::WaveObject{P,N}) where {T,P,N}
  if (size(ψ) != size(ϕ))
    throw((DimensionMismatch, "wave functions must be of same size."));
  end
end
function check_same_size(ψ::WaveObject{T,1}, ϕ::WaveObject{P,1}) where {T,P}
  if (length(ψ) != length(ϕ))
    throw(DimensionMismatch("wave functions must be of same size."));
  end
end

"""
    check_ip_compatable_dim(ψ::WaveObject{T,1}, ϕ::WaveObject{P,2}) where {T,P}

Check that the length of the one-dimensional `WaveObject` is the same as the
length of a column of a two-dimensional `WaveObject`.
"""
function check_ip_compatable_dim(ψ::WaveObject{T,1}, ϕ::WaveObject{P,2}) where {T,P}
  if (length(ψ) != length(ϕ[:,1]))
    throw(DimensionMismatch("wave function must be same length as matrix column."));
  end
end
function check_ip_compatable_dim(ϕ::WaveObject{P,2}, ψ::WaveObject{T,1}) where {T,P}
  check_ip_compatable_dim(ψ, ϕ);
end

"""
    ip(ψ::O, ϕ::O) where {T,S,G<:UniformGrid,O<:WaveFunction{T,1,S,G}}

Inner product between two wave functions on the same one-dimensional
`UniformGrid` system.
"""
function ip(ψ::O, ϕ::O) where {T,S,G<:UniformGrid,O<:WaveFunction{T,1,S,G}}
  check_ip_compatable(ψ, ϕ);
  return do_ip(ψ, ϕ, ψ.Δ);
end

"""
    ip(ψ::WaveFunction{T,1,S,G}, ϕ::WaveMatrix{T,S,G}, pos::Int) where
            {T,S,G<:UniformGrid}

Inner product between `WaveFunction` and wave function in column `pos` of the
matrix-like object `ϕ`.

Objects are defined on a one-dimensional `UniformGrid`.
"""
function ip(ψ::WaveFunction{T,1,S,G}, ϕ::WaveMatrix{T,S,G}, pos::Int) where
            {T,S,G<:UniformGrid}
  check_ip_compatable(ψ, ϕ);
  return do_ip(ψ, ϕ[:,pos], ψ.Δ);
end
function ip(ϕ::WaveMatrix{T,S,G}, pos::Int, ψ::WaveFunction{T,1,S,G}) where
            {T,S,G<:UniformGrid}
  ip(ψ, ϕ, pos);
end

"""
    ip(ψ::O, pos::Int, ϕ::O, pos2::Int) where
          {T,S,G<:UniformGrid,O<:WaveMatrix{T,2,S,1,G}}

Inner product between columns of two wave matrices.
"""
function ip(ψ::O, pos::Int, ϕ::O, pos2::Int) where
          {T,S,G<:UniformGrid,O<:WaveMatrix{T,S,G}}
  check_ip_compatable(ψ, ϕ);
  return do_ip(ψ[:,pos], ϕ[:,pos2], ψ.Δ);
end

"""
    norm2(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid}

Calculate norm² of a one-dimensional wave function on a `UniformGrid`.
"""
function norm2(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid}
  return sum(abs2.(ψ) * ψ.Δ);
end

"""
    norm(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid}

Calculate norm of a one-dimensional wave function on a `UniformGrid`.
"""
norm(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid} = √norm2(ψ);

LinearAlgebra.norm(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid} =
  QuantumSimulation.norm(ψ);

"""
    norm2(ψ::WaveMatrix{T,S,G}, pos::Int) where {T,S,G<:UniformGrid}

Calculate norm² of a matrix `WaveMatrix` element in column `pos` defined on a
one-dimensional `UniformGrid`.
"""
function norm2(ϕ::WaveMatrix{T,S,G}, pos::Int) where {T,S,G<:UniformGrid}
  return sum(abs2.(ϕ[:,pos]) * ϕ.Δ);
end

"""
    norm(ψ::WaveMatrix{T,S,G}, pos::Int) where {T,S,G<:UniformGrid}

Calculate norm of a matrix `WaveMatrix` element in column `pos` defined on a
one-dimensional `UniformGrid`.
"""
function norm(ϕ::WaveMatrix{T,S,G}, pos::Int) where {T,S,G<:UniformGrid}
  return √norm2(ϕ, pos);
end

"""
    normalise!(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid}

Normalise the passed wave function defined on a one-dimensional `UniformGrid`.
"""
function normalise!(ψ::WaveFunction{T,1,S,G}) where {T,S,G<:UniformGrid}
  ψ ./= norm(ψ);
end

"""
    orthog!(ψ::WaveFunction{T,1,S,G}, ϕ::WaveMatrix{T,S,G},
                 max::Int=length(ϕ[1,:])) where {T,S,G<:UniformGrid}

Orthogonalise `ψ` with respect to the first `max` unit normalised wave functions
in ϕ.

Both `ψ` and `ϕ` are defined on a one-dimensional `UniformGrid`.

See also: [`orthog`](@ref)
"""
function orthog!(ψ::WaveFunction{T,1,S,G}, ϕ::WaveMatrix{T,S,G},
                 max::Int=length(ϕ[1,:])) where {T,S,G<:UniformGrid}
  # Implement modified Gram Schmidt to prevent accumulation of round-off errors.
  for i ∈ 1:max
    ψ[:] -= ip(ψ, ϕ, i) * ϕ[:,i];
  end
end
"""
    orthog(ψ::WaveFunction{T,1,S,G}, ϕ::WaveMatrix{T,S,G},
                 max::Int=length(ϕ[1,:])) where {T,S,G<:UniformGrid}

Return a copy of `ψ` orthogonalised with respect to the first `max` unit
normalised wave functions in ϕ.

Both `ψ` and `ϕ` are defined on a one-dimensional `UniformGrid`.

See also: [`orthog!`](@ref)
"""
function orthog(ψ::WaveFunction{T,1,S,G}, ϕ::WaveMatrix{T,S,G},
                 max::Int=length(ϕ[1,:])) where {T,S,G<:UniformGrid}
  ψcopy = deepcopy(ψ);
  orthog!(ψcopy, ϕ, max);
  return ψcopy;
end


text:
#!/usr/bin/env julia

using Pkg

# add the packages to the current environment, defined by JULIA_PROJECT
Pkg.add(PackageSpec(url = "https://github.com/iancze/DiskJockey.git"))
Pkg.add("ArgParse")
Pkg.add("HDF5")
Pkg.add("PyPlot")
Pkg.add("LaTeXStrings")
Pkg.add("YAML")
Pkg.add("Printf")
Pkg.add("PyCall")
Pkg.add("Images")
Pkg.add("NaNMath")
Pkg.add("NPZ")
Pkg.add("Statistics")

Pkg.instantiate()