text:
# BCF Record
# ==========

type BCFRecord
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # fields
    sharedlen::UInt32
    indivlen::UInt32
end

function BCFRecord()
    return BCFRecord(UInt8[], 1:0, 0, 0)
end

function isfilled(record::BCFRecord)
    return !isempty(record.filled)
end

function datarange(record::BCFRecord)
    return record.filled
end

function checkfilled(record::BCFRecord)
    if !isfilled(record)
        throw(ArgumentError("unfilled BCF record"))
    end
end

function Base.:(==)(record1::BCFRecord, record2::BCFRecord)
    if isfilled(record1) == isfilled(record2) == true
        r1 = datarange(record1)
        r2 = datarange(record2)
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function Base.copy(record::BCFRecord)
    return BCFRecord(record.data[record.filled], record.filled, record.sharedlen, record.indivlen)
end

function BCFRecord(base::BCFRecord;
                   chromosome=nothing, position=nothing, quality=nothing,
                   identifier=nothing, reference=nothing, alternate=nothing,
                   filter=nothing, information=nothing, genotype=nothing)
    checkfilled(base)
    data = base.data[1:24]

    if chromosome != nothing
        store!(data, 0, convert(Int32, chromosome - 1))
    end
    if position != nothing
        store!(data, 4, convert(Int32, position - 1))
    end
    if quality != nothing
        store!(data, 12, convert(Float32, quality))
    end

    offset = boffset = 24
    if identifier == nothing
        offset, boffset = copyvec!(data, offset, base.data, boffset)
    else
        offset = storestr!(data, offset, string(identifier))
        boffset = skipvec(base.data, boffset)
    end

    if reference == nothing
        offset, boffset = copyvec!(data, offset, base.data, boffset)
    else
        reference = string(reference)
        offset = storestr!(data, offset, reference)
        if sizeof(reference) != rlen(base)
            store!(data, 8, convert(Int32, sizeof(reference)))
        end
        boffset = skipvec(base.data, boffset)
    end

    n = n_allele(base)
    if alternate == nothing
        for _ in 1:n-1
            offset, boffset = copyvec!(data, offset, base.data, boffset)
        end
    else
        if !isa(alternate, Vector)
            alternate = [alternate]
        end
        for x in alternate
            offset = storestr!(data, offset, string(x))
        end
        for _ in 1:n-1
            boffset = skipvec(base.data, boffset)
        end
        if length(alternate) != n
            n_allele_info = load(UInt32, data, 16)[1]
            n_allele_info = (n_allele_info & 0x0000ffff) | (UInt32(length(alternate) + 1) << 16)
            store!(data, 16, n_allele_info)
        end
    end

    if filter == nothing
        offset, boffset = copyvec!(data, offset, base.data, boffset)
    else
        if !isa(filter, Vector)
            filter = [filter]
        end
        offset = storevec!(data, offset, convert(Vector{Int8}, filter - 1))
        boffset = skipvec(base.data, boffset)
    end

    n = n_info(base)
    if information == nothing
        for _ in 1:n
            # copy key and value(s)
            offset, boffset = copyvec!(data, offset, base.data, boffset)
            offset, boffset = copyvec!(data, offset, base.data, boffset)
        end
    else
        if !isa(information, Associative)
            throw(ArgumentError("info must be an associative object"))
        end
        keyvec = Vector{Int8}(1)
        for (key, val) in information
            if !isa(key, Integer)
                throw(ArgumentError("info key must be an integer"))
            end
            keyvec[1] = key
            offset = storevec!(data, offset, keyvec)
            if isa(val, String)
                offset = storestr!(data, offset, val)
            elseif !isa(val, Vector)
                offset = storevec!(data, offset, [val])
            else
                offset = storevec!(data, offset, val)
            end
        end
        if length(information) != n
            if length(information) > typemax(UInt16)
                throw(ArgumentError("too many info fields"))
            end
            n_allele_info = load(UInt32, data, 16)[1]
            n_allele_info = (n_allele_info & 0xffff0000) | UInt32(length(information))
            store!(data, 16, n_allele_info)
        end
        for _ in 1:n
            # skip key and value(s)
            boffset = skipvec(base.data, boffset)
            boffset = skipvec(base.data, boffset)
        end
    end
    sharedlen = offset

    n = n_format(base)
    if genotype == nothing
        N = n_sample(base)
        for _ in 1:n
            # copy key
            offset, boffset = copyvec!(data, offset, base.data, boffset)
            # copy value(s)
            head, boffset = loadvechead(base.data, boffset)
            offset = storevechead!(data, offset, head)
            for _ in 1:N
                offset, boffset = copyvecbody!(data, offset, base.data, boffset, head)
            end
        end
    else
        error("modifying genotype is yet supported")
    end

    return BCFRecord(data, 1:endof(data), sharedlen, offset - sharedlen)
end

function Base.show(io::IO, record::BCFRecord)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   chromosome: ", chromosome(record))
        println(io, "     position: ", leftposition(record))
        println(io, "   identifier: ", identifier(record))
        println(io, "    reference: ", reference(record))
        println(io, "    alternate: ", join(alternate(record), ' '))
        println(io, "      quality: ", quality(record))
        println(io, "       filter: ", join(filter_(record), ' '))
          print(io, "  information: ", information(record))
    else
        print(io, " <not filled>")
    end
end

function chromosome(rec::BCFRecord)
    checkfilled(rec)
    return load(Int32, rec.data, 0)[1] % Int + 1
end

function leftposition(rec::BCFRecord)
    checkfilled(rec)
    return load(Int32, rec.data, 4)[1] % Int + 1
end

function rlen(rec::BCFRecord)
    checkfilled(rec)
    return load(Int32, rec.data, 8)[1] % Int
end

function quality(rec::BCFRecord)
    checkfilled(rec)
    # 0x7F800001 is a missing value.
    return load(Float32, rec.data, 12)[1]
end

function n_allele(rec::BCFRecord)
    checkfilled(rec)
    return (load(Int32, rec.data, 16)[1] >> 16) % Int
end

function n_info(rec::BCFRecord)
    checkfilled(rec)
    return (load(Int32, rec.data, 16)[1] & 0x0000ffff) % Int
end

function n_format(rec::BCFRecord)
    checkfilled(rec)
    return (load(UInt32, rec.data, 20)[1] >> 24) % Int
end

function n_sample(rec::BCFRecord)
    checkfilled(rec)
    return (load(UInt32, rec.data, 20)[1] & 0x000000ff) % Int
end

function identifier(rec::BCFRecord)
    checkfilled(rec)
    offset = 24
    return loadstr(rec.data, offset)[1]
end

function reference(rec::BCFRecord)
    checkfilled(rec)
    # skip ID
    offset = 24
    len, offset = loadveclen(rec.data, offset)
    # load REF
    return loadstr(rec.data, offset + len)[1]
end

function alternate(rec::BCFRecord)
    checkfilled(rec)
    # skip ID and REF
    offset = 24
    len, offset = loadveclen(rec.data, offset)
    len, offset = loadveclen(rec.data, offset + len)
    # load ALTs
    N = n_allele(rec) - 1
    alt = Vector{String}(N)
    for n in 1:N
        str, offset = loadstr(rec.data, offset + len)
        alt[n] = str
        len = sizeof(str)
    end
    return alt
end

function filter_(rec::BCFRecord)
    checkfilled(rec)
    # skip ID, REF and ALTs
    offset = 24
    len = 0
    for _ in 1:n_allele(rec)+1
        len, offset = loadveclen(rec.data, offset + len)
    end
    # load FILTER
    return loadvec(rec.data, offset + len)[1] .+ 1
end

function information(rec::BCFRecord; simplify::Bool=true)
    checkfilled(rec)
    # skip ID, REF, ALTs and FILTER
    offset::Int = 24
    len = 0
    for _ in 1:n_allele(rec)+2
        len, offset = loadveclen(rec.data, offset + len)
    end
    offset += len
    # load INFO
    ret = Vector{Tuple{Int,Any}}(n_info(rec))
    for i in 1:endof(ret)
        key, offset = loadvec(rec.data, offset)
        @assert length(key) == 1
        val, offset = loadvec(rec.data, offset)
        if simplify
            if isa(val, Vector) && length(val) == 1
                val = val[1]
            elseif isa(val, Vector{Void}) && isempty(val)
                val = nothing
            end
        end
        ret[i] = (key[1], val)
    end
    return ret
end

function information(rec::BCFRecord, key::Integer; simplify::Bool=true)
    checkfilled(rec)
    # skip ID, REF, ALTs and FILTER
    offset::Int = 24
    len = 0
    for _ in 1:n_allele(rec)+2
        len, offset = loadveclen(rec.data, offset + len)
    end
    offset += len
    # load INFO
    for _ in 1:n_info(rec)
        k, offset = loadvec(rec.data, offset)
        @assert length(k) == 1
        if k[1] == key
            val, offset = loadvec(rec.data, offset)
            if simplify
                if isa(val, Vector) && length(val) == 1
                    val = val[1]
                elseif isa(val, Vector{Void}) && isempty(val)
                    val = nothing
                end
            end
            return val
        else
            offset = skipvec(rec.data, offset)
        end
    end
    throw(KeyError(key))
end

function genotype(rec::BCFRecord)
    checkfilled(rec)
    offset::Int = rec.sharedlen
    N = n_sample(rec)
    ret = Tuple{Int,Any}[]
    for j in 1:n_format(rec)
        key, offset = loadvec(rec.data, offset)
        @assert length(key) == 1
        head, offset = loadvechead(rec.data, offset)
        vals = Vector{Any}[]
        for n in 1:N
            val, offset = loadvecbody(rec.data, offset, head)
            push!(vals, val)
        end
        push!(ret, (key[1], vals))
    end
    return ret
end

function genotype(rec::BCFRecord, index::Integer)
    return [(k, geno[index]) for (k, geno) in genotype(rec)]
end

function genotype(rec::BCFRecord, index::Integer, key::Integer)
    checkfilled(rec)
    N = n_sample(rec)
    offset::Int = rec.sharedlen
    for j in 1:n_format(rec)
        k, offset = loadvec(rec.data, offset)
        @assert length(k) == 1
        head, offset = loadvechead(rec.data, offset)
        for n in 1:N
            if k[1] == key && n == index
                return loadvecbody(rec.data, offset, head)[1]
            else
                offset = skipvecbody(rec.data, offset, head)
            end
        end
    end
    throw(KeyError(key))
end

function genotype{T<:Integer}(rec::BCFRecord, indexes::AbstractVector{T}, key::Integer)
    checkfilled(rec)
    N = n_sample(rec)
    offset::Int = rec.sharedlen
    for j in 1:n_format(rec)
        k, offset = loadvec(rec.data, offset)
        @assert length(k) == 1
        head, offset = loadvechead(rec.data, offset)
        if k[1] == key
            vals = Vector{Any}[]
            for n in 1:N
                val, offset = loadvecbody(rec.data, offset, head)
                push!(vals, val)
            end
            return vals[indexes]
        else
            for n in 1:N
                offset = skipvecbody(rec.data, offset, head)
            end
        end
    end
    throw(KeyError(key))
end

function genotype(rec::BCFRecord, ::Colon, key::Integer)
    return genotype(rec, 1:n_sample(rec), key)
end

function gt(x::Int8)
    allele = (x >> 1) - 1
    phased = x & 1 != 0
    return allele, phased
end

# Resolve type encoding.
function bcftype(b::UInt8)
    b &= 0x0f
    if b == 0x00
        # this is not defined in the BCF specs v2.2 but used in htslib
        return :void
    elseif b == 0x01
        return :int8
    elseif b == 0x02
        return :int16
    elseif b == 0x03
        return :int32
    elseif b == 0x05
        return :float32
    elseif b == 0x07
        return :character
    else
        return :unused
    end
end

function bcftypetag(typ::Symbol, len::Int)
    b = 0x00
    if typ == :void
        b |= 0x00
        @assert len == 0
    elseif typ == :int8
        b |= 0x01
    elseif typ == :int16
        b |= 0x02
    elseif typ == :int32
        b |= 0x03
    elseif typ == :float32
        b |= 0x05
    elseif typ == :character
        b |= 0x07
    else
        throw(ArgumentError("invalid type name: $(typ)"))
    end

    if len < 0
        throw(ArgumentError("length must be non-negative"))
    elseif len < 15
        b |= (len % UInt8) << 4
    else
        b |= (15 % UInt8)  << 4
    end
    return b
end

function bcftypesize(typ::Symbol)
    return typ == :int8      ? 1 :
           typ == :int16     ? 2 :
           typ == :int32     ? 4 :
           typ == :float32   ? 4 :
           typ == :character ? 1 :
           typ == :void      ? 0 :
           error("size unknown")
end

function load{T}(::Type{T}, data::Vector{UInt8}, offset::Int)
    if offset + sizeof(T) > length(data)
        throw(BoundsError())
    end
    return unsafe_load(convert(Ptr{T}, pointer(data, offset + 1))), offset + sizeof(T)
