text:
function __init__()
    register(DataDep("WordNet 3.0",
        """
        Dataset: WordNet 3.0
        Website: https://wordnet.princeton.edu/wordnet

        George A. Miller (1995). WordNet: A Lexical Database for English.
        Communications of the ACM Vol. 38, No. 11: 39-41.

        Christiane Fellbaum (1998, ed.) WordNet: An Electronic Lexical Database.
        Cambridge, MA: MIT Press.

        License:
        WordNet Release 3.0 This software and database is being provided to you, the LICENSEE, by Princeton University under the following license. By obtaining, using and/or copying this software and database, you agree that you have read, understood, and will comply with these terms and conditions.: Permission to use, copy, modify and distribute this software and database and its documentation for any purpose and without fee or royalty is hereby granted, provided that you agree to comply with the following copyright notice and statements, including the disclaimer, and that the same appear on ALL copies of the software, database and documentation, including modifications that you make for internal use or for distribution. WordNet 3.0 Copyright 2006 by Princeton University. All rights reserved. THIS SOFTWARE AND DATABASE IS PROVIDED "AS IS" AND PRINCETON UNIVERSITY MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. BY WAY OF EXAMPLE, BUT NOT LIMITATION, PRINCETON UNIVERSITY MAKES NO REPRESENTATIONS OR WARRANTIES OF MERCHANT- ABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE, DATABASE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS. The name of Princeton University or Princeton may not be used in advertising or publicity pertaining to distribution of the software and/or database. Title to copyright in this software, database and any associated documentation shall at all times remain with Princeton University and LICENSEE agrees to preserve same.
        """,
        "http://wordnetcode.princeton.edu/3.0/WNdb-3.0.tar.gz",
        "658b1ba191f5f98c2e9bae3e25c186013158f30ef779f191d2a44e5d25046dc8";
        post_fetch_method = unpack
    ))
end


text:
include("FDI/FDI.jl")
include("multicopters/multicopters.jl")
include("allocators/allocators.jl")
include("integrated_environments/integrated_environments.jl")

ext: jl


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function spatial_varying_tangent_elastic(mu,m,n,h,type=1)
    spatial_varying_tangent_elastic_ = load_op_and_grad("./build/libSpatialVaryingTangentElastic","spatial_varying_tangent_elastic")
    mu,m_,n_,h,type = convert_to_tensor([mu,m,n,h,type], [Float64,Int64,Int64,Float64,Int64])
    H = spatial_varying_tangent_elastic_(mu,m_,n_,h,type)
    set_shape(H, (4*m*n, 2, 2))
end

# TODO: specify your input parameters
m = 10
n = 5
h = 0.1
mu = rand(4*m*n)
H = zeros(4*m*n, 2, 2)

for i = 1:4*m*n
    H[i,:,:] = mu[i] * diagm(0=>ones(2))
end
type = 1
u = spatial_varying_tangent_elastic(mu,m,n,h,type)
sess = Session(); init(sess)
# @show run(sess, u)-H

m = 10
n = 5
h = 0.1
mu = [rand(4*m*n); rand(4*m*n)]
H = zeros(4*m*n, 2, 2)

for i = 1:4*m*n
    H[i,:,:] = [mu[i] 0.0; 0.0 mu[i+4*m*n]]
end
type = 2
u = spatial_varying_tangent_elastic(mu,m,n,h,type)
sess = Session(); init(sess)
@show run(sess, u)-H

m = 10
n = 5
h = 0.1
mu = rand(4*m*n*3)
H = zeros(4*m*n, 2, 2)

for i = 1:4*m*n
    H[i,:,:] = [mu[i] mu[i+8*m*n]; mu[i+8*m*n] mu[i+4*m*n]]
end
type = 3
u = spatial_varying_tangent_elastic(mu,m,n,h,type)
sess = Session(); init(sess)
@show run(sess, u)-H

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m_)
    return sum(spatial_varying_tangent_elastic(m_,m,n,h,type)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(4*m*n*3))
v_ = rand(4*m*n*3)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using Phylo.API
using Compat: mapreduce
import LightGraphs: src, dst, indegree, outdegree, degree

# AbstractTree/Node/Branch type methods
"""
    treetype(::Type{AbstractTree})

Returns tree number (OneTree, ManyTrees) from a tree type.
"""
treetype(::Type{T}) where {TT <: TreeType, RT, NL, N, B,
    T <: AbstractTree{TT, RT, NL, N, B}} = TT

"""
    roottype(::Type{AbstractTree})
    roottype(::Type{AbstractNode})
    roottype(::Type{AbstractBranch})

Returns root type from a tree type.
"""
roottype(::Type{T}) where {TT, RT <: Rootedness, NL, N, B,
                           T <: AbstractTree{TT, RT, NL, N, B}} = RT
roottype(::Type{N}) where {RT <: Rootedness, NL,
                           N <: AbstractNode{RT, NL}} = RT
roottype(::Type{B}) where {RT <: Rootedness, NL,
                           B <: AbstractBranch{RT, NL}} = RT

"""
    nodenametype(::Type{AbstractTree})
    nodenametype(::Type{AbstractNode})
    nodenametype(::Type{AbstractBranch})

Returns type of node names from a tree type.
"""
nodenametype(::Type{T}) where {TT, RT, NL, N, B,
    T <: AbstractTree{TT, RT, NL, N, B}} = NL
nodenametype(::Type{N}) where {RT <: Rootedness, NL,
                               N <: AbstractNode{RT, NL}} = NL
nodenametype(::Type{B}) where {RT <: Rootedness, NL,
                               B <: AbstractBranch{RT, NL}} = NL

"""
    nodetype(::Type{AbstractTree})

Returns type of nodes from a tree type.
"""
nodetype(::Type{T}) where {TT, RT, NL, N <: AbstractNode, B,
    T <: AbstractTree{TT, RT, NL, N, B}} = N

"""
    branchnametype(::AbstractTree)

Returns type of branch names from a branch type.
"""
branchnametype(::Type{<: AbstractTree}) = Int
branchnametype(::Type{<: AbstractNode}) = Int
branchnametype(::Type{<: AbstractBranch}) = Int

"""
    branchtype(::Type{AbstractTree})

Returns type of branches from a tree type.
"""
branchtype(::Type{T}) where {TT, RT, NL, N, B <: AbstractBranch,
    T <: AbstractTree{TT, RT, NL, N, B}} = B

"""
    treenametype(::Type{AbstractTree})

Returns the name type for a tree type.
"""
treenametype(treetype::Type{<: AbstractTree}) = _treenametype(treetype)

"""
    leafinfotype(::Type{<: AbstractTree})

retrieve the leaf info type of a tree.
"""
leafinfotype(::Type{<: AbstractTree}) = _leafinfotype(T)

"""
    noderecordtype(::Type{<: AbstractTree})

retrieve the node info type of a tree.
"""
noderecordtype(::Type{<: AbstractTree}) = _noderecordtype(T)

"""
    branchrecordtype(::Type{<: AbstractTree})

retrieve the branch info type of a tree.
"""
branchrecordtype(::Type{<: AbstractTree}) = _branchrecordtype(T)

# AbstractTree methods
"""
    ntrees(tree::AbstractTree)

Returns the number of trees in a tree object, 1 for a OneTree tree type, and
the count of trees for a ManyTrees type.
"""
ntrees(tree::AbstractTree) = _ntrees(tree)

"""
    gettrees(tree::AbstractTree)

Returns an iterable of trees.
"""
gettrees(tree::AbstractTree) = _gettrees(tree)

"""
    gettree(tree::AbstractTree, label)

Returns a single OneTree object tree corresponding to the label given.
"""
function gettree end
gettree(tree::AbstractTree{OneTree}) = tree
function gettree(tree::AbstractTree{ManyTrees})
    @assert _ntrees(tree) == 1 "Must be only one tree (found $(_ntrees(tree)))"
    return first(_gettrees(tree))
end
gettree(tree::AbstractTree, label) = _gettree(tree, label)

"""
    gettreenames(tree::AbstractTree)

Returns the names of the trees.
"""
function gettreenames end
gettreenames(tree::AbstractTree{OneTree}) = [_gettreename(tree)]
gettreenames(tree::AbstractTree{ManyTrees}) = _gettreenames(tree)

"""
    gettreename(tree::AbstractTree)

Returns the name of the single tree.
"""
gettreename(tree::AbstractTree{OneTree}) = _gettreename(tree)
function gettreename(tree::AbstractTree{ManyTrees})
    @assert _ntrees(tree) == 1 "Must be only one tree (found $(_ntrees(tree)))"
    return first(_gettreenames(tree))
end

"""
    nleaves(::AbstractTree)

Returns the number of leaves (tips) in a tree.
"""
nleaves(tree::AbstractTree) = _nleaves(tree)

"""
    nroots(::AbstractTree)

Returns the number of roots in a tree. For OneTree types, Unrooted trees will
return 0, OneRoot trees should return 1, and manyroots tree (ones with multiple
subtrees) will return the number of subtrees. Manytrees types will return a
Dict of counts of the number of roots for each tree in the set.
"""
function nroots end
nroots(::AbstractTree{OneTree, Unrooted}) = 0
nroots(tree::AbstractTree{OneTree, <: Rooted}) = _nroots(tree)
nroots(tree::AbstractTree{ManyTrees}) =
    Dict(name => nroots(gettree(tree, name)) for name in _gettreenames(tree))

"""
    getroots(::AbstractTree)
    getroots(::AbstractTree, id)

Returns a vector containing the root(s) of a single (OneTree) tree
or a set of (ManyTrees) trees.
"""
function getroots end
getroots(tree::AbstractTree{OneTree}) = _getroots(tree)
getroots(tree::AbstractTree{OneTree, Unrooted, NL, N, B}) where {NL, N, B} = N[]
getroots(tree::AbstractTree{ManyTrees}) =
    Dict(name => getroots(gettree(tree, name)) for name in _gettreenames(tree))

"""
    getroot(::AbstractTree)

Returns the root of a single tree (must be only one tree for a ManyTrees tree).
"""
function getroot end
getroot(tree::AbstractTree{OneTree, <: Rooted}) = _getroot(tree)
getroot(tree::AbstractTree{ManyTrees, <: Rooted}) = getroot(gettree(tree))

"""
    getnodes(::AbstractTree)

Returns the vector of nodes of a single tree, or a Dict of vectors of nodes
for multiple trees.
"""
function getnodes end
getnodes(tree::AbstractTree{OneTree}) = _getnodes(tree)
getnodes(tree::AbstractTree{ManyTrees}) =
    Dict(name => getnodes(gettree(tree, name)) for name in _gettreenames(tree))

"""
    nnodes(::AbstractTree)

Returns the number of nodes of a single tree, or a Dict of numbers of nodes
for multiple trees.
"""
function nnodes end
nnodes(tree::AbstractTree{OneTree}) = _nnodes(tree)
nnodes(tree::AbstractTree{ManyTrees}) =
    Dict(name => nnodes(gettree(tree, name)) for name in _gettreenames(tree))

"""
    nbranches(::AbstractTree)

Returns the number of branches of a single tree, or a Dict of numbers of
branches for multiple trees.
"""
function nbranches end
nbranches(tree::AbstractTree{OneTree}) = _nbranches(tree)
nbranches(tree::AbstractTree{ManyTrees}) =
    Dict(name => nbranches(gettree(tree, name)) for name in _gettreenames(tree))

"""
    getnodenames(tree::AbstractTree)

Return a vector of node names of a single tree (identified by id for a
ManyTrees tree), or a Dict of vectors of node names for multiple trees.
"""
function getnodenames end
getnodenames(tree::AbstractTree{OneTree}) = _getnodenames(tree)
getnodenames(tree::AbstractTree{ManyTrees}) =
    Dict(name => getnodenames(gettree(tree, name))
         for name in _gettreenames(tree))

"""
    getbranches(::AbstractTree)

Returns the vector of branches of a single tree, or a Dict of vectors of
branches for multiple trees.
"""
function getbranches end
getbranches(tree::AbstractTree{OneTree}) = _getbranches(tree)
getbranches(tree::AbstractTree{ManyTrees}) =
    Dict(name => getbranches(gettree(tree, name))
         for name in _gettreenames(tree))

"""
    getbranchnames(tree::AbstractTree)

Return a vector of branch names of a single tree, or a Dict of vectors of
branch names for multiple trees.
"""
function getbranchnames end
getbranchnames(tree::AbstractTree{OneTree}) = _getbranchnames(tree)
getbranchnames(tree::AbstractTree{ManyTrees}) =
    Dict(name => getbranchnames(gettree(tree, name))
         for name in _gettreenames(tree))

"""
    createbranch!(tree::AbstractTree, source, destination[, length::Float64];
                  data)

Add a branch from `source` to `destination` on `tree` with optional length and
data. `source` and `destination` can be either nodes or nodenames.
"""
function createbranch!(tree::AbstractTree{OneTree, <: Rooted},
                       source, destination, length::Float64 = NaN;
                       name::Int = _newbranchlabel(tree), data = missing)
    sn = _getnode(tree, source)
    dn = _getnode(tree, destination)
    _hasnode(tree, sn) ||
        error("Tree does not have an available source node called " *
              _getnodename(tree, source))
    _hasoutboundspace(tree, sn) ||
        error(_getnodename(tree, source) *
              " already has maximum number of outbound connections ($(_outdegree(tree, sn)))")
    _hasnode(tree, dn) ||
        error("Tree does not have a destination node called " *
              _getnodename(tree, destination))
    _hasinboundspace(tree, dn) ||
        error("Tree does not have an available destination node called " *
              _getnodename(tree, destination))
    sn !== dn || error("Branch must connect different nodes")

    return ismissing(data) ? _createbranch!(tree, sn, dn, length, name) :
        _createbranch!(tree, sn, dn, length, name, data)
end
function createbranch!(tree::AbstractTree{OneTree, Unrooted},
                       source, destination, length::Float64 = NaN;
                       name::Int = _newbranchlabel(tree), data = missing)
    sn = _getnode(tree, source)
    dn = _getnode(tree, destination)
    _hasnode(tree, sn) ||
        error("Tree does not have an available source node called " *
              _getnodename(tree, source))
    _hasspace(tree, sn) ||
        error(_getnodename(tree, source) *
              " already has maximum number of connections " *
              "$(_degree(tree, sn)))")
    _hasnode(tree, dn) ||
        error("Tree does not have a destination node called " *
              _getnodename(tree, destination))
    _hasspace(tree, dn) ||
        error(_getnodename(tree, destination) *
              " already has maximum number of connections " *
              "$(_degree(tree, dn)))")
    sn !== dn || error("Branch must connect different nodes")

    return ismissing(data) ? _createbranch!(tree, sn, dn, length, name) :
        _createbranch!(tree, sn, dn, length, name, data)
end

"""
    deletebranch!(tree::AbstractTree, branch)
    deletebranch!(tree::AbstractTree, src, dst)

Delete the branch `branch` from `tree`, or branch connecting `src` node to
`dst` node.
"""
function deletebranch! end
function deletebranch!(tree::AbstractTree{OneTree}, branch)
    _hasbranch(tree, branch) ||
        error("Tree does not have a branch called " *
              "$(_getbranchname(tree, branch))")
    return _deletebranch!(tree, branch)
end
function deletebranch!(tree::AbstractTree{OneTree, RT, NL, N, B},
                       source::Union{N, NL}, dest::Union{N, NL}) where
    {RT, NL, N, B}
    _hasbranch(tree, source, dest) ||
        error("Tree does not have a branch between " *
              "$(_getnodename(tree, source)) and $(_getnodename(tree, dest))")
    return _deletebranch!(tree, _getbranch(tree, source, dest))
end

"""
    createnode!(tree::AbstractTree[, nodename]; data)

Create a node on a tree with optional node info.
"""
function createnode!(tree::AbstractTree{OneTree, RT, NL, N, B},
                     nodename::NL = _newnodelabel(tree);
                     data = missing) where {RT, NL, N, B}
    !_hasnode(tree, nodename) ||
        error("Node $nodename already present in tree")
    return ismissing(data) ? _createnode!(tree, nodename) :
        _createnode!(tree, nodename, data)
end

"""
    createnodes!(tree::AbstractTree, count::Integer)
    createnodes!(tree::AbstractTree, nodenames)
    createnodes!(tree::AbstractTree, nodedict)

Add a number of nodes, a vector with given names, or a Dict with node names
and associated node info to a tree.
"""
function createnodes! end
createnodes!(tree::AbstractTree{OneTree}, count::Int) =
    [createnode!(tree) for count in Base.OneTo(count)]
function createnodes!(tree::T, nodenames::V) where
    {RT, NL, N, B,
     T <: AbstractTree{OneTree, RT, NL, N, B},
     V <: AbstractVector{NL}}
    here = [name for name in nodenames if _hasnode(tree, name)]
    isempty(here) ||
        error("Nodes $here already present in tree")
    return [_createnode!(tree, name) for name in nodenames]
end
function createnodes!(tree::T, nodedata::NDS) where
    {RT, NL, N, B, T <: AbstractTree{OneTree, RT, NL, N, B},
     ND, NDS <: Dict{NL, ND}}
    here = [name for name in keys(nodedata) if _hasnode(tree, name)]
    isempty(here) || error("Nodes $here already present in tree")
    noderecordtype(T) == ND || error("Node info types does not match")
    return [_createnode!(tree, info.first, info.second) for info in nodedata]
end

"""
    deletenode!(tree::AbstractTree, node)

Delete a node (or a name) from a tree
"""
function deletenode!(tree::AbstractTree{OneTree}, node)
    _hasnode(tree, node) || error("Tree does not contain node $node")
    return _deletenode!(tree, _getnode(tree, node))
end

"""
    hasnode(tree::AbstractTree, node)

Returns whether a tree has a given node (or node name) or not.
"""
hasnode(tree::AbstractTree{OneTree}, node) = _hasnode(tree, node)

"""
    getnode(tree::AbstractTree, nodename)

Returns a node from a tree.
"""
function getnode(tree::AbstractTree{OneTree}, node)
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getnode(tree, node)
end

"""
    getnodename(::AbstractTree, node)

Returns the node name associated with a node from a tree. For some
node types, it will be able to extract the node name without reference to
the tree.
"""
function getnodename end
function getnodename(tree::AbstractTree{OneTree}, node)
    _hasnode(tree, node) || error("Node $node does not exist")
    return _getnodename(tree, node)
end

"""
    hasbranch(tree::AbstractTree, branch)
    hasbranch(tree::AbstractTree, source, dest)

Does `tree` have a branch `branch` or a branch from `source` to `dest`?
"""
function hasbranch end
hasbranch(tree::AbstractTree{OneTree}, branch) = _hasbranch(tree, branch)
function hasbranch(tree::AbstractTree{OneTree, RT, NL, N, B},
          source::Union{N, NL}, dest::Union{N, NL}) where {RT, NL, N, B}
    _hasnode(tree, source) || error("Node $source does not exist")
    _hasnode(tree, dest) || error("Node $dest does not exist")
    return _hasbranch(tree, source, dest)
end

"""
    getbranch(tree::AbstractTree, branch)
    getbranch(tree::AbstractTree, source, dest)

Returns a branch from a tree by name or by source and destination node.
"""
function getbranch(tree::AbstractTree{OneTree, RT, NL, N, B},
                   branch::Union{Int, B}) where {RT, NL, N, B}
    _hasbranch(tree, branch) || error("Branch $branch does not exist")
    return _getbranch(tree, branch)
end
function getbranch(tree::AbstractTree{OneTree, RT, NL, N, B},
                   source::Union{N, NL}, dest::Union{N, NL}) where {RT, NL, N, B}
    _hasbranch(tree, source, dest) || error("Branch $branch does not exist")
    return _getbranch(tree, source, dest)
end

"""
    getbranchname(::AbstractTree, branch)
    getbranchname(branch)

Returns the branch name associated with a branch from a tree. For some
branch types, it will be able to extract the branch name without reference to
the tree.
"""
function getbranchname end
getbranchname(tree::T, branch::Union{Int, B}) where
    {TT, RT, NL, N, B <: AbstractBranch{RT, NL},
     T <: AbstractTree{TT, RT, NL, N, B}} = _getbranchname(tree, branch)

"""
    hasrootheight(tree::AbstractTree)

"""
function hasrootheight end
hasrootheight(tree::AbstractTree{OneTree, <: Rooted}) = _hasrootheight(tree)
hasrootheight(tree::AbstractTree{OneTree, Unrooted}) = false
hasrootheight(tree::AbstractTree{ManyTrees}) =
    all(hasrootheight(tree) for tree in gettrees(tree))

"""
    getrootheight(tree::AbstractTree)

"""
function getrootheight(tree::AbstractTree{OneTree, <: Rooted})
    return _getrootheight(tree)
end

"""
    setrootheight!(tree::AbstractTree, height)

"""
function setrootheight!(tree::AbstractTree{OneTree, <: Rooted}, height)
    return _setrootheight!(tree, height)
end

# AbstractNode methods
"""
    isleaf(tree::AbstractTree, node)

Is the node (referenced by name or node object) a leaf of the tree?
"""
function isleaf end
isleaf(tree::AbstractTree{OneTree}, node) = _isleaf(tree, node)

"""
    isroot(tree::AbstractTree, node)

Is the node (referenced by name or node object) a root of the tree?
"""
function isroot end
isroot(tree::AbstractTree{OneTree}, node) = _isroot(tree, node)

"""
    isinternal(tree::AbstractTree, node)

Is the node (referenced by name or node object) internal to the tree (neither
root nor leaf)?
"""
function isinternal end
isinternal(tree::AbstractTree{OneTree}, node) = _isinternal(tree, node)

"""
    isunattached(tree::AbstractTree, node)

Is the node (referenced by name or node object) unattached (i.e. not connected
to other nodes)?
"""
function isunattached end
isunattached(tree::AbstractTree{OneTree}, node) = _isunattached(tree, node)

"""
    indegree(tree::AbstractTree, node)

Return in degree of node.
"""
function indegree end
indegree(tree::AbstractTree{OneTree}, node) = _indegree(tree, node)

"""
    outdegree(tree::AbstractTree, node)

Return out degree of node.
"""
function outdegree end
outdegree(tree::AbstractTree{OneTree}, node) = _outdegree(tree, node)

"""
    degree(tree::AbstractTree, node)

Return the degree of a node including all connections.
"""
function degree end
degree(tree::AbstractTree{OneTree}, node) = _degree(tree, node)

"""
    hasoutboundspace(tree::AbstractTree, node)

Does the node have space for an[other] outbound connection?
"""
function hasoutboundspace end
hasoutboundspace(tree::AbstractTree{OneTree}, node) =
    _hasoutboundspace(tree, node)

"""
    hasinbound(tree::AbstractTree, node)

Does the node have an inbound connection?
"""
function hasinbound end
hasinbound(tree::AbstractTree{OneTree}, node) = _hasinbound(tree, node)

"""
    hasinboundspace(tree::AbstractTree, node)

Does the node have space for an inbound connection?
"""
function hasinboundspace end
hasinboundspace(tree::AbstractTree{OneTree}, node) =
    _hasinboundspace(tree, node)

"""
    getinbound(tree::AbstractTree, node)

return the inbound branch to this node (returns name for node name, branch for
node).
"""
function getinbound end
getinbound(tree::AbstractTree{OneTree, <:Rooted}, node) =
    _getinbound(tree, node)

"""
    getparent(tree::AbstractTree, node)

Return [the name of] the parent node for this node [name]. Second method may
not be implemented for some node types.
"""
function getparent end
getparent(tree::AbstractTree{OneTree, <:Rooted}, node) =
    _getparent(tree, node)

"""
    getoutbounds(tree::AbstractTree, nodename)

Return the names of the outbound branches from this node.
"""
function getoutbounds end
getoutbounds(tree::AbstractTree{OneTree, <:Rooted}, node) =
    _getoutbounds(tree, node)

"""
    getchildren(tree::AbstractTree, node)

Return the [name(s) of] the child node(s) for this node [name].
"""
function getchildren end
getchildren(tree::AbstractTree{OneTree, <:Rooted}, node) =
    _getchildren(tree, node)

"""
    getconnections(tree::AbstractTree, nodee)

Returns all of the branches connected to a node.
"""
function getconnections end
getconnections(tree::AbstractTree{OneTree}, node) =
    _getconnections(tree, node)

"""
    getsiblings(tree::AbstractTree, node)

Returns all of the siblings of a node. Must be implemented for any unrooted
AbstractNode subtype, can be inferred from _getparent and _getchildren for
a rooted node.
"""
function getsiblings end
getsiblings(tree::AbstractTree{OneTree}, node) = _getsiblings(tree, node)

"""
    hasheight(tree::AbstractTree, node)

Does the node have a height defined?
"""
function hasheight end
hasheight(tree::AbstractTree, node) = false
function hasheight(tree::AbstractTree{OneTree, <:Rooted}, node)
    return _hasheight(tree, node) ||
        (_hasrootheight(tree) &&
         mapreduce(b -> haslength(tree, b), &, branchhistory(tree, node);
         init = _hasrootheight(tree)))
end

"""
    getheight(tree::AbstractTree, node)

Return the height of the node.
"""
function getheight(tree::AbstractTree{OneTree, <:Rooted}, node)
    return _hasheight(tree, node) ? _getheight(tree, node) :
        mapreduce(b -> getlength(tree, b), +, branchhistory(tree, node);
                  init = getrootheight(tree))
end

"""
    setheight!(tree::AbstractTree, nodename, height)

Set the height of the node.
"""
function setheight!(tree::AbstractTree{OneTree, <:Rooted}, nodename, height)
    return _setheight!(tree, nodename, height)
end

# AbstractBranch methods
"""
    src(tree::AbstractTree, branch)

Return the source node for this branch.
"""
function src end
src(tree::AbstractTree{OneTree, <:Rooted}, branch) = _src(tree, branch)

"""
    dst(tree::AbstractTree, branch)

Return the destination node for this branch.
"""
function dst end
dst(tree::AbstractTree{OneTree, <:Rooted}, branch) = _dst(tree, branch)

"""
    conns(tree::AbstractTree, branch)

Return the nodes connected to `branch`.
"""
function conns end
conns(tree::AbstractTree{OneTree}, branch) = _conns(tree, branch)