end

# Load a string; this is basically same as `loadvec` but type-stable.
function loadstr(data::Vector{UInt8}, offset::Int)
    (typ, len), offset = loadvechead(data, offset)
    @assert typ == :character
    return String(data[offset+1:offset+len]), offset
end

function loadvec(data::Vector{UInt8}, offset::Int)
    head, offset = loadvechead(data, offset)
    return loadvecbody(data, offset, head)
end

function loadveclen(data::Vector{UInt8}, offset::Int)
    (_, len), offset = loadvechead(data, offset)
    return len, offset
end

function loadvechead(data::Vector{UInt8}, offset::Int)
    b = data[offset+1]
    typ = bcftype(b)
    len::Int = b >> 4
    if len < 15
        return (typ, len), offset + 1
    end
    lenvec, offset = loadvec(data, offset + 1)
    @assert length(lenvec) == 1
    return (typ, Int(lenvec[1])), offset
end

function loadvecbody(data::Vector{UInt8}, offset::Int, head::Tuple{Symbol,Int})
    t, len = head
    T = t == :int8 ? Int8 :
        t == :int16 ? Int16 :
        t == :int32 ? Int32 :
        t == :float32 ? Float32 :
        t == :character ? UInt8 : Void
    ret = Vector{T}(len)
    for i in 1:len
        ret[i], offset = load(T, data, offset)
    end
    if t == :character
        return String(ret), offset
    else
        return ret, offset
    end
end

function skipvec(data::Vector{UInt8}, offset::Int)
    (t, len), offset = loadvechead(data, offset)
    return offset + bcftypesize(t) * len
end

function skipvecbody(data::Vector{UInt8}, offset::Int, head::Tuple{Symbol,Int})
    t, len = head
    return offset + bcftypesize(t) * len
end

function store!{T}(data::Vector{UInt8}, offset::Int, val::T)
    if offset + sizeof(T) > length(data)
        throw(BoundsError())
    end
    unsafe_store!(convert(Ptr{T}, pointer(data, offset + 1)), val)
    return offset + sizeof(T)
end

function storestr!(data::Vector{UInt8}, offset::Int, str::String)
    if !isascii(str)
        throw(ArgumentError("string must be ASCII"))
    end
    len = sizeof(str)
    offset = storevechead!(data, offset, (:character, len))
    if length(data) < offset + len
        resize!(data, offset + len)
    end
    memcpy(pointer(data, offset + 1), pointer(str), len)
    return offset + len
end

function storevec!{T}(data::Vector{UInt8}, offset::Int, vec::Vector{T})
    t = T == Int8 ? :int8 :
        T == Int16 ? :int16 :
        T == Int32 ? :int32 :
        T == Float32 ? :float32 :
        error("unsupported type: $(T)")
    len = length(vec)
    offset = storevechead!(data, offset, (t, len))
    return storevecbody!(data, offset, vec)
end

function storevechead!(data::Vector{UInt8}, offset::Int, head::Tuple{Symbol,Int})
    t, len = head
    if length(data) < offset + 1
        resize!(data, offset + 1)
    end
    data[offset+1] = bcftypetag(t, len)
    if len < 15
        return offset + 1
    elseif len < typemax(Int8)
        return storevec!(data, offset + 1, Int8[len])
    elseif len < typemax(Int16)
        return storevec!(data, offset + 1, Int16[len])
    elseif len < typemax(Int32)
        return storevec!(data, offset + 1, Int32[len])
    else
        error("too long vector")
    end
end

function storevecbody!{T}(data::Vector{UInt8}, offset::Int, vec::Vector{T})
    len = length(vec)
    n = len * sizeof(T)
    if length(data) < offset + n
        resize!(data, offset + n)
    end
    memcpy(pointer(data, offset + 1), pointer(vec), n)
    return offset + n
end

function copyvec!(dst::Vector{UInt8}, doffset::Int, src::Vector{UInt8}, soffset::Int)
    head, soffset = loadvechead(src, soffset)
    doffset = storevechead!(dst, doffset, head)
    return copyvecbody!(dst, doffset, src, soffset, head)
end

function copyvecbody!(dst::Vector{UInt8}, doffset::Int, src::Vector{UInt8}, soffset::Int, head::Tuple{Symbol,Int})
    t, len = head
    n = bcftypesize(t) * len
    if length(dst) < doffset + n
        resize!(dst, doffset + n)
    end
    copy!(dst, doffset + 1, src, soffset + 1, n)
    return doffset + n, soffset + n
end

function memcpy(dst::Ptr, src::Ptr, n::Int)
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Csize_t), dst, src, n)
    return nothing
end


text:
#=

=#

using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
reset_default_graph()

include("nnutil.jl")
include("NNTrussPull_Domain.jl")
testtype = "NeuralNetwork1D"
prop = Dict("name"=> testtype, "rho"=> 8000.0, "E"=> 200e3, "nu"=> 0.45,
           "sigmaY"=>0.3e3, "K"=>1/9*200e3, "B"=> 0.0, "A0"=> 0.005, "nn"=>post_nn)

n_data = [1,2,4,5]
nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(n_data[1])
elements = []
for i = 1:nx
    elnodes = [i, i+1]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))
end
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)

losses = Array{PyObject}(undef, length(n_data))
for (i, tid) in enumerate(n_data)
    state_history, fext_history = read_data("$(@__DIR__)/Data/$tid.dat")
    losses[i] = DynamicMatLawLoss(domain, globdat, state_history, fext_history, nn,Δt)
end
loss = sum(losses)

sess = Session(); init(sess)
@show run(sess, loss)

if !isdir("Data/$(nntype)")
    mkdir("Data/$(nntype)")
end

start_id = 3
ADCME.load(sess, "Data/$(nntype)/NNPreLSfit_nn$(idx)_ite$(start_id).mat")
for i = 1:10
    println("************************** Outer Iteration = $i ************************** ")
    BFGS!(sess, loss, 1000)
    ADCME.save(sess, "Data/$(nntype)/trained_nn$(idx)_from$(start_id)_ite$(i).mat")
end

# X, Y = prepare_strain_stress_data1D(domain)
# x = (constant(X[:,1]), constant(X[:,2]), constant(X[:,3]))
# y = squeeze(nn(x...))
# close("all")
# out = run(sess, y)
# plot(X[:,1], out,"+", label="NN")
# plot(X[:,1], Y, ".", label="Exact")
# legend()

ext: jl


text:
using TravelingSalesmanExact, GLPK
using TravelingSalesmanHeuristics
using TravelingSalesmanExact: ATT, euclidean_distance
using TravelingSalesmanBenchmarks
using Plots
using Random
gr(fmt=:svg)

Random.seed!(758)
N = 50
cost = 100*rand(N, N)
cost = (cost + transpose(cost) )/2

t_exact, c_exact = get_optimal_tour(cost, with_optimizer(GLPK.Optimizer))

c(q) = solve_tsp(cost; quality_factor = q)[2]

qs = range(10, stop = 100, step = 10)

plot(qs, c, xlabel="quality", ylabel="Cost", label="solve_tsp 1", title="random symmetric cost matrix")
for j = 2:5
    plot!(qs, c, label="solve_tsp $j")
end
hline!([c_exact], label="Exact cost")

file = @isdefined(WEAVE_ARGS) ? WEAVE_ARGS[:file] : nothing
TravelingSalesmanBenchmarks.bench_footer(file)

max_stars_repo_name: ericphanson/TravelingSalesmanBenchmarks


text:
@doc raw"""
    AbstractGroupOperation

Abstract type for smooth binary operations $∘$ on elements of a Lie group $\mathcal{G}$:
```math
∘ : \mathcal{G} × \mathcal{G} → \mathcal{G}
```
An operation can be either defined for a specific group manifold over
number system `𝔽` or in general, by defining for an operation `Op` the following methods:

    identity_element!(::AbstractDecoratorManifold, q, q)
    inv!(::AbstractDecoratorManifold, q, p)
    _compose!(::AbstractDecoratorManifold, x, p, q)

Note that a manifold is connected with an operation by wrapping it with a decorator,
[`AbstractDecoratorManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/decorator.html#ManifoldsBase.AbstractDecoratorManifold) using the [`IsGroupManifold`](@ref) to specify the operation.
For a concrete case the concrete wrapper [`GroupManifold`](@ref) can be used.
"""
abstract type AbstractGroupOperation end

"""
    IsGroupManifold{O<:AbstractGroupOperation} <: AbstractTrait

A trait to declare an [`AbstractManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifold)  as a manifold with group structure
with operation of type `O`.

Using this trait you can turn a manifold that you implement _implictly_ into a Lie group.
If you wish to decorate an existing manifold with one (or different) [`AbstractGroupAction`](@ref)s,
see [`GroupManifold`](@ref).

# Constructor

    IsGroupManifold(op)
"""
struct IsGroupManifold{O<:AbstractGroupOperation} <: AbstractTrait
    op::O
end

"""
    AbstractInvarianceTrait <: AbstractTrait

A common supertype for anz [`AbstractTrait`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/decorator.html#ManifoldsBase.AbstractTrait) related to metric invariance
"""
abstract type AbstractInvarianceTrait <: AbstractTrait end

"""
    HasLeftInvariantMetric <: AbstractInvarianceTrait

Specify that a certain the metric of a [`GroupManifold`](@ref) is a left-invariant metric
"""
struct HasLeftInvariantMetric <: AbstractInvarianceTrait end

direction(::HasLeftInvariantMetric) = LeftAction()
direction(::Type{HasLeftInvariantMetric}) = LeftAction()

"""
    HasRightInvariantMetric <: AbstractInvarianceTrait

Specify that a certain the metric of a [`GroupManifold`](@ref) is a right-invariant metric
"""
struct HasRightInvariantMetric <: AbstractInvarianceTrait end

direction(::HasRightInvariantMetric) = RightAction()
direction(::Type{HasRightInvariantMetric}) = RightAction()

"""
    HasBiinvariantMetric <: AbstractInvarianceTrait

Specify that a certain the metric of a [`GroupManifold`](@ref) is a bi-invariant metric
"""
struct HasBiinvariantMetric <: AbstractInvarianceTrait end
function parent_trait(::HasBiinvariantMetric)
    return ManifoldsBase.TraitList(HasLeftInvariantMetric(), HasRightInvariantMetric())
end

"""
    is_group_manifold(G::GroupManifold)
    is_group_manifoldd(G::AbstractManifold, o::AbstractGroupOperation)

returns whether an [`AbstractDecoratorManifold`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/decorator.html#ManifoldsBase.AbstractDecoratorManifold) is a group manifold with
[`AbstractGroupOperation`](@ref) `o`.
For a [`GroupManifold`](@ref) `G` this checks whether the right operations is stored within `G`.
"""
is_group_manifold(::AbstractManifold, ::AbstractGroupOperation) = false

@trait_function is_group_manifold(M::AbstractDecoratorManifold, op::AbstractGroupOperation)
function is_group_manifold(
    ::TraitList{<:IsGroupManifold{<:O}},
    ::AbstractDecoratorManifold,
    ::O,
) where {O<:AbstractGroupOperation}
    return true
end
@trait_function is_group_manifold(M::AbstractDecoratorManifold)
is_group_manifold(::AbstractManifold) = false
function is_group_manifold(
    t::TraitList{<:IsGroupManifold{<:AbstractGroupOperation}},
    M::AbstractDecoratorManifold,
)
    return is_group_manifold(M, t.head.op)
end

"""
    ActionDirection

Direction of action on a manifold, either [`LeftAction`](@ref) or [`RightAction`](@ref).
"""
abstract type ActionDirection end

"""
    LeftAction()

Left action of a group on a manifold.
"""
struct LeftAction <: ActionDirection end

"""
    RightAction()

Right action of a group on a manifold.
"""
struct RightAction <: ActionDirection end

"""
    switch_direction(::ActionDirection)

Returns a [`RightAction`](@ref) when given a [`LeftAction`](@ref) and vice versa.
"""
switch_direction(::ActionDirection)
switch_direction(::LeftAction) = RightAction()
switch_direction(::RightAction) = LeftAction()

@doc raw"""
    Identity{O<:AbstractGroupOperation}

Represent the group identity element ``e ∈ \mathcal{G}`` on a Lie group ``\mathcal G``
with [`AbstractGroupOperation`](@ref) of type `O`.

Similar to the philosophy that points are agnostic of their group at hand, the identity
does not store the group `g` it belongs to. However it depends on the type of the [`AbstractGroupOperation`](@ref) used.

See also [`identity_element`](@ref) on how to obtain the corresponding [`AbstractManifoldPoint`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.AbstractManifoldPoint) or array representation.

# Constructors

    Identity(G::AbstractDecoratorManifold{𝔽})
    Identity(o::O)
    Identity(::Type{O})

create the identity of the corresponding subtype `O<:`[`AbstractGroupOperation`](@ref)
"""
struct Identity{O<:AbstractGroupOperation} end