"""
    conn(tree::AbstractTree, branch, exclude)

Return the other node connected to `branch` that is not `exclude`.
"""
function conn end
conn(tree::AbstractTree{OneTree}, branch, exclude) =
    _conn(tree, branch, exclude)

"""
    getlength(tree::AbstractTree, branch)

Return the length of this branch.
"""
function getlength end
getlength(tree::AbstractTree{OneTree}, branch) = _getlength(tree, branch)

"""
    getleafnames(::AbstractTree)

Retrieve the leaf names from the tree.
"""
getleafnames(tree::AbstractTree) = _getleafnames(tree)

"""
    getleaves(::AbstractTree)

Retrieve the leaves from the tree.
"""
getleaves(tree::AbstractTree) = _getleaves(tree)

"""
    getleafinfo(::AbstractTree, label)

retrieve the leaf info for a leaf of the tree.
"""
getleafinfo(tree::AbstractTree) = _getleafinfo(tree)

"""
    setleafinfo!(::AbstractTree, table)

Set the leaf info for the leaves of the tree.
"""
setleafinfo!(tree::AbstractTree, table) = _setleafinfo!(tree, table)

"""
    getnoderecord(::AbstractTree, label)

retrieve the node info for a leaf of the tree.
"""
getnoderecord(tree::AbstractTree{OneTree}, label) = _getnoderecord(tree, label)

"""
    setnoderecord!(::AbstractTree, label, value)

Set the node info for a node of the tree.
"""
setnoderecord!(tree::AbstractTree{OneTree}, label, value) =
    _setnoderecord!(tree, label, value)

"""
    getbranchinfo(::AbstractTree, label)

retrieve the branch info for a leaf of the tree.
"""
getbranchinfo(tree::AbstractTree{OneTree}, label) = _getbranchinfo(tree, label)

"""
    setbranchinfo!(::AbstractTree, label, value)

Set the branch info for a branch of the tree.
"""
setbranchinfo!(tree::AbstractTree{OneTree}, label, value) =
    _setbranchinfo!(tree, label, value)

"""
    resetleaves!(::AbstractTree)

Reset the leaf records to the current leaves, deleting all leaf records.
"""
resetleaves!(tree::AbstractTree{OneTree}) = _resetleaves!(tree)

"""
    validate(tree::AbstractTree)

Validate the tree by making sure that it is connected up correctly.
"""
function validate(tree::T) where
    {TT, RT, NL, N, B, T <: AbstractTree{TT, RT, NL, N, B}}
    nodes = _getnodes(tree)
    nodenames = _getnodenames(tree)
    branches = _getbranches(tree)
    branchnames = _getbranchnames(tree)
    if !isempty(nodes) || !isempty(branches)
        # We need to validate the connections
        if Set(_getinbound(tree, node) for node in nodes
               if _hasinbound(tree, node)) != Set(branches)
            warn("Inbound branches must exactly match Branch labels")
            return false
        end

        if Set(mapreduce(node -> _getoutbounds(tree, node), append!,
                         nodes; init = B[])) != Set(branches)
            warn("Node outbound branches must exactly match Branch labels")
            return false
        end

        if !(Set(_src(tree, branch) for branch in branches) ⊆ Set(nodenames))
            warn("Branch sources must be node labels")
            return false
        end

        if !(Set(_dst(tree, branch) for branch in branches) ⊆ Set(nodenames))
            warn("Branch destinations must be node labels")
            return false
        end
    end

    return _validate(tree)
end

"""
    traversal(::AbstractTree, ::TraversalOrder)
    traversal(::AbstractTree, ::TraversalOrder, init)

Return an iterable object for a tree containing nodes in given order -
preorder, inorder, postorder or breadthfirst - optionally starting from init.
"""
function traversal end
traversal(tree::AbstractTree, order::TraversalOrder = preorder) =
    _traversal(tree, order)
traversal(tree::AbstractTree{TT, RT, NL, N},
          order::TraversalOrder, init::Union{N, NL}) where {TT, RT, NL, N} =
    _traversal(tree, order, [init])

"""
    getancestors(tree::AbstractTree, nodename)

Return the name of all of the nodes that are ancestral to this node.
"""
function getancestors(tree::AbstractTree{OneTree, <:Rooted}, nodename)
    return _treepast(tree, nodename)[2][2:end]
end

"""
    getdescendants(tree::AbstractTree, nodename)

Return the names of all of the nodes that descend from this node.
"""
function getdescendants(tree::AbstractTree{OneTree, <:Rooted}, nodename)
    return _treefuture(tree, nodename)[2][2:end]
end


text:
module MusicXML

# Alias
export MX
const MX = MusicXML

# Type Importer Macro
export @importMX
"""
    @importMX

Imports all the MusicXML types in case you are sure no conflicts happens.

# Just put it in the code:
```julia
@importMX # gets translated to: import MusicXML: ScorePartwise, Part, Measure, Note, Chord, Unpitched, Rest, Pitch, Grace, Attributes, Time, Transpose, Clef, Key, PartList, ScorePart, MidiInstrument, MidiDevice, ScoreInstrument
```
"""
macro importMX()
    out = quote
        import MusicXML: ScorePartwise, Part, Measure, Note, Chord, Unpitched, Rest, Pitch, Grace, Attributes, Time, Transpose, Clef, Key, PartList, ScorePart, MidiInstrument, MidiDevice, ScoreInstrument
    end
    return out
end

# Type Helper Macro
export @MX
"""
    @MX

A macro that adds `MX.` to the name of the MusicXML types.

Since MusicXML's types are not exported from the package to avoid conflicts with the similarly named types from other libraries (such as `Dates.Time`, `MIDI.Note`), `@MX` is given that facilitates using MusicXML types.

`@MX` assumes that all the types with these kinds of names inside in front of a macro, between `()`, or between `begin end` are MusicXML types unless explicitly written (e.g. `MIDI.Note`).

!!! warn
    It will replace the name of the types if they are used inside strings.

# Examples
```julia
@MX begin
    Note(pitch = Pitch(step = "C", alter = 0, octave = 4), duration =  4))
    # all the code ...
end
```
is translated to:
```julia
    MX.Note(pitch = Pitch(step = "C", alter = 0, octave = 4), duration =  4))
```

Type is untouched if it is called explicitly from another library:
```julia
@MX Dates.Time(12,53,40)
```
remains as:
```julia
Dates.Time(12,53,40)
```
"""
macro MX(expr)
    expr_string = string(expr)
    expr_string = reduce(replace,
        (
        r"(?<!\.)\bScorePartwise\b" => "MX.ScorePartwise",
        r"(?<!\.)\bPart\b" => "MX.Part",
        r"(?<!\.)\bMeasure\b" => "MX.Measure",
        r"(?<!\.)\bNote\b" => "MX.Note",
        r"(?<!\.)\bChord\b" => "MX.Chord",
        r"(?<!\.)\bUnpitched\b" => "MX.Unpitched",
        r"(?<!\.)\bRest\b" => "MX.Rest",
        r"(?<!\.)\bPitch\b" => "MX.Pitch",
        r"(?<!\.)\bAttributes\b" => "MX.Attributes",
        r"(?<!\.)\bTime\b" => "MX.Time",
        r"(?<!\.)\bTranspose\b" => "MX.Transpose",
        r"(?<!\.)\bClef\b" => "MX.Clef",
        r"(?<!\.)\bKey\b" => "MX.Key",
        r"(?<!\.)\bPartList\b" => "MX.PartList",
        r"(?<!\.)\bScorePart\b" => "MX.ScorePart",
        r"(?<!\.)\bMidiInstrument\b" => "MX.MidiInstrument",
        r"(?<!\.)\bMidiDevice\b" => "MX.MidiDevice",
        r"(?<!\.)\bScoreInstrument\b" => "MX.ScoreInstrument",
        ),
        init = expr_string)
    return esc(Meta.parse(expr_string))
end

using AcuteML

# I/O functions
include("io.jl")
# Types:
include("types.jl")

# Utilities
include("utilities.jl")

################################################################
# include("show.jl")

include("precompile_includer.jl")
end


text:
## context

## A context stores objects where named values are looked up.
mutable struct Context
    view ## of what? a Dict, Module, CompositeKind, DataFrame.parent.
    parent ## a context or nothing
    _cache::Dict
end
function Context(view, parent=nothing)
    d = Dict()
    d["."] = isa(view, AnIndex) ? view.value : view
    Context(view, parent, d)
end

## used by renderTokens
function ctx_push(ctx::Context, view)
    Context(view, ctx) ## add context as a parent
end

## Lookup value by key in the context
function lookup(ctx::Context, key)

    if haskey(ctx._cache, key)
        value = ctx._cache[key]
    else
        context = ctx
        value = nothing
        while value == nothing && context != nothing
            ## does name have a .?
            if occursin(r"\.", key)
                ## do something with "."
                ## we use .[ind] to refer to value in parent of given index;
                m = match(r"^\.\[(.*)\]$", key)
                m == nothing && error("Not implemented. Can use Composite Kinds in the view.")

                idx = m.captures[1]
                vals = context.parent.view

                if isa(vals, Vector)
                    if idx == "end"
                        value = AnIndex(:end, vals[end])
                    else
                        ind = Base.parse(Int, idx[1:end])
                        value = AnIndex(ind, vals[ind])
                    end
                    break
                end
            else
                ## strip leading, trailing whitespace in key
                value = lookup_in_view(context.view, stripWhitepace(key))
            end
            context = context.parent

        end

        ## cache
        ctx._cache[key] = value
    end

    if value === Function
        value = value()
    end

    return(value)
end

## Lookup value in an object by key
## This of course varies based on the view.
## we special case dataframes here, so that we don't have to assume package is loaded
function lookup_in_view(view, key)
    if is_dataframe(view)
        if occursin(r":", key)  key = key[2:end] end
        key = Symbol(key)
        out = nothing
        if haskey(view, key)
            out = view[1, key] ## first element only
        end
        return out
    else
        _lookup_in_view(view, key)
    end
end

function _lookup_in_view(view::Dict, key)

    ## is it a symbol?
    if occursin(r"^:", key)
        key = Symbol(key[2:end])
    end

    out = nothing
    if haskey(view, key)
        out = view[key]
    end

    out
end

function _lookup_in_view(view::Module, key)

    hasmatch = false
    re = Regex("^$key\$")
    for i in names(view, all=true)
        if occursin(re, string(i))
            hasmatch = true
            break
        end
    end

    out = nothing
    if hasmatch
        out = getfield(view, Symbol(key))  ## view.key
    end
    out

end

## Default is likely not great, but we use CompositeKind
function _lookup_in_view(view, key)
    nms = fieldnames(typeof(view))
    re = Regex(key)
    has_match = false
    for i in nms
        if occursin(Regex(key), string(i))
            has_match=true
            break
        end
    end

    out = nothing
    if has_match
        out = getfield(view, Symbol(key))  ## view.key
    end

    out
end