@trait_function Identity(M::AbstractDecoratorManifold)
function Identity(
    ::TraitList{<:IsGroupManifold{O}},
    ::AbstractDecoratorManifold,
) where {O<:AbstractGroupOperation}
    return Identity{O}()
end
Identity(::O) where {O<:AbstractGroupOperation} = Identity(O)
Identity(::Type{O}) where {O<:AbstractGroupOperation} = Identity{O}()

# To ensure allocate_result_type works in general if idenitty apears in the tuple
number_eltype(::Identity) = Bool

@doc raw"""
    identity_element(G)

Return a point representation of the [`Identity`](@ref) on the [`IsGroupManifold`](@ref) `G`.
By default this representation is the default array or number representation.
It should return the corresponding default representation of ``e`` as a point on `G` if
points are not represented by arrays.
"""
identity_element(G::AbstractDecoratorManifold)
@trait_function identity_element(G::AbstractDecoratorManifold)
function identity_element(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold)
    q = allocate_result(G, identity_element)
    return identity_element!(G, q)
end

@trait_function identity_element!(G::AbstractDecoratorManifold, p)

function allocate_result(G::AbstractDecoratorManifold, ::typeof(identity_element))
    return zeros(representation_size(G)...)
end

@doc raw"""
    identity_element(G::AbstractDecoratorManifold, p)

Return a point representation of the [`Identity`](@ref) on the [`IsGroupManifold`](@ref) `G`,
where `p` indicates the type to represent the identity.
"""
identity_element(G::AbstractDecoratorManifold, p)
@trait_function identity_element(G::AbstractDecoratorManifold, p)
function identity_element(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p)
    q = allocate_result(G, identity_element, p)
    return identity_element!(G, q)
end
function check_size(
    ::TraitList{<:IsGroupManifold{<:O}},
    M::AbstractDecoratorManifold,
    ::Identity{<:O},
) where {O<:AbstractGroupOperation}
    return nothing
end
function check_size(::EmptyTrait, M::AbstractDecoratorManifold, e::Identity)
    return DomainError(0, "$M seems to not be a group manifold with $e.")
end

@doc raw"""
    is_identity(G::AbstractDecoratorManifold, q; kwargs)

Check whether `q` is the identity on the [`IsGroupManifold`](@ref) `G`, i.e. it is either
the [`Identity`](@ref)`{O}` with the corresponding [`AbstractGroupOperation`](@ref) `O`, or
(approximately) the correct point representation.
"""
is_identity(G::AbstractDecoratorManifold, q)
@trait_function is_identity(G::AbstractDecoratorManifold, q; kwargs...)
function is_identity(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    q;
    kwargs...,
)
    return isapprox(G, identity_element(G), q; kwargs...)
end
function is_identity(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    ::Identity{O};
    kwargs...,
) where {O<:AbstractGroupOperation}
    return true
end
function is_identity(
    ::TraitList{<:IsGroupManifold},
    ::AbstractDecoratorManifold,
    ::Identity;
    kwargs...,
)
    return false
end

@inline function isapprox(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p::Identity{O},
    q;
    kwargs...,
) where {O<:AbstractGroupOperation}
    return is_identity(G, q; kwargs...)
end
@inline function isapprox(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p,
    q::Identity{O};
    kwargs...,
) where {O<:AbstractGroupOperation}
    return is_identity(G, p; kwargs...)
end
function isapprox(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p::Identity{O},
    q::Identity{O};
    kwargs...,
) where {O<:AbstractGroupOperation}
    return true
end
function isapprox(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p::Identity{O},
    q::Identity;
    kwargs...,
) where {O<:AbstractGroupOperation}
    return false
end
function isapprox(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p::Identity,
    q::Identity{O};
    kwargs...,
) where {O<:AbstractGroupOperation}
    return false
end