text:
@testset "LocalField" begin

  @testset "Creation" begin
    Qx, x = FlintQQ["x"]
    f = x^2-2*x+2
    K, a = local_field(f, 2, 10, "a", Hecke.EisensteinLocalField, check = false)
    @test precision(K) == 20
    @test characteristic(K) == 0
    @test prime(K) == 2

    Kt, t = PolynomialRing(K, "t")
    g = t^2+2
    L, b = local_field(g, "b", Hecke.EisensteinLocalField, check = false)
    @test precision(L) == 2*precision(g)
    @test precision(b) == precision(L)
    @test degree(L) == 2
    @test absolute_degree(L) == 4
    @test prime(L) == 2

    Q2 = PadicField(2, 10)
    Q2s, s = PolynomialRing(Q2, "s")
    f = s^2+s+1
    Ku, c = local_field(f, "s", Hecke.UnramifiedLocalField, check = false)
    @test precision(Ku) == precision(f)
    @test precision(f) == precision(c)
    @test degree(Ku) == 2
    @test absolute_degree(Ku) == 2
    @test inertia_degree(Ku) == 2
    @test ramification_index(Ku) == 1
    @test absolute_ramification_index(Ku) == 1
    @test absolute_inertia_degree(Ku) == 2

    Lu, u = PolynomialRing(L, "u")
    Lu, d = local_field(u^2+u+1, "s", Hecke.UnramifiedLocalField, check = false)
    @test absolute_degree(Lu) == 8
    @test ramification_index(Lu, K) == 2
    @test inertia_degree(Lu, K) == 2

  end

  @testset "Norm" begin
    K = QadicField(3, 4, 10)[1]
    Kx, x = PolynomialRing(K, "x")
    L = eisenstein_extension(x^20+3)[1]
    b = @inferred basis(L)
    for i = 1:10
      r = 1+2*uniformizer(L)^i * sum([rand(1:10)*b[i] for i in 1:5])
      M = @inferred representation_matrix(r)
      n = @inferred norm(r)
      @test n == det(M)
      @test trace(r) == trace(M)
    end
  end

  @testset "Completions" begin
    K, gK = cyclotomic_field(15)
    OK = maximal_order(K)
    lP = prime_decomposition(OK, 3)
    P = lP[1][1]
    Kp, mKp = @inferred Hecke.generic_completion(K, P, 100)
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test valuation(mKp(elem_in_nf(uniformizer(P)))) == fmpq(1, 2)
    @test valuation(mKp(gen(K)))*2 == valuation(gen(K), P)
    @test iszero(defining_polynomial(K)(mKp(gen(K))))

    P = prime_decomposition(OK, 7)[1][1]
    Kp, mKp = @inferred Hecke.unramified_completion(K, P, 100)
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test isone(valuation(mKp(elem_in_nf(uniformizer(P)))))
    @test valuation(mKp(gen(K))) == valuation(gen(K), P)
    @test iszero(defining_polynomial(K)(mKp(gen(K))))

    K, gK = cyclotomic_field(9)
    OK = maximal_order(K)
    lP = prime_decomposition(OK, 3)
    P = lP[1][1]
    Kp, mKp = @inferred Hecke.totally_ramified_completion(K, P, 100)
    @test valuation(mKp(elem_in_nf(uniformizer(P)))) == fmpq(1, 6)
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test valuation(mKp(gen(K)))*2 == valuation(gen(K), P)
    @test iszero(defining_polynomial(K)(mKp(gen(K))))

    K, gK = cyclotomic_field(55)
    OK = maximal_order(K)
    lP = prime_decomposition(OK, 11)
    P = lP[1][1]
    Kp, mKp = @inferred Hecke.generic_completion(K, P, 100)
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test valuation(mKp(elem_in_nf(uniformizer(P)))) == fmpq(1, 10)
    @test valuation(mKp(gen(K)))*10 == valuation(gen(K), P)
    @test iszero(defining_polynomial(K)(mKp(gen(K))))

    lP = prime_decomposition(OK, 7)
    P = lP[1][1]
    Kp, mKp = @inferred Hecke.unramified_completion(K, P, 100)
    @test isone(valuation(mKp(elem_in_nf(uniformizer(P)))))
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test valuation(mKp(gen(K))) == valuation(gen(K), P)
    @test iszero(defining_polynomial(K)(mKp(gen(K))))

    Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
    f = x^10 + 5*x^9 - 633*x^8 + 38157*x^7 + 1360601*x^6 - 39808345*x^5 - 464252491*x^4 - 17622187401*x^3 + 2826886632714*x^2 - 96335539805104*x + 1313743135204448
    K, a = number_field(f, check = false, cached = false)
    OK = maximal_order(K)
    P = prime_decomposition(OK, 2)[1][1]
    Kp, mKp = @inferred Hecke.unramified_completion(K, P, 100)
    @test isone(valuation(mKp(elem_in_nf(uniformizer(P)))))
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test valuation(mKp(gen(K))) == valuation(gen(K), P)
    @test iszero(defining_polynomial(K)(mKp(gen(K))))

    lP = prime_decomposition(OK, 11)
    if lP[1][2] == 5
      P = lP[1][1]
    else
      P = lP[2][1]
    end
    Kp, mKp = @inferred Hecke.totally_ramified_completion(K, P, 100)
    @test valuation(mKp(elem_in_nf(uniformizer(P)))) == fmpq(1, 5)
    @test isone(valuation(mKp\(uniformizer(Kp)), P))
    @test valuation(mKp(gen(K)))*5 == valuation(gen(K), P)
    #@test iszero(defining_polynomial(K)(mKp(gen(K))))

  end

  @testset "Exp and Log" begin
    K = PadicField(2, 100)
    Kx, x = PolynomialRing(K, "x", cached = false)
    L, b = Hecke.eisenstein_extension(x^7+2, "a")
    pi = uniformizer(L)
    @test iszero(log(pi))
    B = basis(L)
    for i = 15:20
      el = sum([rand(FlintZZ, 0:10)*B[j] for j = 1:7])*pi^i
      explog = exp(log(1+el))
      logexp = log(exp(el))
      @test iszero(setprecision(explog, precision(explog)-35) - 1 - el)
      @test iszero(setprecision(logexp, precision(logexp)-35)-el)
    end
  end

end


text:
const doc = """Example of program which uses [options] shortcut in pattern.

Usage:
  options_shortcut_example.jl [options] <port>

Options:
  -h --help                show this help message and exit
  --version                show version and exit
  -n, --number N           use N as a number
  -t, --timeout TIMEOUT    set timeout TIMEOUT seconds
  --apply                  apply changes to database
  -q                       operate in quiet mode

"""

using DocOpt

arguments = docopt(doc, version="1.0.0rc2")
dump(arguments)


text:
# Unary operations on fields and arrays of fields

for op in (:+,:-,:tr, :transpose, :adjoint, :symmetic_part)
  @eval begin

    function ($op)(f::Field)
      apply_kernel_to_field(bcast($op),f)
    end

    function ($op)(f::AbstractArray{<:Field})
      apply_to_field_array(bcast($op),f)
    end

  end

end

# Binary operations on fields and arrays of fields

for op in (:+,:-,:*,:inner,:outer)
  @eval begin

    function ($op)(f::Field, g::Field)
      apply_kernel_to_field(bcast($op),f,g)
    end

    function ($op)(f::AbstractArray{<:Field},g::AbstractArray{<:Field})
      apply_to_field_array(bcast($op),f,g)
    end

  end

end

# Define the gradient of some operations

@inline apply_kernel_gradient(k::BCasted{typeof(+)},a) = field_gradient(a)

@inline apply_kernel_gradient(k::BCasted{typeof(-)},a) = apply_kernel_to_field(k,field_gradient(a))

for op in (:+,:-)
  @eval begin
    @inline function apply_kernel_gradient(k::BCasted{typeof($op)},a,b)
      ga = field_gradient(a)
      gb = field_gradient(b)
      apply_kernel_to_field(k,ga,gb)
    end
    @inline function apply_kernel_gradient(k::BCasted{typeof($op)},f...)
      apply_kernel_to_field(k,field_gradients(f...)...)
    end
  end
end

for op in (:+,:-)
  @eval begin
    function apply_gradient(k::Valued{BCasted{typeof($op)}},f...)
      g = field_array_gradients(f...)
      apply(k,g...)
    end
  end
end


text:
using DataFrames
using Statistics
using StatsBase
using PrettyTables
using PrettyTables.Crayons

const latex_booktabs = LatexTableFormat(
    top_line       = "\\toprule",
    header_line    = "\\midrule",
    mid_line       = "\\midrule",
    bottom_line    = "\\bottomrule",
    left_vline     = "|",
    mid_vline      = "|",
    right_vline    = "|",
    header_envs    = ["textbf"],
    subheader_envs = ["texttt"]
)

"""
	print_table(df::DataFrame, metric_col=:tst_auc)

Prints dataframe with columns [:modelname, :dataset] and one scalar variable column
given by `metric_col` argument. By default highlights maximum value in each row.
Last row of the dataframe contains average rank of each model.
"""
function rank_table(df::DataFrame, metric_col=:tst_auc)
	# check if column names are present
	(!(String(metric_col) in names(df)) || !("modelname" in names(df)) || !("dataset" in names(df))) && error("Incorrect column names.")

	std_column = String(metric_col)*"_std"
	(!(std_column in names(df))) && error("DataFrame does not contain std for the given metric.")

	top_10_std_column = String(metric_col)*"_top_10_std"
	(!(top_10_std_column in names(df))) && error("DataFrame does not contain std for the given metric.")

	# get all the models that are present in the dataframe
	all_models = unique(df.modelname)

	# transposing the dataframe
	results = []
	for (dkey, dg) in pairs(groupby(df, :dataset))
		models = dg[:modelname]
		metric = dg[metric_col]
		row = merge(
			(dataset = dkey.dataset,),
			(;zip(Symbol.(models), metric)...))

		# when there are no evaluation files for some models on particular dataset
		if nrow(dg) < length(all_models)
			mm = setdiff(all_models, models)
			row = merge(row, (;zip(Symbol.(mm), fill(NaN, length(mm)))...))
		end
		push!(results, DataFrame([row]))
	end

	rt = vcat(results...)
	sort!(rt, :dataset)

	# add average std of the winning hyperparameter
	mean_std = [mean(
		filter(x -> (x.modelname == m), df)[std_column]) for m in names(rt[:,2:end])]
	push!(rt, ["MEAN_STD", mean_std...])

	# add top 10 std (more precisely mean over dataset)
	top_10_std = [mean(
		filter(x -> (x.modelname == m), df)[top_10_std_column]) for m in names(rt[:,2:end])]
	push!(rt, ["MEAN_TOP10_STD", top_10_std...])

	# round all values in the table before computing rank
	rt[:, 2:end] .= round.(rt[:, 2:end], digits=2)

	# add average rank
	mask_nan_max = (x) -> (isnan(x) ? -Inf : x)
	rs = zeros(size(rt, 2) - 1)
	for row in eachrow(rt[1:end-2,:])
		rs .+= StatsBase.competerank(mask_nan_max.(Vector(row[2:end])), rev = true)
	end
	rs ./= (size(rt, 1) - 2)
	rs = round.(rs, digits=1)
	push!(rt, ["RANK", rs...])

	rt
end

"""
	print_rank_table(rt::DataFrame; backend=:text)
	print_rank_table(io::IO, rt::DataFrame; backend=:text)

Pretty prints the rank table created by `rank_table` function either into given io or to stdout.
There are three backends to choose from `:txt (default)`, `:tex` and `:html`.
"""
print_rank_table(rt::DataFrame; backend=:txt) = print_rank_table(stdout, rt)

function print_rank_table(io::IO, rt::DataFrame; backend=:txt)
	mask_nan_max = (x) -> (isnan(x) ? -Inf : x)

	# horizontal lines to separate derived statistics from the rest of the table
	hlines = [size(rt, 1) - 3, size(rt, 1) - 2, size(rt, 1) - 1]

	# highlight maximum values of metric in each row (i.e. per dataset)
	f_hl_best = (data, i, j) -> (i < size(rt, 1)) && (data[i,j]  == maximum(mask_nan_max, rt[i, 2:end]))

	# highlight minimum rank in last row
	f_hl_best_rank = (data, i, j) -> i == size(rt, 1) && (data[i,j] == minimum(rt[i, 2:end]))

	# format metrics to 2 digits and rank to 1
	f_float = (v, i, j) -> (i == size(rt, 1)) ? ft_printf("%.1f")(v,i,j) : ft_printf("%.2f")(v,i,j)

	if backend == :html
		hl_best = HTMLHighlighter(f_hl_best, HTMLDecoration(color = "blue", font_weight = "bold"))
		hl_best_rank = HTMLHighlighter(f_hl_best_rank, HTMLDecoration(color = "red", font_weight = "bold"))

		pretty_table(
			io,	rt,
			formatters=f_float,
			highlighters=(hl_best, hl_best_rank),
			nosubheader=true,
			tf=html_minimalist
		)
	elseif backend == :tex
		hl_best = LatexHighlighter(f_hl_best, ["color{blue}","textbf"])
		hl_best_rank = LatexHighlighter(f_hl_best_rank,	["color{red}","textbf"])

		pretty_table(
			io,	rt,
			backend=:latex,
			formatters=f_float,
			highlighters=(hl_best, hl_best_rank),
			nosubheader=true,
			hlines=hlines,
			tf=latex_booktabs
		)
	else
		hl_best = Highlighter(f=f_hl_best, crayon=crayon"yellow bold")
		hl_best_rank = Highlighter(f=f_hl_best_rank, crayon=crayon"green bold")

		pretty_table(
			io, rt,
			formatters=f_float,
			highlighters=(hl_best, hl_best_rank),
			nosubheader=true,
			body_hlines=hlines
		)
	end
end


text:
################################################################################
#
#   linfactors.jl : functions for computing linear factors
#
################################################################################

export linearDlog

"""
    linearDlog{T <: PolyElem}(basis::T, degExt::Integer, h0::T, h1::T,
                              card::Nemo.fmpz, defPol::T)

Compute the discrete logarithm of the elements of type ``X + μ`` and of `h1`.

This algorithm works only if the provided basis is *linear*.
"""
function linearDlog{T <: PolyElem}(basis::T, degExt::Integer, h0::T, h1::T,
                                   card::Nemo.fmpz, defPol::T)

    # We set some constants, arrays, matrices
    F = base_ring(h0)
    charac::Int = characteristic(F)
    x = gen(F)
    X = gen(parent(h0))
    j = 0
    n = 0
    o = -coeff(basis, 0)
    ind::Int = coeff(o, 0) + coeff(o, 1)*charac + 1

    S = MatrixSpace(ZZ, charac^2+2,charac^3+charac+1)
    M = zero(S)
    Pq = pglCosets(x)

    # We iterate over Pq = PGL(P_1(F_q²))/PGL(P_1(F_q)) to create new equations
    # involving X and its translations X - μ with μ in F_q², we keep only the
    # one with a smooth left side and we fill a matrix to remember which
    # translations were used. Contrary to the descent algorithm, we also keep
    # trace of the linear factors in the left side.
    for m in Pq
        N = makeEquation(m, X, h0, h1)

        if isSmooth(N, 1)
            j += 1
            unit = fillMatrixBGJT!(M, j, m, F)
            fact = factor(N)
            for f in fact
                cst = -coeff(f[1], 0)
                index::Int = coeff(cst, 0) + coeff(cst, 1)*charac + 1
                M[index, j] -= f[2]
            end
            M[charac^2+2, j] = 1
        end
    end

    M = subMatrix(M, charac^2 + 2, j)
    M = transpose(M)

    # We compute the logarithm of the linear elements and h1
    # modulo the small factors using Pohling Hellman algorithm

    # We first precompute parameters that will not change during the q² calls to
    # Pohlig-Hellman algorithm, in order not to do the same computations many
    # times
    k, d, mid, B = pohligHellmanParam(card, basis, defPol)
    arr, inverses = pohligHellmanRoots(B, k, mid, defPol)

    dlogs = Array(Nemo.fmpz, charac^2+1)
    i = 0
    for y in F
        i += 1
        dlogs[i] = pohligHellmanPrec(card, mid, X-y, defPol, k , d, B, arr,
                                     inverses)[1]
    end

    i += 1
    dlogs[i], n = pohligHellmanPrec(card, basis, h1, defPol, k, d, B, arr,
                                    inverses)

    # We next look at the big factors of ``card-1``, we compute the result
    # modulo each factor and reconstruct the result using Chinese remaindering
    # theorem
    v = div(card-1, n)
    bigFactors = factor(v)

    # We compute the one-dimensionnal kernel, and we choose the vector with a 1
    # at index `ind`, because this is the index of our basis (we chose it linear
    # for that purpose)
    for p in bigFactors
        Mred = rrefMod(M, p[1])[1]
        ker = Nemo.fmpz[Mred[i, charac^2+2] for i in 1:charac^2]
        s = gcdinv(ker[ind], p[1])[2]
        for i in 1:charac^2
            dlogs[i] = crt(dlogs[i], n, ker[i]*s, p[1])
        end
        dlogs[charac^2+1] = crt(dlogs[charac^2+1], n, -s, p[1])
        n *= p[1]
    end

    # And we put the results in a dictionnary, to be more easy to use
    dico = Dict()
    i = 0
    for y in F
        i += 1
        dico[X-y] = dlogs[i]
    end
    dico[h1] = dlogs[i+1]

    return dico
end

function linearDlog{T <: PolyElem}(basis:: Nemo.RingElem, degExt::Integer,
                                   h0::T, h1::T, card::Integer, defPol::T)
    return linearDlog(basis, degExt, h0, h1, Nemo.fmpz(card), defPol)
end


text:
# utility methods

isinitialized(obj::Any) = isfilled(obj)

set_field!(obj::Any, fld::Symbol, val) = (setfield!(obj, fld, val); fillset(obj, fld); nothing)
@deprecate set_field(obj::Any, fld::Symbol, val) set_field!(obj, fld, val)

get_field(obj::Any, fld::Symbol) = isfilled(obj, fld) ? getfield(obj, fld) : error("uninitialized field $fld")

clear = fillunset

has_field(obj::Any, fld::Symbol) = isfilled(obj, fld)

function copy!{T}(to::T, from::T)
    fillunset(to)
    fill = filled(from)
    fnames = fld_names(T)
    for idx in 1:length(fnames)
        if fill[1, idx]
            name = fnames[idx]
            set_field!(to, name, getfield(from, name))
        end
    end
    nothing
end

function add_field!(obj::Any, fld::Symbol, val)
    typ = typeof(obj)
    attrib = meta(typ).symdict[fld]
    (attrib.occurrence != 2) && error("$(typ).$(fld) is not a repeating field")

    ptyp = attrib.ptyp
    jtyp = WIRETYPES[ptyp][4]
    (ptyp == :obj) && (jtyp = attrib.meta.jtype)

    !isdefined(obj, fld) && setfield!(obj, fld, jtyp[])
    push!(getfield(obj, fld), val)
    nothing
end
@deprecate add_field(obj::Any, fld::Symbol, val) add_field!(obj, fld, val)

protobuild{T}(::Type{T}, nv::Dict{Symbol}=Dict{Symbol,Any}()) = _protobuild(T(), collect(nv))

function _protobuild{T}(obj::T, nv)
    for (n,v) in nv
        fldtyp = fld_type(obj, n)
        set_field!(obj, n, isa(v, fldtyp) ? v : convert(fldtyp, v))
    end
    obj
end

# hash method that considers fill status of types
function protohash(v)
    h = 0
    fill = filled(v)
    fnames = fld_names(typeof(v))
    for idx in 1:length(fnames)
        fill[1, idx] && (h += hash(getfield(v, fnames[idx])))
    end
    hash(h)
end

# equality method that considers fill status of types
function protoeq{T}(v1::T, v2::T)
    fillv1 = filled(v1)
    fillv2 = filled(v2)
    fnames = fld_names(T)
    for idx in 1:length(fnames)
        (fillv1[1,idx] == fillv2[1,idx]) || (return false)
        if fillv1[1,idx]
            f = fnames[idx]
            (getfield(v1,f) == getfield(v2,f)) || (return false)
        end
    end
    true
end

# isequal method that considers fill status of types
function protoisequal{T}(v1::T, v2::T)
    fillv1 = filled(v1)
    fillv2 = filled(v2)
    fnames = fld_names(T)
    for idx in 1:length(fnames)
        (fillv1[1,idx] == fillv2[1,idx]) || (return false)
        if fillv1[1,idx]
            f = fnames[idx]
            isequal(getfield(v1,f), getfield(v2,f)) || (return false)
        end
    end
    true
end

function enumstr(enumname, t::Int32)
    for name in fld_names(typeof(enumname))
        (getfield(enumname, name) == t) && (return string(name))
    end
    error(string("Invalid enum value ", t, " for ", typeof(enumname)))
end


text:
function SpacecraftState(arg0::AbsolutePVCoordinates)
    return SpacecraftState((AbsolutePVCoordinates,), arg0)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::Attitude)
    return SpacecraftState((AbsolutePVCoordinates, Attitude), arg0, arg1)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::Attitude, arg2::Map)
    return SpacecraftState((AbsolutePVCoordinates, Attitude, Map), arg0, arg1, arg2)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::Attitude, arg2::jdouble)
    return SpacecraftState((AbsolutePVCoordinates, Attitude, jdouble), arg0, arg1, arg2)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::Attitude, arg2::jdouble, arg3::Map)
    return SpacecraftState((AbsolutePVCoordinates, Attitude, jdouble, Map), arg0, arg1, arg2, arg3)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::Map)
    return SpacecraftState((AbsolutePVCoordinates, Map), arg0, arg1)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::jdouble)
    return SpacecraftState((AbsolutePVCoordinates, jdouble), arg0, arg1)
end

function SpacecraftState(arg0::AbsolutePVCoordinates, arg1::jdouble, arg2::Map)
    return SpacecraftState((AbsolutePVCoordinates, jdouble, Map), arg0, arg1, arg2)
end

function SpacecraftState(arg0::Orbit)
    return SpacecraftState((Orbit,), arg0)
end

function SpacecraftState(arg0::Orbit, arg1::Attitude)
    return SpacecraftState((Orbit, Attitude), arg0, arg1)
end

function SpacecraftState(arg0::Orbit, arg1::Attitude, arg2::Map)
    return SpacecraftState((Orbit, Attitude, Map), arg0, arg1, arg2)
end

function SpacecraftState(arg0::Orbit, arg1::Attitude, arg2::jdouble)
    return SpacecraftState((Orbit, Attitude, jdouble), arg0, arg1, arg2)
end

function SpacecraftState(arg0::Orbit, arg1::Attitude, arg2::jdouble, arg3::Map)
    return SpacecraftState((Orbit, Attitude, jdouble, Map), arg0, arg1, arg2, arg3)
end

function SpacecraftState(arg0::Orbit, arg1::Map)
    return SpacecraftState((Orbit, Map), arg0, arg1)
end

function SpacecraftState(arg0::Orbit, arg1::jdouble)
    return SpacecraftState((Orbit, jdouble), arg0, arg1)
end

function SpacecraftState(arg0::Orbit, arg1::jdouble, arg2::Map)
    return SpacecraftState((Orbit, jdouble, Map), arg0, arg1, arg2)
end

function add_additional_state(obj::SpacecraftState, arg0::JString, arg1::Vector{jdouble})
    return jcall(obj, "addAdditionalState", SpacecraftState, (JString, Vector{jdouble}), arg0, arg1)
end

function ensure_compatible_additional_states(obj::SpacecraftState, arg0::SpacecraftState)
    return jcall(obj, "ensureCompatibleAdditionalStates", void, (SpacecraftState,), arg0)
end

function get_a(obj::SpacecraftState)
    return jcall(obj, "getA", jdouble, ())
end

function get_abs_pva(obj::SpacecraftState)
    return jcall(obj, "getAbsPVA", AbsolutePVCoordinates, ())
end

function get_additional_state(obj::SpacecraftState, arg0::JString)
    return jcall(obj, "getAdditionalState", Vector{jdouble}, (JString,), arg0)
end

function get_additional_states(obj::SpacecraftState)
    return jcall(obj, "getAdditionalStates", Map, ())
end

function get_attitude(obj::SpacecraftState)
    return jcall(obj, "getAttitude", Attitude, ())
end

function get_date(obj::SpacecraftState)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::SpacecraftState)
    return jcall(obj, "getE", jdouble, ())
end

function get_equinoctial_ex(obj::SpacecraftState)
    return jcall(obj, "getEquinoctialEx", jdouble, ())
end

function get_equinoctial_ey(obj::SpacecraftState)
    return jcall(obj, "getEquinoctialEy", jdouble, ())
end

function get_frame(obj::SpacecraftState)
    return jcall(obj, "getFrame", Frame, ())
end

function get_hx(obj::SpacecraftState)
    return jcall(obj, "getHx", jdouble, ())
end

function get_hy(obj::SpacecraftState)
    return jcall(obj, "getHy", jdouble, ())
end

function get_i(obj::SpacecraftState)
    return jcall(obj, "getI", jdouble, ())
end

function get_keplerian_mean_motion(obj::SpacecraftState)
    return jcall(obj, "getKeplerianMeanMotion", jdouble, ())
end

function get_keplerian_period(obj::SpacecraftState)
    return jcall(obj, "getKeplerianPeriod", jdouble, ())
end

function get_le(obj::SpacecraftState)
    return jcall(obj, "getLE", jdouble, ())
end

function get_lm(obj::SpacecraftState)
    return jcall(obj, "getLM", jdouble, ())
end

function get_lv(obj::SpacecraftState)
    return jcall(obj, "getLv", jdouble, ())
end

function get_mass(obj::SpacecraftState)
    return jcall(obj, "getMass", jdouble, ())
end

function get_mu(obj::SpacecraftState)
    return jcall(obj, "getMu", jdouble, ())
end

function get_orbit(obj::SpacecraftState)
    return jcall(obj, "getOrbit", Orbit, ())
end

function get_pv_coordinates(obj::SpacecraftState)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, ())
end

function get_pv_coordinates(obj::SpacecraftState, arg0::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (Frame,), arg0)
end

function has_additional_state(obj::SpacecraftState, arg0::JString)
    return jcall(obj, "hasAdditionalState", jboolean, (JString,), arg0)
end