@inline function isapprox(
    ::TraitList{IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p::Identity{O},
    X,
    Y;
    kwargs...,
) where {𝔽,O<:AbstractGroupOperation}
    return isapprox(G, identity_element(G), X, Y; kwargs...)
end
function Base.isapprox(
    ::TraitList{<:IsGroupManifold},
    ::AbstractDecoratorManifold,
    ::Identity,
    ::Identity;
    kwargs...,
)
    return false
end

function Base.show(io::IO, ::Identity{O}) where {O<:AbstractGroupOperation}
    return print(io, "Identity($O)")
end

function check_point(
    ::TraitList{<:IsGroupManifold{O1}},
    G::AbstractDecoratorManifold,
    e::Identity{O2};
    kwargs...,
) where {O1<:AbstractGroupOperation,O2<:AbstractGroupOperation}
    return DomainError(
        e,
        "The Identity $e does not lie on $G, since its the identity with respect to $O2 and not $O1.",
    )
end

function is_point(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    e::Identity,
    te=false;
    kwargs...,
)
    ie = is_identity(G, e; kwargs...)
    (te && !ie) && throw(DomainError(e, "The provided identity is not a point on $G."))
    return ie
end

function is_vector(
    t::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    e::Identity,
    X,
    te=false,
    cbp=true;
    kwargs...,
)
    if cbp
        ie = is_point(G, e; kwargs...)
        (!te) && return ie
    end
    return is_vector(next_trait(t), G, identity_element(G), X, te, false; kwargs...)
end

@doc raw"""
    adjoint_action(G::AbstractDecoratorManifold, p, X)

Adjoint action of the element `p` of the Lie group `G` on the element `X`
of the corresponding Lie algebra.

It is defined as the differential of the group authomorphism ``Ψ_p(q) = pqp⁻¹`` at
the identity of `G`.

The formula reads
````math
\operatorname{Ad}_p(X) = dΨ_p(e)[X]
````
where $e$ is the identity element of `G`.

Note that the adjoint representation of a Lie group isn't generally faithful.
Notably the adjoint representation of SO(2) is trivial.
"""
adjoint_action(G::AbstractDecoratorManifold, p, X)
@trait_function adjoint_action(G::AbstractDecoratorManifold, p, Xₑ)
function adjoint_action(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p, Xₑ)
    Xₚ = translate_diff(G, p, Identity(G), Xₑ, LeftAction())
    Y = inverse_translate_diff(G, p, p, Xₚ, RightAction())
    return Y
end

@trait_function adjoint_action!(G::AbstractDecoratorManifold, Y, p, Xₑ)
function adjoint_action!(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    Y,
    p,
    Xₑ,
)
    Xₚ = translate_diff(G, p, Identity(G), Xₑ, LeftAction())
    inverse_translate_diff!(G, Y, p, p, Xₚ, RightAction())
    return Y
end

@doc raw"""
    inv(G::AbstractDecoratorManifold, p)

Inverse $p^{-1} ∈ \mathcal{G}$ of an element $p ∈ \mathcal{G}$, such that
$p \circ p^{-1} = p^{-1} \circ p = e ∈ \mathcal{G}$, where $e$ is the [`Identity`](@ref)
element of $\mathcal{G}$.
"""
inv(::AbstractDecoratorManifold, ::Any...)
@trait_function Base.inv(G::AbstractDecoratorManifold, p)
function Base.inv(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p)
    q = allocate_result(G, inv, p)
    return inv!(G, q, p)
end

function Base.inv(
    ::TraitList{IsGroupManifold{O}},
    ::AbstractDecoratorManifold,
    e::Identity{O},
) where {O<:AbstractGroupOperation}
    return e
end

@trait_function inv!(G::AbstractDecoratorManifold, q, p)
function inv!(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, q, p)
    return inv!(G.manifold, q, p)
end

function inv!(
    ::TraitList{IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    q,
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return identity_element!(G, q)
end

function Base.copyto!(
    ::TraitList{IsGroupManifold{O}},
    ::AbstractDecoratorManifold,
    e::Identity{O},
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return e
end
function Base.copyto!(
    ::TraitList{IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    p,
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return identity_element!(G, p)
end

@doc raw"""
    compose(G::AbstractDecoratorManifold, p, q)

Compose elements ``p,q ∈ \mathcal{G}`` using the group operation ``p \circ q``.

For implementing composition on a new group manifold, please overload `_compose`
instead so that methods with [`Identity`](@ref) arguments are not ambiguous.
"""
compose(::AbstractDecoratorManifold, ::Any...)

@trait_function compose(G::AbstractDecoratorManifold, p, q)
function compose(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, p, q)
    return _compose(G, p, q)
end
function compose(
    ::AbstractDecoratorManifold,
    ::Identity{O},
    p,
) where {O<:AbstractGroupOperation}
    return p
end
function compose(
    ::AbstractDecoratorManifold,
    p,
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return p
end
function compose(
    ::AbstractDecoratorManifold,
    e::Identity{O},
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return e
end

function _compose(G::AbstractDecoratorManifold, p, q)
    x = allocate_result(G, compose, p, q)
    return _compose!(G, x, p, q)
end

@trait_function compose!(M::AbstractDecoratorManifold, x, p, q)

function compose!(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, x, q, p)
    return _compose!(G, x, q, p)
end
function compose!(
    G::AbstractDecoratorManifold,
    q,
    p,
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return copyto!(G, q, p)
end
function compose!(
    G::AbstractDecoratorManifold,
    q,
    ::Identity{O},
    p,
) where {O<:AbstractGroupOperation}
    return copyto!(G, q, p)
end
function compose!(
    G::AbstractDecoratorManifold,
    q,
    ::Identity{O},
    e::Identity{O},
) where {O<:AbstractGroupOperation}
    return identity_element!(G, q)
end
function compose!(
    ::AbstractDecoratorManifold,
    e::Identity{O},
    ::Identity{O},
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return e
end

transpose(e::Identity) = e

@doc raw"""
    hat(M::AbstractDecoratorManifold{𝔽,O}, ::Identity{O}, Xⁱ) where {𝔽,O<:AbstractGroupOperation}

Given a basis $e_i$ on the tangent space at a the [`Identity`}(@ref) and tangent
component vector ``X^i``, compute the equivalent vector representation
``X=X^i e_i**, where Einstein summation notation is used:

````math
∧ : X^i ↦ X^i e_i
````

For array manifolds, this converts a vector representation of the tangent
vector to an array representation. The [`vee`](@ref) map is the `hat` map's
inverse.
"""
function hat(
    M::AbstractDecoratorManifold,
    ::Identity{O},
    X,
) where {O<:AbstractGroupOperation}
    return get_vector_lie(M, X, VeeOrthogonalBasis())
end
function hat!(
    M::AbstractDecoratorManifold,
    Y,
    ::Identity{O},
    X,
) where {O<:AbstractGroupOperation}
    return get_vector_lie!(M, Y, X, VeeOrthogonalBasis())
end
function hat(M::AbstractManifold, e::Identity, ::Any)
    return throw(ErrorException("On $M there exsists no identity $e"))
end
function hat!(M::AbstractManifold, ::Any, e::Identity, ::Any)
    return throw(ErrorException("On $M there exsists no identity $e"))
end

@doc raw"""
    vee(M::AbstractManifold, p, X)

Given a basis $e_i$ on the tangent space at a point `p` and tangent
vector `X`, compute the vector components $X^i$, such that $X = X^i e_i$, where
Einstein summation notation is used:

````math
\vee : X^i e_i ↦ X^i
````

For array manifolds, this converts an array representation of the tangent
vector to a vector representation. The [`hat`](@ref) map is the `vee` map's
inverse.
"""
function vee(
    M::AbstractDecoratorManifold,
    ::Identity{O},
    X,
) where {O<:AbstractGroupOperation}
    return get_coordinates_lie(M, X, VeeOrthogonalBasis())
end
function vee!(
    M::AbstractDecoratorManifold,
    Y,
    ::Identity{O},
    X,
) where {O<:AbstractGroupOperation}
    return get_coordinates_lie!(M, Y, X, VeeOrthogonalBasis())
end
function vee(M::AbstractManifold, e::Identity, ::Any)
    return throw(ErrorException("On $M there exsists no identity $e"))
end
function vee!(M::AbstractManifold, ::Any, e::Identity, ::Any)
    return throw(ErrorException("On $M there exsists no identity $e"))
end

"""
    lie_bracket(G::AbstractDecoratorManifold, X, Y)

Lie bracket between elements `X` and `Y` of the Lie algebra corresponding to
the Lie group `G`, cf. [`IsGroupManifold`](@ref).

This can be used to compute the adjoint representation of a Lie algebra.
Note that this representation isn't generally faithful. Notably the adjoint
representation of 𝔰𝔬(2) is trivial.
"""
lie_bracket(G::AbstractDecoratorManifold, X, Y)
@trait_function lie_bracket(M::AbstractDecoratorManifold, X, Y)

@trait_function lie_bracket!(M::AbstractDecoratorManifold, Z, X, Y)

_action_order(p, q, ::LeftAction) = (p, q)
_action_order(p, q, ::RightAction) = (q, p)

@doc raw"""
    translate(G::AbstractDecoratorManifold, p, q, conv::ActionDirection=LeftAction()])

Translate group element $q$ by $p$ with the translation $τ_p$ with the specified
`conv`ention, either left ($L_p$) or right ($R_p$), defined as
```math
\begin{aligned}
L_p &: q ↦ p \circ q\\
R_p &: q ↦ q \circ p.
\end{aligned}
```
"""
translate(::AbstractDecoratorManifold, ::Any...)
@trait_function translate(
    G::AbstractDecoratorManifold,
    p,
    q,
    conv::ActionDirection=LeftAction(),
)
function translate(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    p,
    q,
    conv::ActionDirection,
)
    return compose(G, _action_order(p, q, conv)...)
end

@trait_function translate!(
    G::AbstractDecoratorManifold,
    X,
    p,
    q,
    conv::ActionDirection=LeftAction(),
)
function translate!(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    X,
    p,
    q,
    conv::ActionDirection,
)
    return compose!(G, X, _action_order(p, q, conv)...)
end

@doc raw"""
    inverse_translate(G::AbstractDecoratorManifold, p, q, conv::ActionDirection=LeftAction())

Inverse translate group element $q$ by $p$ with the inverse translation $τ_p^{-1}$ with the
specified `conv`ention, either left ($L_p^{-1}$) or right ($R_p^{-1}$), defined as
```math
\begin{aligned}
L_p^{-1} &: q ↦ p^{-1} \circ q\\
R_p^{-1} &: q ↦ q \circ p^{-1}.
\end{aligned}
```
"""
inverse_translate(::AbstractDecoratorManifold, ::Any...)
@trait_function inverse_translate(
    G::AbstractDecoratorManifold,
    p,
    q,
    conv::ActionDirection=LeftAction(),
)
function inverse_translate(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    p,
    q,
    conv::ActionDirection,
)
    return translate(G, inv(G, p), q, conv)
end

@trait_function inverse_translate!(
    G::AbstractDecoratorManifold,
    X,
    p,
    q,
    conv::ActionDirection=LeftAction(),
)
function inverse_translate!(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    X,
    p,
    q,
    conv::ActionDirection,
)
    return translate!(G, X, inv(G, p), q, conv)
end

@doc raw"""
    translate_diff(G::AbstractDecoratorManifold, p, q, X, conv::ActionDirection=LeftAction())

For group elements $p, q ∈ \mathcal{G}$ and tangent vector $X ∈ T_q \mathcal{G}$, compute
the action of the differential of the translation $τ_p$ by $p$ on $X$, with the specified
left or right `conv`ention. The differential transports vectors:
```math
(\mathrm{d}τ_p)_q : T_q \mathcal{G} → T_{τ_p q} \mathcal{G}\\
```
"""
translate_diff(::AbstractDecoratorManifold, ::Any...)
@trait_function translate_diff(
    G::AbstractDecoratorManifold,
    p,
    q,
    X,
    conv::ActionDirection=LeftAction(),
)
function translate_diff(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    p,
    q,
    X,
    conv::ActionDirection,
)
    Y = allocate_result(G, translate_diff, X, p, q)
    translate_diff!(G, Y, p, q, X, conv)
    return Y
end
@trait_function translate_diff!(
    G::AbstractDecoratorManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection=LeftAction(),
)

@doc raw"""
    inverse_translate_diff(G::AbstractDecoratorManifold, p, q, X, conv::ActionDirection=LeftAction())

For group elements $p, q ∈ \mathcal{G}$ and tangent vector $X ∈ T_q \mathcal{G}$, compute
the action on $X$ of the differential of the inverse translation $τ_p$ by $p$, with the
specified left or right `conv`ention. The differential transports vectors:
```math
(\mathrm{d}τ_p^{-1})_q : T_q \mathcal{G} → T_{τ_p^{-1} q} \mathcal{G}\\
```
"""
inverse_translate_diff(::AbstractDecoratorManifold, ::Any...)
@trait_function inverse_translate_diff(
    G::AbstractDecoratorManifold,
    p,
    q,
    X,
    conv::ActionDirection=LeftAction(),
)
function inverse_translate_diff(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return translate_diff(G, inv(G, p), q, X, conv)
end

@trait_function inverse_translate_diff!(
    G::AbstractDecoratorManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection=LeftAction(),
)
function inverse_translate_diff!(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    Y,
    p,
    q,
    X,
    conv::ActionDirection,
)
    return translate_diff!(G, Y, inv(G, p), q, X, conv)
end

@doc raw"""
    exp_lie(G, X)
    exp_lie!(G, q, X)

Compute the group exponential of the Lie algebra element `X`. It is equivalent to the
exponential map defined by the [`CartanSchoutenMinus`](@ref) connection.

Given an element $X ∈ 𝔤 = T_e \mathcal{G}$, where $e$ is the [`Identity`](@ref) element of
the group $\mathcal{G}$, and $𝔤$ is its Lie algebra, the group exponential is the map

````math
\exp : 𝔤 → \mathcal{G},
````
such that for $t,s ∈ ℝ$, $γ(t) = \exp (t X)$ defines a one-parameter subgroup with the
following properties:

````math
\begin{aligned}
γ(t) &= γ(-t)^{-1}\\
γ(t + s) &= γ(t) \circ γ(s) = γ(s) \circ γ(t)\\
γ(0) &= e\\
\lim_{t → 0} \frac{d}{dt} γ(t) &= X.
\end{aligned}
````

!!! note
    In general, the group exponential map is distinct from the Riemannian exponential map
    [`exp`](@ref).

For example for the [`MultiplicationOperation`](@ref) and either `Number` or `AbstractMatrix`
the Lie exponential is the numeric/matrix exponential.

````math
\exp X = \operatorname{Exp} X = \sum_{n=0}^∞ \frac{1}{n!} X^n.
````

Since this function also depends on the group operation, make sure to implement
the corresponding trait version `exp_lie(::TraitList{<:IsGroupManifold}, G, X)`.
"""
exp_lie(G::AbstractManifold, X)
@trait_function exp_lie(M::AbstractDecoratorManifold, X)
function exp_lie(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, X)
    q = allocate_result(G, exp_lie, X)
    return exp_lie!(G, q, X)
end

@trait_function exp_lie!(M::AbstractDecoratorManifold, q, X)

@doc raw"""
    log_lie(G, q)
    log_lie!(G, X, q)

Compute the Lie group logarithm of the Lie group element `q`. It is equivalent to the
logarithmic map defined by the [`CartanSchoutenMinus`](@ref) connection.

Given an element $q ∈ \mathcal{G}$, compute the right inverse of the group exponential map
[`exp_lie`](@ref), that is, the element $\log q = X ∈ 𝔤 = T_e \mathcal{G}$, such that
$q = \exp X$

!!! note
    In general, the group logarithm map is distinct from the Riemannian logarithm map
    [`log`](@ref).

    For matrix Lie groups this is equal to the (matrix) logarithm:

````math
\log q = \operatorname{Log} q = \sum_{n=1}^∞ \frac{(-1)^{n+1}}{n} (q - e)^n,
````

where $e$ here is the [`Identity`](@ref) element, that is, $1$ for numeric $q$ or the
identity matrix $I_m$ for matrix $q ∈ ℝ^{m × m}$.

Since this function also depends on the group operation, make sure to implement
either
* `_log_lie(G, q)` and `_log_lie!(G, X, q)` for the points not being the [`Identity`](@ref)
* the trait version `log_lie(::TraitList{<:IsGroupManifold}, G, e)`, `log_lie(::TraitList{<:IsGroupManifold}, G, X, e)` for own implementations of the identity case.
"""
log_lie(::AbstractDecoratorManifold, q)
@trait_function log_lie(G::AbstractDecoratorManifold, q)
function log_lie(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, q)
    return _log_lie(G, q)
end
function log_lie(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return zero_vector(G, identity_element(G))
end
# though identity was taken care of – as usual restart decorator dispatch
function _log_lie(G::AbstractDecoratorManifold, q)
    X = allocate_result(G, log_lie, q)
    return log_lie!(G, X, q)
end

@trait_function log_lie!(G::AbstractDecoratorManifold, X, q)
function log_lie!(::TraitList{<:IsGroupManifold}, G::AbstractDecoratorManifold, X, q)
    return _log_lie!(G, X, q)
end
function log_lie!(
    ::TraitList{<:IsGroupManifold{O}},
    G::AbstractDecoratorManifold,
    X,
    ::Identity{O},
) where {O<:AbstractGroupOperation}
    return zero_vector!(G, X, identity_element(G))
end

"""
    GroupExponentialRetraction{D<:ActionDirection} <: AbstractRetractionMethod

Retraction using the group exponential [`exp_lie`](@ref) "translated" to any point on the
manifold.

For more details, see
[`retract`](@ref retract(::GroupManifold, p, X, ::GroupExponentialRetraction)).

# Constructor

    GroupExponentialRetraction(conv::ActionDirection = LeftAction())
"""
struct GroupExponentialRetraction{D<:ActionDirection} <: AbstractRetractionMethod end

function GroupExponentialRetraction(conv::ActionDirection=LeftAction())
    return GroupExponentialRetraction{typeof(conv)}()
end

"""
    GroupLogarithmicInverseRetraction{D<:ActionDirection} <: AbstractInverseRetractionMethod

Retraction using the group logarithm [`log_lie`](@ref) "translated" to any point on the
manifold.

For more details, see
[`inverse_retract`](@ref inverse_retract(::GroupManifold, p, q ::GroupLogarithmicInverseRetraction)).

# Constructor

    GroupLogarithmicInverseRetraction(conv::ActionDirection = LeftAction())
"""
struct GroupLogarithmicInverseRetraction{D<:ActionDirection} <:
       AbstractInverseRetractionMethod end

function GroupLogarithmicInverseRetraction(conv::ActionDirection=LeftAction())
    return GroupLogarithmicInverseRetraction{typeof(conv)}()
end

direction(::GroupExponentialRetraction{D}) where {D} = D()
direction(::GroupLogarithmicInverseRetraction{D}) where {D} = D()

@doc raw"""
    retract(
        G::AbstractDecoratorManifold,
        p,
        X,
        method::GroupExponentialRetraction{<:ActionDirection},
    )

Compute the retraction using the group exponential [`exp_lie`](@ref) "translated" to any
point on the manifold.
With a group translation ([`translate`](@ref)) $τ_p$ in a specified direction, the
retraction is

````math
\operatorname{retr}_p = τ_p \circ \exp \circ (\mathrm{d}τ_p^{-1})_p,
````

where $\exp$ is the group exponential ([`exp_lie`](@ref)), and $(\mathrm{d}τ_p^{-1})_p$ is
the action of the differential of inverse translation $τ_p^{-1}$ evaluated at $p$ (see
[`inverse_translate_diff`](@ref)).
"""
function retract(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    p,
    X,
    method::GroupExponentialRetraction,
)
    conv = direction(method)
    Xₑ = inverse_translate_diff(G, p, p, X, conv)
    pinvq = exp_lie(G, Xₑ)
    q = translate(G, p, pinvq, conv)
    return q
end

function retract!(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    q,
    p,
    X,
    method::GroupExponentialRetraction,
)
    conv = direction(method)
    Xₑ = inverse_translate_diff(G, p, p, X, conv)
    pinvq = exp_lie(G, Xₑ)
    return translate!(G, q, p, pinvq, conv)
end

@doc raw"""
    inverse_retract(
        G::AbstractDecoratorManifold,
        p,
        X,
        method::GroupLogarithmicInverseRetraction{<:ActionDirection},
    )

Compute the inverse retraction using the group logarithm [`log_lie`](@ref) "translated"
to any point on the manifold.
With a group translation ([`translate`](@ref)) $τ_p$ in a specified direction, the
retraction is

````math
\operatorname{retr}_p^{-1} = (\mathrm{d}τ_p)_e \circ \log \circ τ_p^{-1},
````

where $\log$ is the group logarithm ([`log_lie`](@ref)), and $(\mathrm{d}τ_p)_e$ is the
action of the differential of translation $τ_p$ evaluated at the identity element $e$
(see [`translate_diff`](@ref)).
"""
function inverse_retract(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    p,
    q,
    method::GroupLogarithmicInverseRetraction,
)
    conv = direction(method)
    pinvq = inverse_translate(G, p, q, conv)
    Xₑ = log_lie(G, pinvq)
    return translate_diff(G, p, Identity(G), Xₑ, conv)
end

function inverse_retract!(
    ::TraitList{<:IsGroupManifold},
    G::AbstractDecoratorManifold,
    X,
    p,
    q,
    method::GroupLogarithmicInverseRetraction,
)
    conv = direction(method)
    pinvq = inverse_translate(G, p, q, conv)
    Xₑ = log_lie(G, pinvq)
    return translate_diff!(G, X, p, Identity(G), Xₑ, conv)
end

@doc raw"""
    get_vector_lie(G::AbstractDecoratorManifold, a, B::AbstractBasis)

Reconstruct a tangent vector from the Lie algebra of `G` from cooordinates `a` of a basis `B`.
This is similar to calling [`get_vector`](@ref) at the `p=`[`Identity`](@ref)`(G)`.
"""
function get_vector_lie(G::AbstractManifold, X, B::AbstractBasis)
    return get_vector(base_manifold(G), identity_element(G), X, B)
end
function get_vector_lie!(G::AbstractManifold, Y, X, B::AbstractBasis)
    return get_vector!(base_manifold(G), Y, identity_element(G), X, B)
end

@doc raw"""
    get_coordinates_lie(G::AbstractManifold, X, B::AbstractBasis)

Get the coordinates of an element `X` from the Lie algebra og `G` with respect to a basis `B`.
This is similar to calling [`get_coordinates`](@ref) at the `p=`[`Identity`](@ref)`(G)`.
"""
function get_coordinates_lie(G::AbstractManifold, X, B::AbstractBasis)
    return get_coordinates(base_manifold(G), identity_element(G), X, B)
end
function get_coordinates_lie!(G::AbstractManifold, a, X, B::AbstractBasis)
    return get_coordinates!(base_manifold(G), a, identity_element(G), X, B)
end


text:
###### readtimearray ############

function readtimearray(source; delim::Char=',', meta=nothing, format::AbstractString="")
    cfile = readdlm(source, delim)

    # remove empty lines if any
    inoempty = find(s -> length(s) > 2, cfile[:,1])
    cfile = cfile[inoempty,:]

    time  = cfile[2:end,1]
    if length(time[1]) < 11
        # assuming Date not DateTime
        format == "" ?
        tstamps = Date[Date(t) for t in time] :
        tstamps = Date[Date(t, format) for t in time]
    else
        format == "" ?
        tstamps = DateTime[DateTime(t) for t in time] :
        tstamps = DateTime[DateTime(t, format) for t in time]
    end

    vals   = insertNaN(cfile[2:end, 2:end])
    cnames = String[]
    for c in cfile[1, 2:end]
        push!(cnames, string(c))
    end
    TimeArray(tstamps, vals, cnames, meta)
end

function insertNaN{N}(aa::Array{Any,N})
    for i in 1:size(aa,1)
        for j in 1:size(aa,2)
            if !isa(aa[i,j], Real)
                aa[i,j] = NaN
            end
        end
    end
    convert(Array{Float64,N},aa)
end

function writetimearray(ta::TimeArray, fname::AbstractString)

    open(fname, "w") do io

        strvals = join(ta.colnames, ",")
        write(io, string("Timestamp,", strvals, "\n"))

        for i in eachindex(ta.timestamp)
            strvals = replace(join(ta.values[i, :], ","), "NaN", "")
            write(io, string(ta.timestamp[i], ",", strvals, "\n"))
        end #for

    end

end #writetimearray


text:
module ChebyshevQuantum

using Plots
using ForwardDiff
using Polynomials
using LinearAlgebra

greet() = print("Hello World!")

include("Interp.jl")
include("EdgeDetect.jl")
include("Cheb.jl")
include("Op.jl")
include("SolveEig.jl")
include("SolveDiff.jl")

end # module


text:
using BadAD
using ForwardDiff
using Test

@testset "BadAD.jl" begin

    a = b = Dual(1, 1)

    S(x) = 1/(1 + exp(-x))
    dS(x) = exp(-x)/(1 + exp(-x))^2

    f(x) = x^2
    g(x) = sin(x)
    h(x) = f(x)/g(x)

    f2(x) = 3^x

    @test -1 * (a + b) == -a + -b
    @test derivative(g ∘ f, 3.14) == 2*3.14*cos(3.14^2)
    @test derivative(h, 3) == (2*3*sin(3) - 3^2 * cos(3))/(sin(3)^2)
    @test derivative(S, 1) == dS(1)
    @test ForwardDiff.derivative(x->3.5^x, 3) ≈ derivative(x->3.5^x, 3)

    f3(x) = SA[log(x[1])]
end


text:
using ModelingToolkit, OrdinaryDiffEq, Test

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

@named lorenz1 = ODESystem(eqs)
@named lorenz2 = ODESystem(eqs)

@parameters γ
@variables a(t) α(t)
connections = [0 ~ lorenz1.x + lorenz2.y + a*γ,
               α ~ 2lorenz1.x + a*γ]
@named sys = ODESystem(connections,t,[a,α],[γ],systems=[lorenz1,lorenz2])
sys_simplified = structural_simplify(sys)

u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0,
      a => 2.0]

p  = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8/3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8/3,
      γ => 2.0]

tspan = (0.0,100.0)
prob = ODEProblem(sys_simplified,u0,tspan,p)
sol = solve(prob,Rodas4())

@test_throws UndefVarError sol[b]
@test_throws UndefVarError sol[b, 1]
@test_throws UndefVarError sol[b, 1:5]
@test_throws UndefVarError sol[b, [1,2,3]]
@test_throws ErrorException sol['a']
@test_throws ErrorException sol['a',  1]
@test_throws ErrorException sol['a',  1:5]
@test_throws ErrorException sol['a',  [1,2,3]]

@test sol[a] isa AbstractVector
@test sol[a,  1] isa Real
@test sol[a,  1:5] isa AbstractVector
@test sol[a,  [1,2,3]] isa AbstractVector

@test sol[1] isa AbstractVector
@test sol[1:2] isa AbstractArray
@test sol[[1,2]] isa AbstractArray

@test sol[lorenz1.x] isa Vector
@test sol[lorenz1.x,2] isa Float64
@test sol[lorenz1.x,:] isa Vector
@test sol[t] isa Vector
@test sol[t,2] isa Float64
@test sol[t,:] isa Vector
@test length(sol[lorenz1.x,1:5]) == 5
@test sol[α] isa Vector
@test sol[α,3] isa Float64
@test length(sol[α,5:10]) == 6

# Check if indexing using variable names from interpolated solution works
interpolated_sol = sol(0.0:1.0:10.0)
@test interpolated_sol[α] isa Vector
@test interpolated_sol[α,:] isa Vector
@test interpolated_sol[α,2] isa Float64
@test length(interpolated_sol[α,1:5]) == 5
@test interpolated_sol[α] ≈ 2interpolated_sol[lorenz1.x] .+ interpolated_sol[a].*2.0
@test collect(interpolated_sol[t]) isa Vector
@test collect(interpolated_sol[t,:]) isa Vector
@test interpolated_sol[t,2] isa Float64
@test length(interpolated_sol[t,1:5]) == 5

sol1 = sol(0.0:1.0:10.0)
@test sol1.u isa Vector
@test first(sol1.u) isa Vector
@test length(sol1.u) == 11
@test length(sol1.t) == 11

sol2 = sol(0.1)
@test sol2 isa Vector
@test length(sol2) == length(states(sys_simplified))
@test first(sol2) isa Real

sol3 = sol(0.0:1.0:10.0, idxs=[lorenz1.x, lorenz2.x])
@test sol3.u isa Vector
@test first(sol3.u) isa Vector
@test length(sol3.u) == 11
@test length(sol3.t) == 11
@test collect(sol3[t]) ≈ sol3.t
@test collect(sol3[t, 1:5]) ≈ sol3.t[1:5]
@test_throws ErrorException sol(0.0:1.0:10.0, idxs=[lorenz1.x, 1])

sol4 = sol(0.1, idxs=[lorenz1.x, lorenz2.x])
@test sol4 isa Vector
@test length(sol4) == 2
@test first(sol4) isa Real
@test_throws ErrorException sol(0.1, idxs=[lorenz1.x, 1])

sol5 = sol(0.0:1.0:10.0, idxs=lorenz1.x)
@test sol5.u isa Vector
@test first(sol5.u) isa Real
@test length(sol5.u) == 11
@test length(sol5.t) == 11
@test collect(sol5[t]) ≈ sol3.t
@test collect(sol5[t, 1:5]) ≈ sol3.t[1:5]
@test_throws ErrorException sol(0.0:1.0:10.0, idxs=1.2)

sol6 = sol(0.1, idxs=lorenz1.x)
@test sol6 isa Real
@test_throws ErrorException sol(0.1, idxs=1.2)

sol7 = sol(0.0:1.0:10.0, idxs=[2,1])
@test sol7.u isa Vector
@test first(sol7.u) isa Vector
@test length(sol7.u) == 11
@test length(sol7.t) == 11
@test collect(sol7[t]) ≈ sol3.t
@test collect(sol7[t, 1:5]) ≈ sol3.t[1:5]

sol8 = sol(0.1, idxs=[2,1])
@test sol8 isa Vector
@test length(sol8) == 2
@test first(sol8) isa Real

sol9 = sol(0.0:1.0:10.0, idxs=2)
@test sol9.u isa Vector
@test first(sol9.u) isa Real
@test length(sol9.u) == 11
@test length(sol9.t) == 11
@test collect(sol9[t]) ≈ sol3.t
@test collect(sol9[t, 1:5]) ≈ sol3.t[1:5]

sol10 = sol(0.1, idxs=2)
@test sol10 isa Real

#=
using Plots
plot(sol,vars=(lorenz2.x,lorenz2.z))
plot(sol,vars=(α,lorenz2.z))
plot(sol,vars=(lorenz2.x,α))
plot(sol,vars=α)
plot(sol,vars=(t,α))
=#


text:
#=
  本文件配置个性化参数，具体参见 initAlarmSegparas! 中的注释
=#

function initParaTuple( ;segID::Int,StaticThreshold::Int , DynamicThreshold_CalTimePeriod::UInt,DynamicThresHold_calCFactor::Int,
                        microMergeTimePeriod::Int,
                        microFilterTimePeriod::Int,
                        microMergeSpaceLength::UInt,
                        microFilterSpaceLength::UInt,
                        macroMaxSpaceLength::UInt,
                        macroMinSpaceLength::UInt,
                        macroMaxTimeInterval::UInt,
                        macroMinTimeInterval::UInt,
            )

end

function initAlarmSegparas!(SegAlarmParasDictObj::Dict )   #####################
    SegParaNameTuple1 =  ( segID = 0 , microtime=Seg_MicroTimeParaStruct() ,
    microSpace=Seg_MicroSpaceParaStruct(),
    macroSpace= Seg_MacroSpaceParaStruct(),
    macroTime= Seg_MacroTimeParaStruct(),
    )
    # 静态阈值
    SegParaNameTuple1.microtime.StaticThreshold = Int(2200) # if equal to the typeMAxInt, the Static TH  does not be in  effect
    #动态阈值
    SegParaNameTuple1.microtime.DynamicThreshold_CalTimePeriod = UInt( 60 ) ## in 0.426 seconds  # calculated from Seconds.seconds ,may be any time period ,between 10seconds and the cached data period(about 10 minutes above )
    #= the  TH is: the Calculated Avg * factor / 100
        如果动态阈值不适用该点，则DynamicThresHold_calCFactor = typemax(Int)
    =#
    SegParaNameTuple1.microtime.DynamicThresHold_calCFactor=Int( typemax(Int) )

    #微观过滤和合并 时间
    SegParaNameTuple1.microtime.microMergeTimePeriod  = Int(2) # in 0.426 seconds # calculated from :seconds * 0.426
    SegParaNameTuple1.microtime.microFilterTimePeriod = Int(2) # in 0.426 seconds # same as above

    #微观过滤和合并 空间
    SegParaNameTuple1.microSpace.microMergeSpaceLength = UInt(1) # in Segment or  grain. Calculated from the space length and the segPermeters
    SegParaNameTuple1.microSpace.microFilterSpaceLength = UInt(3) # in Segment   or  grain

    #宏观 空间 限制范围
    SegParaNameTuple1.macroSpace.macroMaxSpaceLength = UInt(100)# in Segment or  grain. Calculated from the space length and the segPermeters
    SegParaNameTuple1.macroSpace.macroMinSpaceLength = UInt(4)# in Segment   or  grain
    #宏观 时间  限制范围
    SegParaNameTuple1.macroTime.macroMaxTimeInterval = UInt(30) # in 0.426 seconds
    SegParaNameTuple1.macroTime.macroMinTimeInterval = UInt(15) # in  0.426 seconds

    SegAlarmParasDictObj[SegParaNameTuple1.segID ] = SegParaNameTuple1
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
@testset "repartition_coalesce" begin

# test repartition and coalesce
sc = SparkContext(master="local")

nums1 = parallelize(sc, 1:30)
@test num_partitions(nums1) == 1
nums3 = repartition(nums1, 3)
@test num_partitions(nums3) == 3
nums2 = Spark.coalesce(nums3, 2)
@test num_partitions(nums2) == 2
nums2 = Spark.coalesce(nums3, 2; shuffle=true)
@test num_partitions(nums2) == 2

pnums1 = cartesian(nums1, nums1)
@test num_partitions(pnums1) == 1
pnums3 = repartition(pnums1, 3)
@test num_partitions(pnums3) == 3
pnums2 = Spark.coalesce(pnums3, 2)
@test num_partitions(pnums2) == 2
pnums2 = Spark.coalesce(pnums3, 2; shuffle=true)
@test num_partitions(pnums2) == 2

close(sc)

end


text:
"""
`MaternKernel([ρ=1.0,[ν=1.0]])`
The matern kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = 2^{1-ν}/Γ(ν)*(√(2ν)‖x-y‖)^ν K_ν(√(2ν)‖x-y‖)
```
For `ν=n+1/2, n=0,1,2,...` it can be simplified and you should instead use [`ExponentialKernel`](@ref) for `n=0`, [`Matern32Kernel`](@ref), for `n=1`, [`Matern52Kernel`](@ref) for `n=2` and [`SqExponentialKernel`](@ref) for `n=∞`.
"""
struct MaternKernel{T,Tr,Tν<:Real} <: Kernel{T,Tr}
    transform::Tr
    metric::Euclidean
    ν::Tν
    function MaternKernel{T,Tr,Tν}(transform::Tr,ν::Tν) where {T,Tr<:Transform,Tν<:Real}
        return new{T,Tr,Tν}(transform,Euclidean(),ν)
    end
end

function MaternKernel(ρ::T₁=1.0,ν::T₂=1.5) where {T₁<:Real,T₂<:Real}
    @check_args(MaternKernel, ν, ν > zero(T₂), "ν > 0")
    MaternKernel{T₁,ScaleTransform{Base.RefValue{T₁}},T₂}(ScaleTransform(ρ),ν)
end

function MaternKernel(ρ::A,ν::T=1.5) where {A<:AbstractVector{<:Real},T<:Real}
    @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
    MaternKernel{eltype(A),ScaleTransform{A},T}(ScaleTransform(ρ),ν)
end

function MaternKernel(t::Tr,ν::T=1.5) where {Tr<:Transform,T<:Real}
    @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
    MaternKernel{eltype(t),Tr,T}(t,ν)
end

params(k::MaternKernel) = (params(transform(k)),k.ν)

@inline kappa(κ::MaternKernel, d::Real) = iszero(d) ? one(d) : exp((1.0-κ.ν)*logtwo-lgamma(κ.ν) + κ.ν*log(sqrt(2κ.ν)*d)+log(besselk(κ.ν,sqrt(2κ.ν)*d)))

"""
`Matern32Kernel([ρ=1.0])`
The matern 3/2 kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = (1+√(3)ρ‖x-y‖)exp(-√(3)ρ‖x-y‖)
```
"""
struct Matern32Kernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::Euclidean
    function Matern32Kernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
        return new{T,Tr}(transform,Euclidean())
    end
end

@inline kappa(κ::Matern32Kernel, d::T) where {T<:Real} = (1+sqrt(3)*d)*exp(-sqrt(3)*d)

"""
`Matern52Kernel([ρ=1.0])`
The matern 5/2 kernel is an isotropic Mercer kernel given by the formula:
```
    κ(x,y) = (1+√(5)ρ‖x-y‖ + 5ρ²‖x-y‖^2/3)exp(-√(5)ρ‖x-y‖)
```
"""
struct Matern52Kernel{T,Tr} <: Kernel{T,Tr}
    transform::Tr
    metric::Euclidean
    function Matern52Kernel{T,Tr}(transform::Tr) where {T,Tr<:Transform}
        return new{T,Tr}(transform,Euclidean())
    end
end

@inline kappa(κ::Matern52Kernel, d::Real) where {T} = (1+sqrt(5)*d+5*d^2/3)*exp(-sqrt(5)*d)


text:
"""
Quadratic Programming Optimal Enriched-Set Cover problem.
"""
struct QuadraticCoverProblem <: AbstractCoverProblem{Float64}
    params::CoverParams

    var2set::Vector{Int}
    var_scores::Vector{Float64}
    varXvar_scores::Matrix{Float64}

    function QuadraticCoverProblem(params::CoverParams,
            var2set::Vector{Int},
            var_scores::Vector{Float64},
            varXvar_scores::Matrix{Float64}
    )
        length(var2set) == length(var_scores) ==
        size(varXvar_scores, 1) == size(varXvar_scores, 2) ||
            throw(ArgumentError("var2set, var_scores and varXvar_scores sizes do not match"))
        new(params, var2set, var_scores, varXvar_scores)
    end
end

struct QuadraticOptimizerParams <: AbstractOptimizerParams{QuadraticCoverProblem}
    solver_factory #=::JuMP.OptimizerFactory=#

    QuadraticOptimizerParams(;
        solver_factory#=::JuMP.OptimizerFactory=# = with_default_quadratic_solver(),
        #iterations::Int = 100,
        #ini_weights::Vector{Float64}
    ) = new(solver_factory)
end

with_quadratic_solver(opt_params::QuadraticOptimizerParams) = opt_params.solver_factory

function QuadraticCoverProblem(mosaic::MaskedSetMosaic, params::CoverParams = CoverParams())
    v2set = var2set(mosaic)
    QuadraticCoverProblem(params, v2set,
                          var_scores(mosaic, v2set, params),
                          varXvar_scores(mosaic, v2set, params, true))
end

"""
Construct JuMP quadratic minimization model with linear contraints for the given OESC problem.
"""
function opt_model(problem::QuadraticCoverProblem, opt_params::QuadraticOptimizerParams)
    m = JuMP.Model(with_quadratic_solver(opt_params))
    nw = nvars(problem)
    #= @variable =# JuMP.variable(m, 0.0 <= w[1:nw] <= 1.0)
    #= @objective =# JuMP.objective(m, :Min, dot(vec(problem.var_scores), w) - dot(w, problem.varXvar_scores * w))
    return m
end

"""
Score of the OESC coverage.

* `w` probabilities of the sets being covered
"""
function score(problem::QuadraticCoverProblem, w::Vector{Float64})
    #pen = fix_cover_weights!(w) # FIXME throw an error?
    dot(problem.var_scores - problem.varXvar_scores * w, w)
end

aggscore(problem::QuadraticCoverProblem, w::Vector{Float64}) = score(w, problem)

"""
    optimize(problem::AbstractCoverProblem, method_params::QuadraticMethodParams)
    optimize(problem::AbstractCoverProblem, [method_kwparams...])

Optimize the cover problem.
"""
function optimize(problem::QuadraticCoverProblem,
                  opt_params::QuadraticOptimizerParams = QuadraticOptimizerParams(),
                  ini_weights::Vector{Float64} = rand(nvars(problem)) # unused
)
    if nvars(problem) == 0
        return CoverProblemResult(problem.var2set, Vector{Float64}(), Vector{Float64}(), 0.0, 0.0)
    end

    # Perform the optimization
    #try
    # using JuMP
    m = opt_model(problem, opt_params)
    solve(m)
    w = copy(getvalue(getindex(m, :w)))
    # remove small non-zero probabilities due to optimization method errors
    for i in eachindex(w)
        @inbounds if w[i] < problem.params.min_weight
            w[i] = 0.0
        end
    end
    s = getobjectivevalue(m)
    return CoverProblemResult(problem.var2set, w, problem.var_scores .* w, s, s)
    #catch x
    #    @warn("Exception in optimize(CoverProblem): $x")
    #    return nothing
    #end
end

function exclude_vars(problem::QuadraticCoverProblem,
                      vars::AbstractVector{Int};
                      penalize_overlaps::Bool = true)
    varmask = fill(true, nvars(problem))
    varmask[vars] = false
    v_scores = problem.var_scores[varmask]
    if penalize_overlaps
        # penalize overlapping sets
        penalty_weights = fill(0.0, nvars(problem))
        penalty_weights[vars] = 1.0
        v_scores .-= (problem.varXvar_scores * penalty_weights)[varmask]
    end
    return QuadraticCoverProblem(problem.params,
                problem.var2set[varmask],
                v_scores, problem.varXvar_scores[varmask, varmask])
end


text:
include("quadrature.jl")
include("../computation/derivatives.jl")
include("../transforms.jl")

"""
``getTensorGrid(train, test, priorθ, priorλ, nodesWeightsθ, nodesWeightsλ, transform, quadtype)``

* quadtype can either be \"Gaussian\" or \"Turan\" (N.B. this refers to integration of θ variable, λ always integrated out using Gaussian quadrature)

Works by defining a function ``f`` from ``R^k`` to ``Mat(n, n)``, such that ``f(z_0)_{ij} = p(z_0|z, θ_i, λ_j)``,
where ``i`` and ``j`` range over the meshgrids over ``θ`` and ``λ``. Takes the dot product of this
matrix function with a weights matrix, ``w_ij = p(θ_i, λ_j|z)`` to obtain PDF and CDF functions.

```julia
getTensorGrid(train, test, priorθ, priorλ, nodesWeightsθ, nodesWeightsλ, boxCoxObj, "Gaussian") #uses BoxCox family of nonlinear transforms and Gaussian quadrature
```

"""
function getTensorGrid(integrand, train::trainingData{T1, T2}, test::testingData{T1}, priorθ, priorλ, nodesWeightsθ::nodesWeights, nodesWeightsλ::nodesWeights, transform, quadtype::String) where T1 <:Array{Float64} where T2<:Array{Float64}
    nodesθ = nodesWeightsθ.nodes
    nodesλ = nodesWeightsλ.nodes
    weightsθ = nodesWeightsθ.weights
    weightsλ = nodesWeightsλ.weights
    X = train.X; X0 = test.X0; z = train.z; n = size(X, 1); p = size(X, 2)
    l1 = length(nodesθ); l2 = length(nodesλ); l3 = size(nodesWeightsθ.weights, 2)
    #initializes buffer of theta-dependent quantities
    function func_fixed(θ::Float64)
        return funcθ(θ, train, test, quadtype)
    end
    #precompute buffer for θ-values
    theta_param_list = Array{Union{θ_params{Array{Float64, 2}, Cholesky{Float64,Array{Float64, 2}}}, θ_param_derivs{Array{Float64, 2}, Cholesky{Float64,Array{Float64, 2}}}}}(undef, l1)
    for i=1:l1
        theta_param_list[i] = func_fixed(nodesθ[i])
    end

    #initialize tensor grid of weights
    n1 = length(nodesθ); n2 = length(nodesλ); n3 = size(weightsθ, 2);
    weightsTensorGrid = zeros(n1, n2, n3) #tensor grid of weights
    for i = 1:n1
        for j = 1:n2
            funcs = posterior_theta(nodesθ[i], nodesλ[j], priorθ, priorλ, train, test, boxCoxObj, theta_param_list[i], quadtype)
            for k = 1:n3
                weightsTensorGrid[i, j, k] = weightsθ[i, k]*weightsλ[j]*funcs[k]
            end
        end
    end

    if quadtype == "Gaussian"
    #|Σθ|, qtilde, and Jac(z) are prone to overflow and underflow. As such, we work strictly with their exponents,
    #and shift them so the largest sum-of-exponents (exponentiating then gives 1). After that, we can safely
    #normalize the weight so that they sum to 1.
        #compute exponents of qtilde^(-(n-p)/2)
        qvals = zeros(n1, n2)
        g = transform.f; dg = transform.df
        gλvals = zeros(length(z), n2)
        lmbda = λ -> g(z, λ)
        for i = 1:n2
            gλvals[:, i] = lmbda(nodesλ[i])
        end
        for i = 1:n1
            for j = 1:n2
                #compute auxiliary quantities
                choleskyXΣX = theta_param_list[i].choleskyXΣX
                choleskyΣθ = theta_param_list[i].choleskyΣθ
                gλz = gλvals[:, j]
                βhat = choleskyXΣX\(X'*(choleskyΣθ\gλz))
                qtilde = (expr = gλz-X*βhat; expr'*(choleskyΣθ\expr))
                qvals[i, j] = -(n-p)/2 * log(qtilde)
            end
        end
        qTensorGrid = repeat(qvals, 1, 1, n3)

        #compute exponents of |Σθ|^(-1/2)
        detvals = zeros(n1, 1)
        for i = 1:n1
            choleskyΣθ = theta_param_list[i].choleskyΣθ
            detvals[i] = sum(log.(diag(choleskyΣθ.U)))
        end
        #we multiply by -1.0 instead of -0.5, because the determinant is the
        #product of the squares of the diagonal entries of choleskyΣθ.U
        detTensorGrid = repeat((- 1.0 .* detvals), 1, n2, n3)

        #compute exponents of (Jac(z)^(1-p/n)
        jacvals = zeros(1, n2)
        df = transform.df
        for i = 1:n2
            jacvals[i] = sum(log.(abs.(map( x-> df(x, nodesλ[i]), z))))
        end
        jacTensorGrid = repeat((1-p/n) .* jacvals, n1, 1, n3)

        powerGrid =  qTensorGrid + detTensorGrid + jacTensorGrid #sum of exponents
        powerGrid = exp.(powerGrid .- maximum(powerGrid)) #linear scaling
        powerGrid = powerGrid .* weightsTensorGrid
        weightsTensorGrid =  powerGrid/sum(powerGrid) #normalized grid of weights

    elseif quadtype== "Turan"
        weightsTensorGrid =  weightsTensorGrid/sum(weightsTensorGrid) #normalized grid of weights
    end

    tgridpdf = Array{Any, 3}(undef, l1, l2, l3)
    #tensor grid for CDF P(z0|theta, lambda, z)
    tgridcdf = Array{Any, 3}(undef, l1, l2, l3)

    for i = 1:l1
        for j = 1:l2
            funcs = integrand(nodesθ[i], nodesλ[j], train, test, boxCoxObj, theta_param_list[i], quadtype)
            for k = 1:l3
                #tgrid[i, j, k] = funcs[k]
                tgridpdf[i, j, k] = funcs[1][k]
                tgridcdf[i, j, k] = funcs[2][k]
            end
        end
    end
    function evalTgrid_pdf(z0)
        res = Array{Float64, 3}(undef, l1, l2, l3)
        for i=1:l1
            for j = 1:l2
                for k =1:l3
                    #println("tgridpdf(z0): ", tgridpdf[i, j, k](z0))
                    res[i, j, k] = tgridpdf[i, j, k](z0)[1]
                end
            end
        end
        return res
    end
    function evalTgrid_cdf(z0)
        res = Array{Float64, 3}(undef, l1, l2, l3)
        for i=1:l1
            for j = 1:l2
                for k =1:l3
                    res[i, j, k] = tgridcdf[i, j, k](z0)
                end
            end
        end
        return res
    end
    #product rule loop for Turan quadrature (k>1). In the k=1 case, this reduces to taking the dot
    # product between the tensor grid and weights grid.
    function pdf(z0)
        grid = evalTgrid_pdf(z0)
        res = 0
        for i = 1:n3
            for j = 1:i
                res = binomial(i, j) * dot(grid[:, :, j],  weightsTensorGrid[:, :, n3-j+1])
            end
        end
        return res
    end
    function cdf(z0)
        grid = evalTgrid_cdf(z0)
        res = 0
        for i = 1:n3
            for j = 1:i
                res = binomial(i, j) * dot(grid[:, :, j], weightsTensorGrid[:, :, n3-j+1])
            end
        end
        return res
    end
    return (pdf, cdf)
    #return z0 -> dot(evalTgrid_pdf(z0), weightsTensorGrid), z0 -> dot(evalTgrid_cdf(z0), weightsTensorGrid)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../../test/useful_routines.jl")

function time_step_size(levelsetspeed, dx; CFL = 0.5)
    s = maximum(abs.(levelsetspeed))
    return CFL * dx / s
end

function step_interface(
    levelset,
    levelsetspeed,
    dt,
    paddedmesh,
    tol,
    boundingradius,
)
    mesh = CutCellDG.background_mesh(levelset)
    cutmesh = CutCellDG.CutMesh(mesh, levelset)

    refseedpoints, refseedcellids =
        CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
    spatialseedpoints =
        CutCellDG.map_to_spatial(refseedpoints, refseedcellids, mesh)

    paddedlevelset = CutCellDG.BoundaryPaddedLevelSet(
        paddedmesh,
        spatialseedpoints,
        refseedcellids,
        levelset,
        tol,
        boundingradius,
    )

    return CutCellDG.step_first_order_levelset(
        paddedlevelset,
        levelsetspeed,
        dt,
    )
end

function signed_distance_coefficients(levelset, mesh, tol, boundingradius)
    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    refseedpoints, refseedcellids =
        CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
    spatialseedpoints =
        CutCellDG.map_to_spatial(refseedpoints, refseedcellids, mesh)
    nodalcoordinates = CutCellDG.nodal_coordinates(mesh)

    signeddistance = CutCellDG.distance_to_zero_levelset(
        nodalcoordinates,
        spatialseedpoints,
        refseedcellids,
        levelset,
        tol,
        boundingradius,
    )

    return signeddistance
end

function run_time_steps(levelset, levelsetspeed, paddedmesh, dt, numiter)
    mesh = CutCellDG.background_mesh(levelset)
    dx = CutCellDG.element_size(mesh)
    polyorder = CutCellDG.order(levelset)
    tol = 1e-2(minimum(dx))^(polyorder + 1)
    boundingradius = 1.5maximum(dx)

    numnodes = CutCellDG.number_of_nodes(mesh)
    newcoeffs = zeros(numnodes)
    itercoeffs = zeros(numnodes, numiter + 1)
    itercoeffs[:, 1] = CutCellDG.coefficients(levelset)

    for iter = 1:numiter
        newcoeffs = step_interface(
            levelset,
            levelsetspeed,
            dt,
            paddedmesh,
            tol,
            boundingradius,
        )
        CutCellDG.update_coefficients!(levelset, newcoeffs)
        dist = signed_distance_coefficients(levelset, mesh, tol, boundingradius)
        CutCellDG.update_coefficients!(levelset, dist)
        newcoeffs = dist
    end
    return newcoeffs
end

function grid_range(mesh)
    x0 = CutCellDG.reference_corner(mesh)
    w = CutCellDG.mesh_widths(mesh)
    nfmside = CutCellDG.nodes_per_mesh_side(mesh)

    x = range(x0[1], stop = x0[1] + w[1], length = nfmside[1])
    y = range(x0[2], stop = x0[2] + w[2], length = nfmside[2])

    return x, y
end

function error_for_nelmts(nelmts)
    x0 = [0.0, 0.0]
    L, W = 1.0, 1.0
    polyorder = 2
    numqp = required_quadrature_order(polyorder)

    xc = [0.5, 0.5]
    radius = 0.45
    speed = -1.0
    stoptime = 0.25

    basis = LagrangeTensorProductBasis(2, polyorder)
    dim,nf = size(interpolation_points(basis))
    quad = tensor_product_quadrature(2, numqp)
    mesh = CutCellDG.CGMesh(x0, [L, W], [nelmts, nelmts], nf)
    paddedmesh = CutCellDG.BoundaryPaddedMesh(mesh, 1)
    numnodes = CutCellDG.number_of_nodes(mesh)
    dx = minimum(CutCellDG.grid_size(paddedmesh))

    levelset = CutCellDG.LevelSet(
        x -> -circle_distance_function(x, xc, radius),
        mesh,
        basis,
    )
    levelsetspeed = speed * ones(numnodes)
    dt = time_step_size(levelsetspeed, dx)
    @assert isinteger(stoptime / dt)
    numiter = round(Int, stoptime / dt)
    itercoeffs =
        run_time_steps(levelset, levelsetspeed, paddedmesh, dt, numiter)

    finalradius = radius - stoptime
    err = uniform_mesh_L2_error(
        itercoeffs',
        x -> -circle_distance_function(x, xc, finalradius)[1],
        basis,
        quad,
        mesh,
    )[1]

    normalizer = integral_norm_on_uniform_mesh(
        x -> -circle_distance_function(x, xc, finalradius)[1],
        quad,
        mesh,
        1,
    )[1]
    normalizederr = err / normalizer

    return normalizederr
end

powers = 2:5
nelmts = 2 .^powers
err = error_for_nelmts.(nelmts)
dx = 1.0 ./ nelmts
rate = convergence_rate(dx,err)

@test all(rate .> 0.9)


text:
using HTTP
using LightXML

headers = ["Ocp-Apim-Subscription-Key" => ENV["MYKEY"]]
r = HTTP.request("POST", string("https://", ENV["MYREGION"], ".api.cognitive.microsoft.com/sts/v1.0/issueToken"), headers)
token = String(r.body)
headers = ["Authorization" => string("Bearer ", token),
           "Content-Type" => "application/ssml+xml",
           "X-Microsoft-OutputFormat" => "riff-24khz-16bit-mono-pcm",
           "User-Agent" => "YOUR_RESOURCE_NAME"]

lang = "en-US"
xdoc = XMLDocument()
xroot = create_root(xdoc, "speak")
set_attributes(xroot, Dict("version" => "1.0",
                            "xmlns" => "http://www.w3.org/2001/10/synthesis",
                            "xmlns:mstts" => "https://www.w3.org/2001/mstts",
                            "xml:lang" => lang))
xs = new_child(xroot, "voice")
set_attributes(xs, Dict("xml:lang" => lang,
                        "xml:gender" => "Female",
                        "name" => "en-US-AriaRUS"))
add_text(xs, "123, 456. 789? 0!")

io = open("synthesized.wav", "w")
r = HTTP.request("POST", string("https://", ENV["MYREGION"], ".tts.speech.microsoft.com/cognitiveservices/v1"), headers, string(xdoc), response_stream=io)
close(io)
println("saved result to synthesized.wav")

max_stars_repo_name: boltomli/Cognitive-Speech-TTS

max_issues_repo_name: boltomli/Cognitive-Speech-TTS

max_forks_repo_name: boltomli/Cognitive-Speech-TTS


text:
using BenchmarkTools
using SimplexQuad

# @fastmath function integrate2d(f, X, W)
#     if size(X, 2) == 2
#         @inbounds sum(W[i] * f(X[i,1], X[i,2]) for i in 1:length(W))
#     elseif size(X, 2) == 2
#         @inbounds sum(W[i] * f(X[i,1], X[i,2], X[i,3]) for i in 1:length(W))
#     else
#         @error "not implemented"
#     end
# end

@fastmath function integrate(f::F, X, W) where {F}
    return @inbounds begin
        if size(X, 2) == 2
            s = zero(W[1] * f(X[1, 1], X[1, 2]))
            for i in 1:length(W)
                s += W[i] * f(X[i, 1], X[i, 2])
            end
            s
        elseif size(X, 2) == 3
            s = zero(W[1] * f(X[1, 1], X[1, 2], X[1, 3]))
            for i in 1:length(W)
                s += W[i] * f(X[i, 1], X[i, 2], X[i, 3])
            end
            s
        else
            @error "not implemented"
        end
    end
end

function bench_integrate(outfile, D, N, f, nflop_f)
    X, W = simplexquad(N, D)
    t = @benchmark integrate($f, $X, $W)
    nflop = (nflop_f + 2) * length(W)
    nflop_sec = nflop / minimum(t).time
    return println(outfile,
                   "D=$D N=N   $nflop Flop   $(round(nflop_sec; digits=1)) GFlop/s")
end

open("figures/benchmarks.txt", "w") do outfile
    f2(x, y) = @fastmath 1 + x + y^2
    nflop_f2 = 3
    bench_integrate(outfile, 2, 4, f2, nflop_f2)
    bench_integrate(outfile, 2, 8, f2, nflop_f2)

    f3(x, y, z) = @fastmath 1 + x + y^2 + z^3
    nflop_f3 = 6
    bench_integrate(outfile, 3, 4, f3, nflop_f3)
    return bench_integrate(outfile, 3, 8, f3, nflop_f3)
end


text:
using Test, WAVI, MAT, JLD2, NetCDF, Dates

function output_test(; dt  = 0.5,
                    end_time = 100.,
                    output_freq = 5.,
                    output_format = "jld2",
                    zip_format = "none",
                    prefix = "outfile",
                    dump_vel = false,
                    output_path = ".",
                    pchkpt_freq = Inf,
                    chkpt_freq = Inf,
                    niter0 = 0,
                    output_start = false)

    grid = Grid() #default grid with nx = 80, ny = 10
    bed = WAVI.mismip_plus_bed
    solver_params = SolverParams(maxiter_picard = 1)
    model = Model(grid = grid, bed_elevation = bed, solver_params = solver_params)
    timestepping_params = TimesteppingParams(niter0 = niter0, dt = dt, end_time = end_time, chkpt_freq = chkpt_freq, pchkpt_freq = pchkpt_freq)

    outputs = (h = model.fields.gh.h,
                u = model.fields.gh.u,
                v = model.fields.gh.v,
                b = model.fields.gh.b) #output velocities and thickness

    output_freq = output_freq
    output_params = OutputParams(outputs = outputs,
                        output_freq = output_freq,
                        output_format = output_format,
                        output_path = output_path,
                        zip_format = zip_format,
                        prefix = prefix,
                        dump_vel = dump_vel,
                        output_start = output_start)

    simulation = Simulation(model = model,
                        timestepping_params = timestepping_params,
                        output_params = output_params)

    #perform the simulation
    run_simulation!(simulation)
    return simulation
end

#########################################################
############### Test flags ##############################
test_basic_outputting = true
test_zipping_output   = true
test_output_after_restart = true
test_output_errors = true
##########################################################

@testset "Outputting" begin
    if test_basic_outputting
        @testset "Output files" begin
            @info "Testing outputting..."

            for output_format in ["mat", "jld2"] #check that both mat and jld2 work
                for folder = ["outputs/", "outputs", "./outputs"] #check that both with and without / works
                    isdir(folder) && rm(folder, force = true, recursive = true) #make a new folder
                    mkdir(folder)

                    #run the simulation
                    sim = output_test(output_path = folder,
                                end_time = 100.,
                                output_freq = 5.,
                                prefix = "testoutfile",
                                output_format = output_format)

                    foldersim = sim.output_params.output_path
                    @test foldersim[end] == '/' #test that we do have the / at end of path
                    files = [string(foldersim,f) for f in  readdir(foldersim) if endswith(f, output_format)]
                    println(files)
                    @test length(files) == 20    #check there are the correct number of output files

                    #check they have the correct suffix
                    prefices = [split(str,".")[1] for str in files]
                    suffices = [split(str,".")[end] for str in files]
                    @test all(suffices .== output_format)

                    #check that variables included and correct size
                    fname = string(foldersim,readdir(folder)[1])
                    if output_format == "mat"
                    dict = matread(fname)
                    elseif output_format == "jld2"
                    dict = load(fname)
                    end
                    @test length(dict["t"]) == 1
                    @test size(dict["x"]) == (80,10)
                    @test size(dict["y"]) == (80,10)
                    @test size(dict["b"]) == (80,10)
                    @test size(dict["u"]) == (80,10)
                    @test size(dict["h"]) == (80,10)
                    @test size(dict["v"]) == (80,10)

                    #check it has the correct prefix and number of zeros
                    fname = readdir(folder)[1]
                    @test fname[1:11]  == "testoutfile"

                    #delete the folder
                    rm(folder, force = true, recursive = true)
                end
            end #end loop over output format

        end
        @testset "Test output start" begin
            for output_start in [true, false]
                folder = "outputs/"
                isdir(folder) && rm(folder, force = true, recursive = true) #make a new folder
                mkdir(folder)
                #test that if we only get an output at the start (e.g. if the output freq > end time), then outputted solution matches the expected thing
                sim = output_test(output_path = folder,
                                end_time = 1.,
                                output_freq = 5.,
                                prefix = "testoutfile",
                                output_format = "jld2",
                                output_start = output_start)

                foldersim = sim.output_params.output_path
                if output_start
                    #check only a single file outputted
                    @test length(readdir(foldersim)) == 1

                    fname = string(foldersim,readdir(folder)[1])
                    dict = load(fname)

                    #check that thickness same as IC and time = 0
                    @test dict["h"] == sim.model.initial_conditions.initial_thickness
                    @test dict["t"] == 0.0

                else
                    #check that no file output
                    @test isempty(readdir(foldersim))
                end

                #delete the folder
                rm(folder, force = true, recursive = true)
            end
        end
    end

    if test_zipping_output
        @testset "Zipping output" begin
            @info "Testing zipping output..."
            for output_format in ["mat", "jld2"]
                for output_start in [true, false]
                    folder = "outputs/"
                    isdir(folder) && rm(folder, force = true, recursive = true) #make a new folder
                    mkdir(folder)

                    sim = output_test(zip_format = "nc",
                                        output_path = folder,
                                        output_format = "mat",
                                        output_start = output_start)

                    @test sim isa Simulation
                    fname = string(folder, sim.output_params.prefix, ".nc")
                    @test isfile(fname) #check the zipped file exists

                    #test variables read from nc file
                    x = ncread(fname, "x")
                    y = ncread(fname, "y")
                    h = ncread(fname, "h")
                    u = ncread(fname, "u")
                    v = ncread(fname, "v")
                    b = ncread(fname, "b")
                    t = ncread(fname, "TIME")
                    if output_start
                        @test ncread(fname, "TIME") == 0.:5.:100.
                    else
                        @test ncread(fname, "TIME") == 5.:5.:100.
                    end

                    @test size(x) == (80,)
                    @test size(y) == (10,)
                    @test size(b) == (length(x),length(y),length(t))
                    @test size(u) == (length(x),length(y),length(t))
                    @test size(v) == (length(x),length(y),length(t))
                    @test size(h) == (length(x),length(y),length(t))

                    #delete the folder
                    rm(folder, force = true, recursive = true)
                end
            end
        end
    end

    if test_output_after_restart
        @testset "Outputting after a restart" begin
            @info "Testing outputting after a restart"

            for output_format in ["mat", "jld2"]

                folder = "outputs/"
                isdir(folder) && rm(folder, force = true, recursive = true) #make a new folder
                mkdir(folder)

                #run simulation first time with no zipping
                output_test(; dt  = 0.5,
                            end_time = 10.,
                            output_freq = 0.5,
                            output_format = output_format,
                            output_path = folder,
                            zip_format = "nc",
                            dump_vel = true,
                            pchkpt_freq = 1.)

                #get the time that the first file was outputted (test for https://github.com/RJArthern/WAVI.jl/issues/35)
                first_file_name = joinpath("outputs",string("outfile0000000001.", output_format));
                @test isfile(first_file_name) #check the first output point exists
                dt1 = Dates.unix2datetime(mtime(first_file_name)) #time of last modification

                #run again with different niter0
                sim = output_test(; dt  = 0.5,
                            end_time = 20.,
                            output_freq = 0.5,
                            niter0 = 20,
                            output_format = output_format,
                            output_path = folder,
                            zip_format = "nc",
                            dump_vel = true,
                            pchkpt_freq = 1.)

                #check that the time that first file modified has not changed -- i.e. simulation has not touched outputs at earlier times than it should have
                dt2 = Dates.unix2datetime(mtime(first_file_name))
                @test dt1 == dt2

                #test variables read from nc file
                fname = string(folder, sim.output_params.prefix, ".nc")
                println(fname)
                x = ncread(fname, "x")
                y = ncread(fname, "y")
                h = ncread(fname, "h")
                u = ncread(fname, "u")
                v = ncread(fname, "v")
                b = ncread(fname, "b")
                t = ncread(fname, "TIME")
                @test ncread(fname, "TIME") == 0.5:0.5:20.
                @test size(x) == (80,)
                @test size(y) == (10,)
                @test size(b) == (length(x),length(y),length(t))
                @test size(u) == (length(x),length(y),length(t))
                @test size(v) == (length(x),length(y),length(t))
                @test size(h) == (length(x),length(y),length(t))

                #check we have the right number of output files
                foldersim = sim.output_params.output_path
                @test foldersim[end] == '/' #test that we do have the / at end of path
                files = [string(foldersim,f) for f in  readdir(foldersim) if endswith(f, output_format)]
                println(files)
                @test length(files) == 40    #check there are the correct number of output files

                #now repeat for just one timestep
                sim = output_test(; dt  = 0.5,
                                end_time = 20.5,
                                output_freq = 0.5,
                                niter0 = 40,
                                output_format = output_format,
                                output_path = folder,
                                zip_format = "nc",
                                dump_vel = true,
                                pchkpt_freq = 1.)

                #test variables read from nc file
                fname = string(folder, sim.output_params.prefix, ".nc")
                println(fname)
                x = ncread(fname, "x")
                y = ncread(fname, "y")
                h = ncread(fname, "h")
                u = ncread(fname, "u")
                v = ncread(fname, "v")
                b = ncread(fname, "b")
                t = ncread(fname, "TIME")
                @test ncread(fname, "TIME") == 0.5:0.5:20.5
                @test size(x) == (80,)
                @test size(y) == (10,)
                @test size(b) == (length(x),length(y),length(t))
                @test size(u) == (length(x),length(y),length(t))
                @test size(v) == (length(x),length(y),length(t))
                @test size(h) == (length(x),length(y),length(t))

                #check we have the right number of files
                foldersim = sim.output_params.output_path
                @test foldersim[end] == '/' #test that we do have the / at end of path
                files = [string(foldersim,f) for f in  readdir(foldersim) if endswith(f, output_format)]
                println(files)
                @test length(files) == 41    #check there are the correct number of output files

                #check we have dumped the velocity
                @test isfile(string(foldersim, "outfile.nc")) #check the zipped file exists

                # delete everything you just made
                rm(folder, force = true, recursive = true)
                foreach(rm, filter(endswith(".mat"), readdir()))
                foreach(rm, filter(endswith(".jld2"), readdir()))
                foreach(rm, filter(endswith(".bin"), readdir()))
            end
        end
    end

    if test_output_errors
        @testset "Outputting errors" begin
            @info "Testing outputting errors"
            #check error for weird output format
            @test_throws ArgumentError output_test(output_format = "incorrect_format")

            #check that non-standard zip_format reverts to none
            sim = output_test(end_time = 1., zip_format = "incorrect_zip_format")
            @test sim.output_params.zip_format == "none"

            @test_throws ArgumentError output_test(output_freq =  -2.)

        end
    end
end


text:
function SimJulia.run( mpSim::MPsim, showInfo::Bool=false;
    saveConfig::Bool=true, seed::Integer=-1, sysEnt::Bool=false,
    nCommits::Int=1 )

    if !verifySimulation!( mpSim )
        error( "Simulation configuration is inconsistent, cannot run." )
        return
    elseif showInfo
        @info "Simulation configuration is consistent, initialising processes."
    end  # if !verifySimulation!( mpSim )

    mpSim.showInfo = showInfo
    seedSimulation!( mpSim, seed, sysEnt )
    DBInterface.execute( mpSim.simDB, "BEGIN TRANSACTION" )

    try
        resetSimulation( mpSim )

        # Initialise the recruitment processes.
        for name in keys( mpSim.recruitmentByName )
            for recruitment in mpSim.recruitmentByName[name]
                @process recruitProcess( mpSim.sim, recruitment, mpSim )
            end  # for recruitment in mpSim.recruitmentByName[name]
        end  # for name in keys( mpSim.recruitmentByName )

        # Initialise the default retirement process.
        @process retireProcess( mpSim.sim, mpSim )

        # Initalise the attrition process.
        @process checkAttritionProcess( mpSim.sim, mpSim )

        # Initialise the transition processes.
        orderTransitions!( mpSim )

        for name in keys( mpSim.transitionsByName )
            for transition in mpSim.transitionsByName[name]
                @process transitionProcess( mpSim.sim, transition, mpSim )
            end  # for transition in mpSim.transitionsByName[name]
        end  # for name in keys( mpSim.transitionsByName )

        mpSim.nCommits = max( nCommits, 1 )

        if mpSim.nCommits > 1
            @process dbCommitProcess( mpSim.sim, mpSim )
        end

        # Execute the simulation.
        run( mpSim.sim )
    catch err
        DBInterface.execute( mpSim.simDB, "COMMIT" )
        rethrow( err )
    end  # try

    DBInterface.execute( mpSim.simDB, "COMMIT" )

    if mpSim.showInfo
        println( "Attrition execution processes took ",
            mpSim.attritionExecTime.value / 1000, " seconds." )
    end  # if mpSim.showInfo

    # Save the configuration if needed.
    if saveConfig
        saveSimulationConfiguration( mpSim )
    end  # if saveConfig

    mpSim.isVirgin = true
    return

end  # run( mpSim, showInfo, saveConfig, seed, sysEnt )

include( joinpath( simPath, "dboperations.jl" ) )
include( joinpath( simPath, "snapshot.jl" ) )
include( joinpath( simPath, "multirunsimulation.jl" ) )
include( joinpath( simPrivPath, "attrition.jl" ) )
include( joinpath( simPrivPath, "attribute.jl" ) )
include( joinpath( simPrivPath, "recruitment.jl" ) )
include( joinpath( simPrivPath, "retirement.jl" ) )
include( joinpath( simPrivPath, "transition.jl" ) )
include( joinpath( simPrivPath, "dboperations.jl" ) )
include( joinpath( simPrivPath, "manpowersimulation.jl" ) )


text:
# latex treats # as a special character, so we have to escape it. See:
# https://github.com/sisl/TikzGraphs.jl/issues/12

# latex_escape(s::String) = replace(s, Dict("#"=>"\\#", "&"=>"\\&"))

function latex_escape(s::String)
    s = replace(s, "#"=>"\\#")
    s = replace(s, "&"=>"\\&")
    return s
end

"Convert a symbol or  into a LaTeX label"
function latex_label(sym)
    sym == :(^) && return "\\textasciicircum"  # TikzGraphs chokes on ^

    return latex_escape(string("\\texttt{", sym, "}"))
end

"""
Return a Tikz representation of a tree object.
The tree object must have fields `g` (the graph) and `labels`.
"""
function tikz_representation(tree::LabelledDiGraph)
    labels = String[latex_label(x) for x in tree.labels]
    return TikzGraphs.plot(tree.g, labels)
end

function Base.show(io::IO, mime::MIME"image/svg+xml", tree::LabelledTree)

    p = tikz_representation(tree)  # TikzPicture object
    show(io, mime, p)

end