function interpolate(obj::SpacecraftState, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", SpacecraftState, (AbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Collection), arg0, arg1)
end

function is_orbit_defined(obj::SpacecraftState)
    return jcall(obj, "isOrbitDefined", jboolean, ())
end

function shifted_by(obj::SpacecraftState, arg0::jdouble)
    return jcall(obj, "shiftedBy", SpacecraftState, (jdouble,), arg0)
end

function to_string(obj::SpacecraftState)
    return jcall(obj, "toString", JString, ())
end

function to_transform(obj::SpacecraftState)
    return jcall(obj, "toTransform", Transform, ())
end


text:
using Documenter, BiobakeryUtils

makedocs(
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "BiobakeryUtils.jl",
    pages = [
        "BiobakeryUtils" => "index.md"
        "General Utilities" => "general.md"
        "Working with HUMAnN2" => "humann2.md"
        "Working with MetaPhlAn2" => "metaphlan2.md"
        "Contributing" => "contributing.md"
    ],
    authors = "Kevin Bonham, PhD"
)

deploydocs(
    repo = "github.com/BioJulia/BiobakeryUtils.jl.git",
    push_preview=true,
    devbranch="main"
)


text:
# Adaptive mel-cepstrum analysis

function acep!(c::StridedVector{Cdouble}, x::Cdouble;
               λ::Float64=0.98,
               step::Float64=0.1,
               τ::Float64=0.9,
               pd::Int=4,
               eps::Float64=1.0e-6)
    assert_pade(pd)
    order = length(c) - 1
    prederr = ccall((:acep, libSPTK), Cdouble,
                    (Cdouble, Ptr{Cdouble}, Cint, Cdouble, Cdouble, Cdouble,
                     Cint, Cdouble),
                    x, c, order, λ, step, τ, pd, eps)
    prederr
end

# γ = -1/stage (-1 <= γ < 0)
function agcep!(c::StridedVector{Cdouble}, x::Cdouble, stage=1;
                λ::Float64=0.98,
                step::Float64=0.1,
                τ::Float64=0.9,
                eps::Float64=1.0e-6)
    stage >= 1 || throw(ArgumentError("stage >= 1 (-1 <= γ < 0)"))
    order = length(c) - 1
    prederr = ccall((:agcep, libSPTK), Cdouble,
                    (Cdouble, Ptr{Cdouble}, Cint, Cint, Cdouble, Cdouble,
                     Cdouble, Cdouble),
                    x, c, order, stage, λ, step, τ, eps)
    prederr
end

function amcep!(b::StridedVector{Cdouble}, x::Cdouble, α=0.41;
                λ::Float64=0.98,
                step::Float64=0.1,
                τ::Float64=0.9,
                pd::Int=4,
                eps::Float64=1.0e-6)
    assert_pade(pd)
    order = length(b) - 1
    prederr = ccall((:amcep, libSPTK), Cdouble,
                    (Cdouble, Ptr{Cdouble}, Cint, Cdouble, Cdouble, Cdouble,
                     Cdouble, Cint, Cdouble),
                    x, b, order, α, λ, step, τ, pd, eps)
    prederr
end


text:
module Tracing #end

export tracetasks, @tracetasks

using TaskHeritableStorage

include("ThreadSafeAccumulators.jl")
using .ThreadSafeAccumulators

# Key for Task Local Storage
struct TraceTLSKey
    id::Int  # Use a unique TLS key for each trace to allow composable tracing
    TraceTLSKey() = new(traceid[] += 1)
end
traceid = Threads.Atomic{Int}(0)
struct TraceTLSValue
    task_edges::Accumulator{Pair{Task,Task}}
    parent_task::Task  # Set to current_task()
    TraceTLSValue() = new(Accumulator([(current_task() => current_task())]), current_task())
    TraceTLSValue(edges, task) = new(edges, task)
end

function TaskHeritableStorage.fork_task_local_storage(k, v::TraceTLSValue)
    # It is thread-safe to push! to an Accumulator
    push!(v.task_edges, v.parent_task => current_task())
    # Fork a new "current task", but share the same edges accumulator:
    (k, TraceTLSValue(v.task_edges, current_task()))
end

"""
    @tracetasks expr

Equivalent to tracetasks() do ; expr; end
"""
macro tracetasks(e)
    :(tracetasks(()->$(esc(e))))
end

"""
    @tracetasks expr
    tracetasks() do ; end
    tracetasks(f::Function)

Run a function and trace all tasks that are spawned during its execution.
Returns an Array of pairs from `parent_task => child_task`. The root task is represented
with a special self-edge `t => t`.
"""
function tracetasks(f::Function)
    tls = task_local_storage()
    key = TraceTLSKey()
    task_local_storage(key, TraceTLSValue()) do
        f()
        return collect(task_local_storage()[key].task_edges)
    end
end

end  # module


text:
## ------------------------------------------------------------------
_set_ipass!(pass::Int) = setindex!(GLOB_STATE, pass, :ipass)
get_ipass() = get!(GLOB_STATE, :ipass, 0)
is_ipass0() = (get_ipass() == 0)
is_ipass1() = (get_ipass() == 1)
is_ipass2() = (get_ipass() == 2)

_keepout_git(path) = (basename(path) == ".git")

_rbfile_suffix() = ".rb.jl"
_is_rbfile(path) = endswith(path, _rbfile_suffix())

_include_order_configkey() = "include-order"
function include_order(bookdir)
    config = read_configfile(bookdir)
    order = get(config, _include_order_configkey(), "")
    order = (order isa String) ? [order] : order
    return bookdir_relpath.(order)
end

const _INCLUDE_MTIME_EVENT = FileMTimeEvent()
function _include_rdfile(path; force = false)
    !isfile(path) && return false
    !_is_rbfile(path) && return false

    need_update = has_event!(_INCLUDE_MTIME_EVENT, path)
    if force || need_update
        # include into main
        Base.include(Main, path)
        return true
    end
    return false
end

function _include_rbfiles!(book::RBook; force = false, nipass = 3)

    # TODO: for now including all files every time
    # A no-action-required mechanism must be provided
    force = true

    if is_ipass0()
        try
            empty!(book)
            bdir = bookdir(book)

            for ipass in 1:nipass

                _set_ipass!(ipass)

                # include in order
                for path in include_order(bdir)
                    _include_rdfile(path; force)
                end

                # walkdown (un-ordered)
                keepout = _keepout_git
                walkdown(bdir; keepout) do path
                    _include_rdfile(path; force)
                    return nothing
                end

            end

        finally
            _set_ipass!(0)
        end
    end

end


text:
using TreeKnit.SRG
using TreeTools
using Test

# Lots of intricated singletons

nwk1 = "(((((A:2,X3)AX3:1,X4)AX4:2,X5)AX5:2,B)AB:0.,X1,X2):0.R"
nwk2 = "((((A:1,X2)AX2:3,X1)AX1:3,B)AB:0.,X3,X4,X5):0.R"

t1 = parse_newick_string(nwk1);
t2 = parse_newick_string(nwk2);

MCCs = [["A","B"],["X1"],["X2"],["X3"],["X4"],["X5"]]

X1, X2 = SRG.shared_nodes(t1, t2, MCCs)

tc1, tc2, Xc1, Xc2 = begin
       tc1 = copy(t1)
       tc2 = copy(t2)
       Xc1 = copy(X1)
       Xc2 = copy(X2)
       SRG.fix_shared_singletons!(tc1, tc2, Xc1, Xc2, MCCs)
       tc1, tc2, Xc1, Xc2
end

@testset "fix_shared_singletons!" begin
	@test Xc1["AX3"][1] == :shared
	@test Xc1["AX3"][2][1:9] == "Singleton"
	@test Xc2[Xc1["AX3"][2]][2] == "AX3"

	@test Xc1["AX4"][1] == :shared
	@test Xc1["AX4"][2][1:9] == "Singleton"
	@test Xc2[Xc1["AX4"][2]][2] == "AX4"

	@test Xc1["AX5"][1] == :shared
	@test Xc1["AX5"][2][1:9] == "Singleton"
	@test Xc2[Xc1["AX5"][2]][2] == "AX5"
	#
	@test Xc2["AX1"][1] == :shared
	@test Xc2["AX1"][2][1:9] == "Singleton"
	@test Xc1[Xc2["AX1"][2]][2] == "AX1"

	@test Xc2["AX2"][1] == :shared
	@test Xc2["AX2"][2][1:9] == "Singleton"
	@test Xc1[Xc2["AX2"][2]][2] == "AX2"

	@test divtime(t1.lnodes["A"], t1.root) == divtime(tc1.lnodes["A"], tc1.root)
	@test divtime(t1.lnodes["AX3"], t1.root) == divtime(tc1.lnodes["AX3"], tc1.root)
	@test divtime(t1.lnodes["AX4"], t1.root) == divtime(tc1.lnodes["AX4"], tc1.root)
	@test divtime(t1.lnodes["AX5"], t1.root) == divtime(tc1.lnodes["AX5"], tc1.root)

	@test divtime(t2.lnodes["A"], t2.root) == divtime(tc2.lnodes["A"], tc2.root)
	@test divtime(t2.lnodes["AX2"], t2.root) == divtime(tc2.lnodes["AX2"], tc2.root)
	@test divtime(t2.lnodes["AX1"], t2.root) == divtime(tc2.lnodes["AX1"], tc2.root)
end

# Weird branch lenghts
nwk1 = "((A:1,B:1)AB:1,C:2)R:0"
nwk2 = "(A:4,(B:0,C:0):4)R:0" # from a resolution

t1 = parse_newick_string(nwk1);
t2 = parse_newick_string(nwk2);

MCCs = [["A"], ["B", "C"]]

X1, X2 = SRG.shared_nodes(t1, t2, MCCs)

tc1, tc2, Xc1, Xc2 = begin
       tc1 = copy(t1)
       tc2 = copy(t2)
       Xc1 = copy(X1)
       Xc2 = copy(X2)
       SRG.fix_shared_singletons!(tc1, tc2, Xc1, Xc2, MCCs)
       tc1, tc2, Xc1, Xc2
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
module TimeParallel

using Distributed
using ArrowMacros
using LinearAlgebra
using NSDEBase
using RungeKutta
using RecipesBase

include("abstract.jl")
include("weights.jl")
include("error.jl")
include("parareal/parareal.jl")
include("parareal/cache.jl")
include("parareal/iterate.jl")
include("parareal/solution.jl")
include("parareal/coarse.jl")
include("parareal/serial.jl")
include("parareal/distributed.jl")
include("parareal/mpi.jl")
include("parareal/solve.jl")
include("solve.jl")
include("plotrecipes.jl")

export AbstractTimeParallelSolver
export AbstractTimeParallelSolution
export Parareal
export PararealSolution
export ErrorControl, ErrorWeights
export ψ₁, ψ₂

end


text:
"""
	text(cmd0::String="", arg1=[]; kwargs...)

Plots text strings of variable size, font type, and orientation. Various map projections are
provided, with the option to draw and annotate the map boundaries.

Full option list at [`pstext`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- **A** : **azimuths** : -- Bool or [] --

    Angles are given as azimuths; convert them to directions using the current projection.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#a)
- **C** : **clearance** : -- Str --

    Sets the clearance between the text and the surrounding box [15%].
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#c)
- **D** : **offset** : -- Str --

    Offsets the text from the projected (x,y) point by dx,dy [0/0].
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#d)
- **F** : **text_attrib** : -- Str or number --

    Specify up to three text attributes (font, angle, and justification).
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#f)
- **G** : **fill** : -- Number or Str --

    Sets the shade or color used for filling the text box [Default is no fill].
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#g)
- $(GMT.opt_Jz)
- **L** : **list** : -- Bool or [] --

    Lists the font-numbers and font-names available, then exits.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#l)
- **N** : **no_clip** : --- Str or [] --

    Do NOT clip text at map boundaries.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#n)
- $(GMT.opt_P)
- **Q** : **change_case** : --- Str --

    Change all text to either lower or upper case.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#q)
- **T** : **text_box** : --- Str --

    Specify the shape of the textbox when using G and/or W.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#t)
- **W** : **line_attribs** : -- Str --

    Sets the pen used to draw a rectangle around the text string.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#w)
- **Z** : **threeD** : -- Str --

    For 3-D projections: expect each item to have its own level given in the 3rd column.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/pstext.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function text(cmd0::String="", arg1=[]; caller=[], K=false, O=false, first=true, kwargs...)

	arg2 = []		# May be needed if GMTcpt type is sent in via G
	N_args = isempty_(arg1) ? 0 : 1

	length(kwargs) == 0 && return monolitic("pstext", cmd0, arg1)	# Speedy mode

    d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

    cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", caller, O, " -JX12c/0")
	cmd, = parse_JZ(cmd, d)
	cmd  = parse_UVXY(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd = parse_params(cmd, d)

    cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, opt_i = read_data(d, cmd0, cmd, arg1, opt_R, "", opt_bi, opt_di)

	cmd, arg1, arg2, N_args = add_opt_cpt(d, cmd, [:C :color], 'C', N_args, arg1, arg2)

	cmd = add_opt(cmd, 'A', d, [:A :horizontal])
	cmd = add_opt(cmd, 'D', d, [:D :annot :annotate])
	cmd = add_opt(cmd, 'F', d, [:F :text_attrib])
    cmd = add_opt(cmd, 'G', d, [:G :fill])
	cmd = add_opt(cmd, 'I', d, [:I :inquire])
	cmd = add_opt(cmd, 'L', d, [:L :pen])
	cmd = add_opt(cmd, 'N', d, [:N :no_clip])
	cmd = add_opt(cmd, 'Q', d, [:Q :change_case])
	cmd = add_opt(cmd, 'T', d, [:T :text_box])
	cmd = add_opt(cmd, 'Z', d, [:Z :threeD])

	opt_W = ""
	pen = build_pen(d)						# Either a full pen string or empty ("")
	if (!isempty(pen))
		opt_W = " -W" * pen
	else
		for sym in [:W :line_attrib]
			if (haskey(d, sym))
				if (isa(d[sym], String))
					opt_W = " -W" * arg2str(d[sym])
				elseif (isa(d[sym], Tuple))	# Like this it can hold the pen, not extended atts
					opt_W = " -W" * parse_pen(d[sym])
				else
					error("Nonsense in W option")
				end
				break
			end
		end
	end

	if (!isempty(opt_W)) 		# We have a rectangle request
		cmd = finish_PS(d, cmd * opt_W, output, K, O)
	else
		cmd = finish_PS(d, cmd, output, K, O)
	end

    return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "pstext", arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
text!(cmd0::String="", arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
    text(cmd0, arg1; caller=caller, K=K, O=O,  first=false, kw...)

text(arg1=[]; caller=[], K=false, O=false, first=true, kw...) =
    text("", arg1; caller=caller, K=K, O=O, first=first, kw...)

text!(arg1=[]; caller=[],  K=true, O=true, first=false, kw...) =
    text("", arg1; caller=caller, K=K, O=O, first=first, kw...)

pstext  = text			# Alias
pstext! = text!			# Alias


text:
#Linear Kalman filter with control input example
#This example closely follows that given on "Greg Czerniak's Website". Namely
#the canonball example on this page: http://greg.czerniak.info/guides/kalman1/
#Basically we shoot a canonball into the air at a given angle of elevation.
#Using observations of the canonball's velocity and position we're going to
#track the position of the canonball.

#To use the Kalman Filter we need three things:
#1) Some observations
#2) The parameters of the Kalman Filter model describing the system process and
#the observation model. In this example we also need the control input
#parameters.
#3) An initial guess of the state.
#

#Let's import the modules required to execute the Kalman Filter
using StateSpace
using Distributions
using Gadfly
using Colors

################################################################################
#Section: Generate noisy Observations
#-------------------------------------------------------------------------------
#First we're going to provide some parameter values required to model the motion
#of the canonball.

#Set the Parameters
elevation_angle = 45.0 #Angle (measured anti-clockwise) between the ground (zero degrees) and the direction the canon ball is fired
muzzle_speed = 100.0 #Speed at which the canonball leaves the muzzle
initial_velocity = [muzzle_speed*cos(deg2rad(elevation_angle)), muzzle_speed*sin(deg2rad(elevation_angle))] #initial x and y components of the velocity
gravAcc = 9.81 #gravitational acceleration
initial_location = [0.0, 0.0] # initial position of the canonball
Δt = 0.1 #time between each measurement

#Functions describing the position of canonball
x_pos(x0::Float64, Vx::Float64, t::Float64) = x0 + Vx*t
y_pos(y0::Float64, Vy::Float64, t::Float64, g::Float64) = y0 + Vy*t - (g * t^2)/2
#Function to describe the evolution of the velocity in the vertical direction
velocityY(Vy::Float64, t::Float64, g::Float64) = Vy - g * t

#Give variances of the observation noise for the position and velocity
x_pos_var = 200.0
y_pos_var = 200.0
Vx_var = 1.0
Vy_var = 1.0

#Set the number of observations and preallocate vectors to store true and noisy
#measurement values
numObs = 145
x_pos_true = Vector{Float64}(numObs)
x_pos_obs = Vector{Float64}(numObs)
y_pos_true = Vector{Float64}(numObs)
y_pos_obs = Vector{Float64}(numObs)

Vx_true = Vector{Float64}(numObs)
Vx_obs = Vector{Float64}(numObs)
Vy_true = Vector{Float64}(numObs)
Vy_obs = Vector{Float64}(numObs)

#Generate the data (true values and noisy observations)
for i in 1:numObs
    x_pos_true[i] = x_pos(initial_location[1], initial_velocity[1], (i-1)*Δt)
    y_pos_true[i] = y_pos(initial_location[2], initial_velocity[2], (i-1)*Δt, gravAcc)
    Vx_true[i] = initial_velocity[1]
    Vy_true[i] = velocityY(initial_velocity[2], (i-1)*Δt, gravAcc)

    x_pos_obs[i] = x_pos_true[i] + randn() * sqrt(x_pos_var)
    y_pos_obs[i] = y_pos_true[i] + randn() * sqrt(y_pos_var)
    Vx_obs[i] = Vx_true[i] + randn() * sqrt(Vx_var)
    Vy_obs[i] = Vy_true[i] + randn() * sqrt(Vy_var)
end
#Create the observations vector for the Kalman filter
observations = [x_pos_obs Vx_obs y_pos_obs Vy_obs]'
#End Section: Generate noisy Observations
################################################################################

################################################################################
#Section: Describe Kalman Filter parameters
#-------------------------------------------------------------------------------

#Describe the system parameters
process_matrix = [[1.0, Δt, 0.0, 0.0] [0.0, 1.0, 0.0, 0.0] [0.0, 0.0, 1.0, Δt] [0.0, 0.0, 0.0, 1.0]]'
process_covariance = 0.01*eye(4)
observation_matrix = eye(4)
observation_covariance = 0.2*eye(4)
control_matrix = [[0.0, 0.0, 0.0, 0.0] [0.0, 0.0, 0.0, 0.0] [0.0, 0.0, 1.0, 0.0] [0.0, 0.0, 0.0, 1.0]]
control_input = [0.0, 0.0, -(gravAcc * Δt^2)/2, -(gravAcc * Δt)]

#Create an instance of the LKF with the control inputs
linCISMM = LinearGaussianCISSM(process_matrix, process_covariance, observation_matrix, observation_covariance, control_matrix, control_input)
#End Section: Describe Kalman Filter parameters
################################################################################

################################################################################
#Section: Set Initial Guess
#-------------------------------------------------------------------------------
initial_guess_state = [0.0, initial_velocity[1], 500.0, initial_velocity[2]]
initial_guess_covariance = eye(4)
initial_guess = MvNormal(initial_guess_state, initial_guess_covariance)
#End Section: Set Initial Guess
################################################################################

################################################################################
#Section: Execute Kalman Filter
#-------------------------------------------------------------------------------
filtered_state = filter(linCISMM, observations, initial_guess)
#End Section: Execute Kalman Filter
################################################################################

################################################################################
#Section: Plot Filtered results
#-------------------------------------------------------------------------------
#Here we are plotting the filtered results with Gadfly. See the Gadfly
#documentation for information about how plotting works if you are unfamiliar.
#Website: http://gadflyjl.org/
#Plot results
x_filt = Vector{Float64}(numObs)
y_filt = Vector{Float64}(numObs)
for i in 1:numObs
    current_state = filtered_state.state[i]
    x_filt[i] = current_state.μ[1]
    y_filt[i] = current_state.μ[3]
end

n = 3
getColors = distinguishable_colors(n, Color[LCHab(70, 60, 240)],
                                   transform=c -> deuteranopic(c, 0.5),
                                   lchoices=Float64[65, 70, 75, 80],
                                   cchoices=Float64[0, 50, 60, 70],
                                   hchoices=linspace(0, 330, 24))

cannonball_plot = plot(
    layer(x=x_pos_true, y=y_pos_true, Geom.line, Theme(default_color=getColors[3])),
    layer(x=[initial_guess_state[1]; x_filt], y=[initial_guess_state[3]; y_filt], Geom.line, Theme(default_color=getColors[1])),
    layer(x=x_pos_obs, y=y_pos_obs, Geom.point, Theme(default_color=getColors[2])),
    Guide.xlabel("X position"), Guide.ylabel("Y position"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Measurement of a Canonball in Flight")
    )
#End Section: Plot Filtered results
################################################################################


text:
function test_𝐼x(actual::Field1D, test::Field1D, input::Field1D)
    𝐼x!(test, input)
    return isapprox(test.data, actual.data; rtol = rtol_interpolation)
end

function test_𝐼x(actual::Field2D, test::Field2D, input::Field2D)
    𝐼x!(test, input)
    return isapprox(test.data, actual.data; rtol = rtol_interpolation)
end

function test_𝐼y(actual::Field2D, test::Field2D, input::Field2D)
    𝐼y!(test, input)
    return isapprox(test.data, actual.data; rtol = rtol_interpolation)
end

function test_∂x(actual::Field1D, test::Field1D, input::Field1D)
    ∂x!(test, input)
    return isapprox(test.data, actual.data; rtol = rtol_derivatives)
end

function test_∂x(actual::Field2D, test::Field2D, input::Field2D)
    ∂x!(test, input)
    return isapprox(test.data, actual.data, rtol = rtol_derivatives)
end

function test_∂y(actual::Field2D, test::Field2D, input::Field2D)
    ∂y!(test, input)
    return isapprox(test.data, actual.data, rtol = rtol_derivatives)
end


text:
export GaussianCopulaCSObs, GaussianCopulaCSModel

"""
    GaussianCopulaCSObs
GaussianCopulaCSObs(y, X, V)
A realization of Gaussian copula variance component data.
"""
struct GaussianCopulaCSObs{T <: BlasReal}
    # data
    n::Int
    p::Int
    y::Vector{T}
    X::Matrix{T}
    V::SymmetricToeplitz{T}
    vec::Vector{T}
    # working arrays
    ∇CSV::SymmetricToeplitz{T}
    ∇β::Vector{T}   # gradient wrt β
    ∇τ::Vector{T}   # gradient wrt τ
    ∇ρ::Vector{T}
    ∇σ2::Vector{T}
    Hβ::Matrix{T}   # Hessian wrt β
    Hτ::Matrix{T}   # Hessian wrt τ
    Hρ::Matrix{T}   # Hessian wrt ρ
    Hσ2::Matrix{T}   # Hessian wrt ρ
    Hρσ2::Matrix{T}  # Hessian wrt ρ, σ2
    Hβσ2::Vector{T}  # Hessian wrt β and σ2
    # Hβρ::Vector{T}
    res::Vector{T}  # residual vector res_i
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    xtx::Matrix{T}  # Xi'Xi
    storage_n::Vector{T}
    storage_p1::Vector{T}
    storage_np::Matrix{T}
    storage_pp::Matrix{T}
    added_term_numerator::Matrix{T}
    added_term2::Matrix{T}
end

function GaussianCopulaCSObs(
    y::Vector{T},
    X::Matrix{T}
    ) where T <: BlasReal
    n, p = size(X, 1), size(X, 2)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    V = SymmetricToeplitz(ones(T, n))
    vec = Vector{T}(undef, n)
    ∇CSV = SymmetricToeplitz(ones(T, n))
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇ρ  = Vector{T}(undef, 1)
    ∇σ2  = Vector{T}(undef, 1)
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    Hρ  = Matrix{T}(undef, 1, 1)
    Hσ2  = Matrix{T}(undef, 1, 1)
    Hρσ2 = Matrix{T}(undef, 1, 1)
    Hβσ2 = zeros(T, p)
    # Hβρ = Vector{T}(undef, p)
    res = Vector{T}(undef, n)
    t   = [tr(V)/2]
    q   = Vector{T}(undef, 1)
    xtx = transpose(X) * X
    storage_n = Vector{T}(undef, n)
    storage_p1 = Vector{T}(undef, p)
    storage_np = Matrix{T}(undef, n, p)
    storage_pp = Matrix{T}(undef, p, p)
    added_term_numerator = Matrix{T}(undef, n, p)
    added_term2 = Matrix{T}(undef, p, p)
    # constructor
    GaussianCopulaCSObs{T}(n, p, y, X, V, vec, ∇CSV, ∇β, ∇τ, ∇ρ, ∇σ2, Hβ, Hτ, Hρ, Hσ2, Hρσ2, Hβσ2,# Hβρ,
   res, t, q, xtx, storage_n, storage_p1, storage_np, storage_pp, added_term_numerator, added_term2)
end

"""
    GaussianCopulaCSModel
GaussianCopulaCSModel(gcs)
Gaussian copula compound symmetric CS model, which contains a vector of
`GaussianCopulaCSObs` as data, model parameters, and working arrays.
"""
struct GaussianCopulaCSModel{T <: BlasReal} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{GaussianCopulaCSObs{T}}
    Ytotal::T
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    ρ::Vector{T}            # autocorrelation parameter
    σ2::Vector{T}           # autoregressive noise parameter
    θ::Vector{T}
    # working arrays
    ∇β::Vector{T}   # gradient of beta from all observations
    ∇τ::Vector{T}   # gradient of tau from all observations
    ∇ρ::Vector{T}           # gradient of rho from all observations
    ∇σ2::Vector{T}          # gradient of sigmasquared from all observations
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    Hβ::Matrix{T}    # Hessian β from all observations
    Hτ::Matrix{T}    # Hessian τ from all observations
    Hρ::Matrix{T}    # Hessian from all observations
    Hσ2::Matrix{T}    # Hessian from all observations
    Hρσ2::Matrix{T}
    Hβσ2::Vector{T}
    Ainv::Matrix{T}
    Aevec::Matrix{T}
    M::Matrix{T}
    vcov::Matrix{T}
    ψ::Vector{T}
    # Hβρ::Vector{T}
    TR::Matrix{T}
    QF::Matrix{T}         # n-by-1 matrix with qik = res_i' Vi res_i
    storage_n::Vector{T}
    storage_m::Vector{T}
    storage_θ::Vector{T}
    penalized::Bool
end

function GaussianCopulaCSModel(gcs::Vector{GaussianCopulaCSObs{T}}; penalized::Bool = false) where T <: BlasReal
    n, p = length(gcs), size(gcs[1].X, 2)
    β   = Vector{T}(undef, p)
    τ   = [1.0]
    ρ = [1.0]
    σ2 = [1.0]
    θ   = Vector{T}(undef, 1) # initial MM update for crude estimate of noise
    ∇β  = Vector{T}(undef, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇ρ  = Vector{T}(undef, 1)
    ∇σ2  = Vector{T}(undef, 1)
    XtX = zeros(T, p, p) # sum_i xi'xi
    Hβ  = Matrix{T}(undef, p, p)
    Hτ  = Matrix{T}(undef, 1, 1)
    Hρ  = Matrix{T}(undef, 1, 1)
    Hσ2  = Matrix{T}(undef, 1, 1)
    Hρσ2 = Matrix{T}(undef, 1, 1)
    Hβσ2 = zeros(T, p)
    Ainv    = zeros(T, p + 3, p + 3)
    Aevec   = zeros(T, p + 3, p + 3)
    M       = zeros(T, p + 3, p + 3)
    vcov    = zeros(T, p + 3, p + 3)
    ψ       = Vector{T}(undef, p + 3)
    # Hβρ = Vector{T}(undef, p)
    TR  = Matrix{T}(undef, n, 1) # collect trace terms
    QF  = Matrix{T}(undef, n, 1)
    Ytotal = 0.0
    ntotal = 0.0
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        Ytotal  += sum(gcs[i].y)
        BLAS.axpy!(one(T), gcs[i].xtx, XtX)
        TR[i, :] = gcs[i].t
    end
    storage_n = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, 1)
    storage_θ = Vector{T}(undef, 1)
    GaussianCopulaCSModel{T}(gcs, Ytotal, ntotal, p, β, τ, ρ, σ2, θ,
        ∇β, ∇τ, ∇ρ, ∇σ2, XtX, Hβ, Hτ, Hρ, Hσ2, Hρσ2, Hβσ2, Ainv, Aevec,  M, vcov, ψ,
        TR, QF, storage_n, storage_m, storage_θ, penalized)
end

"""
    initialize_model!(gcm)
Initialize the linear regression parameters `β` and `τ=σ0^{-2}` by the least
squares solution.
"""
function initialize_model!(
    gcm::GaussianCopulaCSModel{T}
    ) where T <: BlasReal
    # accumulate sufficient statistics X'y
    xty = zeros(T, gcm.p)
    @inbounds for i in eachindex(gcm.data)
        BLAS.gemv!('T', one(T), gcm.data[i].X, gcm.data[i].y, one(T), xty)
    end
    # least square solution for β
    ldiv!(gcm.β, cholesky(Symmetric(gcm.XtX)), xty)
    @show gcm.β
    # accumulate residual sum of squares
    rss = zero(T)
    for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rss += abs2(norm(gcm.data[i].res))
    end
    println("initializing dispersion using residual sum of squares")
    gcm.τ[1] = gcm.ntotal / rss
    @show gcm.τ
    println("initializing CS noise paramter using MM-algorithm")
    fill!(gcm.ρ, 1.0)
    fill!(gcm.θ, 1.0)
    update_θ!(gcm)
    copyto!(gcm.σ2, gcm.θ)
    nothing
end

"""
    update_res!(gc, β)
Update the residual vector according to `β`.
"""
function update_res!(
    gc::GaussianCopulaCSObs{T},
    β::Vector{T}
    ) where T <: BlasReal
    copyto!(gc.res, gc.y)
    BLAS.gemv!('N', -one(T), gc.X, β, one(T), gc.res)
    gc.res
end

function update_res!(
    gcm::GaussianCopulaCSModel{T}
    ) where T <: BlasReal
    @inbounds for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
    end
    nothing
end

function standardize_res!(
    gc::GaussianCopulaCSObs{T},
    σinv::T
    ) where T <: BlasReal
    gc.res .*= σinv
end

function standardize_res!(
    gcm::GaussianCopulaCSModel{T}
    ) where T <: BlasReal
    σinv = sqrt(gcm.τ[1])
    # standardize residual
    @inbounds for i in eachindex(gcm.data)
        standardize_res!(gcm.data[i], σinv)
    end
    nothing
end

"""
    update_quadform!(gc)
Update the quadratic forms `(r^T V[k] r) / 2` according to the current residual `r`.
"""
function update_quadform!(gc::GaussianCopulaCSObs)
    gc.q[1] = dot(gc.res, mul!(gc.storage_n, gc.V, gc.res)) / 2
    gc.q
end

function update_θ_jensen!(
    gcm::GaussianCopulaCSModel{T},
    maxiter::Integer=50000,
    reltol::Number=1e-6,
    verbose::Bool=false) where T <: BlasReal
    rsstotal = zero(T)
    @inbounds for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rsstotal += abs2(norm(gcm.data[i].res))
        update_quadform!(gcm.data[i])
        gcm.QF[i, :] = gcm.data[i].q
    end
    # MM iteration
    for iter in 1:maxiter
        # store previous iterate
        copyto!(gcm.storage_θ, gcm.θ)
        # update τ
        mul!(gcm.storage_n, gcm.QF, gcm.θ) # gcm.storage_n[i] = q[i]
        gcm.τ[1] = update_τ(gcm.τ[1], gcm.storage_n, gcm.ntotal, rsstotal, 1)
        # numerator in the multiplicative update
        gcm.storage_n .= inv.(inv(gcm.τ[1]) .+ gcm.storage_n) # use newest τ to update θ
        mul!(gcm.storage_m, transpose(gcm.QF), gcm.storage_n)
        gcm.θ .*= gcm.storage_m
        # denominator in the multiplicative update
        mul!(gcm.storage_n, gcm.TR, gcm.storage_θ)
        gcm.storage_n .= inv.(1 .+ gcm.storage_n)
        mul!(gcm.storage_m, transpose(gcm.TR), gcm.storage_n)
        gcm.θ ./= gcm.storage_m
        # monotonicity diagnosis
        verbose && println(sum(log, 1 .+ gcm.τ[1] .* (gcm.QF * gcm.θ)) -
            sum(log, 1 .+ gcm.TR * gcm.θ) +
            gcm.ntotal / 2 * (log(gcm.τ[1]) - log(2π)) -
            rsstotal / 2 * gcm.τ[1])
        # convergence check
        gcm.storage_m .= gcm.θ .- gcm.storage_θ
        # norm(gcm.storage_m) < reltol * (norm(gcm.storage_θ) + 1) && break
        if norm(gcm.storage_m) < reltol * (norm(gcm.storage_θ) + 1)
            verbose && println("iters=$iter")
            break
        end
        verbose && iter == maxiter && @warn "maximum iterations $maxiter reached"
    end
    gcm.θ
end

# loglik_obs(::Normal, y, μ, wt, ϕ) = wt*GLM.logpdf(Normal(μ, sqrt(abs(ϕ))), y)

# this gets the loglikelihood from the glm.jl package for the component density
"""
    component_loglikelihood(gc::GaussianCopulaCSObs{T, D, Link}, τ, logl)
Calculates the loglikelihood of observing `y` given mean `μ`, a distribution
`d` using the GLM.jl package.
"""
function component_loglikelihood(gc::GaussianCopulaCSObs{T}, β::Vector{T}, τ::T) where {T <: BlasReal}
  logl = zero(T)
  μ = zeros(Integer(gc.n))
  mul!(μ, gc.X, β)
    @inbounds for j in eachindex(gc.y)
      logl += QuasiCopula.loglik_obs(Normal(), gc.y[j], μ[j], 1.0, inv(τ))
  end
  logl
end

function loglikelihood!(
    gc::GaussianCopulaCSObs{T},
    β::Vector{T},
    τ::T,
    ρ::T,
    σ2::T,
    needgrad::Bool = false,
    needhess::Bool = false;
    penalized::Bool = false
    ) where {T <: BlasReal}
    # n, p = size(gc.X, 1), size(gc.X, 2)
    needgrad = needgrad || needhess
    if needgrad
        fill!(gc.∇β, 0.0)
        fill!(gc.∇τ, 0.0)
        fill!(gc.∇σ2, 0.0)
        fill!(gc.∇ρ, 0.0)
        get_∇V!(gc)
    end
    needhess && fill!(gc.Hβ, 0.0)
    sqrtτ = sqrt(abs(τ))
    update_res!(gc, β)
    standardize_res!(gc, sqrtτ)
    rss  = abs2(norm(gc.res)) # RSS of standardized residual
    get_V!(ρ, gc)

    #evaluate copula loglikelihood
    # mul!(gc.storage_n, gc.V, gc.res) # storage_n = V[k] * res
    mul!(gc.storage_n, gc.V, gc.res, one(T), zero(T)) # storage_n = V[k] * res

    if needgrad
        BLAS.gemv!('T', σ2, gc.X, gc.storage_n, 1.0, gc.∇β) # stores ∇resβ*Γ*res (standardized residual)
    end
    q = dot(gc.res, gc.storage_n)
    tsum = 0.5 * gc.n * σ2
    qsum = 0.5 * σ2 * q
    # @show q
    c1 = 1 + tsum # 1 + tsum
    c2 = 1 + qsum # 1 + qsum
    # loglikelihood
    # logl = QuasiCopula.component_loglikelihood(gc, β, τ) #### make this function for normal
    logl = -log(c1) - (gc.n * log(2π) -  gc.n * log(abs(τ)) + rss) / 2
    # @show logl
    logl += log(c2)
    # add L2 ridge penalty
    if penalized
        logl -= 0.5 * (σ2)^2
    end
    # @show logl
    if needgrad
        inv1pq = inv(c2)
        # gradient with respect to rho
        # mul!(gc.storage_n, gc.∇CSV, gc.res) # storage_n = ∇CSV * res
        mul!(gc.storage_n, gc.∇CSV, gc.res, one(T), zero(T)) # storage_n = ∇ARV * res
        q2 = dot(gc.res, gc.storage_n) #
        gc.∇ρ .= inv(c2) * 0.5 * σ2 * q2
        # gradient with respect to sigma2
        gc.∇σ2 .= -(0.5 * gc.n * inv(c1)) .+ (inv(c2) * 0.5 * q)
        if penalized
            gc.∇σ2 .-= σ2
        end
      if needhess
            gc.Hρ .= - (inv(c2) * 0.5 * σ2 * q2)^2

            # hessian for sigma2
            gc.Hσ2 .= 0.25 * gc.n^2 * inv(c1)^2 - inv(c2)^2 * (0.25 * q^2)

            # hessian cross term for rho and sigma2
            gc.Hρσ2 .= 0.5 * q2 * inv1pq - 0.25 * σ2 * q * q2 * inv1pq^2

            # hessian cross term for beta and sigma2
            gc.Hβσ2 .= inv1pq * gc.∇β - 0.5 * q * inv1pq^2 * σ2 * gc.∇β

            BLAS.syrk!('L', 'N', -abs2(inv1pq), gc.∇β, 0.0, gc.Hβ) # only lower triangular
            gc.Hτ[1, 1] = - abs2(q * inv1pq / τ)
      end
      BLAS.gemv!('T', one(T), gc.X, gc.res, -inv1pq, gc.∇β)
      gc.∇β .*= sqrtτ
      gc.∇τ  .= (gc.n - rss + 2qsum * inv1pq) / 2τ
    end
    logl
end

function loglikelihood!(
    gcm::GaussianCopulaCSModel{T},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where {T <: BlasReal}
    if needgrad
        fill!(gcm.∇β, 0.0)
        fill!(gcm.∇τ, 0.0)
        fill!(gcm.∇ρ, 0.0)
        fill!(gcm.∇σ2, 0.0)
    end
    if needhess
        gcm.Hβ .= - gcm.XtX
        gcm.Hτ .= - gcm.ntotal / 2abs2(gcm.τ[1])
        fill!(gcm.Hρ, 0.0)
        fill!(gcm.Hσ2, 0.0)
        fill!(gcm.Hρσ2, 0.0)
        fill!(gcm.Hβσ2, 0.0)
        # @show gcm.Hβσ2
        # fill!(gcm.Hβρ, 0)
    end
    logl = zeros(Threads.nthreads())
    Threads.@threads for i in eachindex(gcm.data)
        @inbounds logl[Threads.threadid()] += loglikelihood!(gcm.data[i], gcm.β,
         gcm.τ[1], gcm.ρ[1], gcm.σ2[1], needgrad, needhess; penalized = gcm.penalized)
     end
     @inbounds for i in eachindex(gcm.data)
         if needgrad
             gcm.∇β .+= gcm.data[i].∇β
             gcm.∇τ .+= gcm.data[i].∇τ
             gcm.∇ρ .+= gcm.data[i].∇ρ
             gcm.∇σ2 .+= gcm.data[i].∇σ2
         end
         if needhess
             gcm.Hβ .+= gcm.data[i].Hβ
             gcm.Hτ .+= gcm.data[i].Hτ
             gcm.Hρ .+= gcm.data[i].Hρ
             gcm.Hσ2 .+= gcm.data[i].Hσ2
             gcm.Hρσ2 .+= gcm.data[i].Hρσ2
             gcm.Hβσ2 .+= gcm.data[i].Hβσ2
         end
     end
    needhess && (gcm.Hβ .*= gcm.τ[1])
    return sum(logl)
end


text:
using DrWatson
@quickactivate
using GenerativeAD
using BSON, FileIO

factor = 4

# downscale from 1024x1024 to 256x256
function downscale(x::AbstractArray{T,4}, factor) where T
	h,w,c,n = size(x)
	(w == h == 1024) ? nothing : error("only implemented for downscaling images from 1024 pixels")
	nw = nh = Int(h/factor)
	y = similar(x, nh, nw, c, n)
	step = factor-1
	for ni in 1:n
		for ic in 1:c
			for iw in 1:nw
				for ih in 1:nh
					y[ih, iw, ic, ni] = sum(x[(factor*ih-step):(factor*ih), (factor*iw-step):(factor*iw), ic, ni])/factor^2
				end
			end
		end
	end
	y
end

outpath = datadir("mvtec_ad/downscaled_data")
mkpath(outpath)

# start with wood
function downscale_data(category, factor)
	normal, anomalous = GenerativeAD.Datasets.load_mvtec_ad_data(category=category)
	dnormal, danomalous = map(x->downscale(x,factor), (normal, anomalous))
	s = Int(1024/factor)
	fname = joinpath(outpath, "$(category)_$(s).bson")
	save(fname, Dict(:normal => dnormal, :anomalous => danomalous))
	@info "saved to $fname"
end

map(f->map(c->downscale_data(c,f), ["wood", "transistor", "grid"]), [4, 8, 16, 32])