text:
struct Dictionary{I, T} <: AbstractDictionary{I, T}
    indices::Indices{I}
    values::Vector{T}

    function Dictionary{I, T}(inds::Indices{I}, values::Vector{T}, ::Nothing) where {I, T}
       @assert length(values) == length(_values(inds))
       return new{I,T}(inds, values)
    end
end

"""
    Dictionary{I,T}(;sizehint = 8)

Construct an empty hash-based dictionary. `I` and `T` default to `Any` if not specified. A
`sizehint` may be specified to set the initial size of the hash table, which may speed up
subsequent `insert!` operations.

# Example

```julia
julia> d = Dictionary{Int, Int}()
0-element Dictionary{Int64,Int64}
```
"""
Dictionary(; sizehint = 8) = Dictionary{Any, Any}(; sizehint = sizehint)
Dictionary{I}(; sizehint = 8) where {I} = Dictionary{I, Any}(; sizehint = sizehint)

function Dictionary{I, T}(; sizehint = 8) where {I, T}
    values = Vector{T}()
    sizehint!(values, sizehint)
    Dictionary{I, T}(Indices{I}(; sizehint = sizehint), values, nothing)
end

"""
    Dictionary(indexable)
    Dictionary{I}(indexable)
    Dictionary{I,T}(indexable)

Construct a hash-based dictionary from an indexable input `indexable`, equivalent to
`Dictionary(keys(indexable), values(indexable))`. The input might not be copied.

Note: to construct a dictionary from `Pair`s use the `dictionary` function. See also the
`index` function.

# Examples

```julia
julia> Dictionary(Dict(:a=>1, :b=>2))
2-element Dictionary{Symbol,Int64}
 :a │ 1
 :b │ 2

julia> Dictionary(3:-1:1)
3-element Dictionary{Int64,Int64}
 1 │ 3
 2 │ 2
 3 │ 1
```
"""
Dictionary(indexable) = Dictionary(keys(indexable), values(indexable))
Dictionary{I}(indexable) where {I} = Dictionary{I}(keys(indexable), values(indexable))
Dictionary{I, T}(indexable) where {I, T} = Dictionary{I, T}(keys(indexable), values(indexable))

Dictionary{I, T}(indexable::Dictionary) where {I, T} = Dictionary{I, T}(convert(Indices{I}, keys(indexable)), convert(Vector{T}, indexable.values), nothing)

"""
    Dictionary(inds, values)
    Dictionary{I}(inds, values)
    Dictionary{I, T}(inds, values)

Construct a hash-based dictionary from two iterable inputs `inds` and `values`. The first
value of `inds` will be the index for the first value of `values`. The input might not be
copied.

Note: the values of `inds` must be distinct. Consider using `dictionary(zip(inds, values))`
if they are not. See also the `index` function.

# Example

julia> Dictionary(["a", "b", "c"], [1, 2, 3])
3-element Dictionary{String,Int64}
 "a" │ 1
 "b" │ 2
 "c" │ 3

julia> Dictionary{String, Float64}(["a", "b", "c"], [1, 2, 3])
3-element Dictionary{String,Float6464}
 "a" │ 1.0
 "b" │ 2.0
 "c" │ 3.0
"""
function Dictionary(inds, values)
    return Dictionary(Indices(inds), values)
end

function Dictionary(inds::Indices{I}, values) where {I}
    return Dictionary{I}(inds, values)
end

function Dictionary{I}(inds, values) where {I}
    return Dictionary{I}(Indices{I}(inds), values)
end

function Dictionary{I}(inds::Indices{I}, values) where {I}
    if Base.IteratorEltype(values) === Base.EltypeUnknown()
        values = collect(values)
    end

    return Dictionary{I, eltype(values)}(inds, values)
end

function Dictionary{I, T}(inds, values) where {I, T}
    return Dictionary{I, T}(Indices{I}(inds), values)
end

function Dictionary{I, T}(inds::Indices{I}, values) where {I, T}
    if _holes(inds) != 0
        # TODO instead constructor a vector with holes in it...
        inds = copy(inds)
    end

    iter_size = Base.IteratorSize(values)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        vs = Vector{T}(undef, length(values))
        @inbounds for (i, v) in enumerate(values)
            vs[i] = v
        end
        return Dictionary{I, T}(inds, vs, nothing)
    else
        vs = Vector{T}()
        for v in values
            push!(vs, v)
        end
        return Dictionary{I, T}(inds, vs, nothing)
    end
end

function Dictionary{I, T}(inds::Indices{I}, values::AbstractVector) where {I, T}
    if _holes(inds) != 0
        # TODO instead constructor a vector with holes in it...
        inds = copy(inds)
    end
    return Dictionary{I, T}(inds, convert(Vector{T}, values), nothing)
end

"""
    Dictionary(indices, undef::UndefInitializer)

Construct a `Dictionary` from an iterable of `indices`, where the values are
undefined/uninitialized.

# Example

```julia
julia> Dictionary{Int, Float64}([1,2,3], undef)
3-element Dictionary{Int64,Float64}
 1 │ 6.9220016379355e-310
 2 │ 6.9220016379426e-310
 3 │ 6.92200163794736e-310
```
"""
function Dictionary{I, T}(inds::Indices{I}, ::UndefInitializer) where {I, T}
    vs = Vector{T}(undef, length(_values(inds)))
    return Dictionary{I, T}(inds, vs, nothing)
end

function Base.convert(::Type{Dictionary{I, T}}, dict::Dictionary) where {I, T}
    return Dictionary{I, T}(convert(Indices{I}, dict.indices), convert(Vector{T}, dict.values), nothing)
end
Base.convert(::Type{T}, dict::T) where {T<:Dictionary} = dict

"""
    dictionary(iter)

Construct a new `AbstractDictionary` from an iterable `iter` of key-value `Pair`s (or other
iterables of two elements, such as a two-tuples). The default container type is
`Dictionary`. If duplicate keys are detected, the first encountered value is retained.

See also the `index` function.

# Examples

```julia
julia> dictionary(["a"=>1, "b"=>2, "c"=>3])
3-element Dictionary{String,Int64}
 "a" │ 1
 "b" │ 2
 "c" │ 3

julia> dictionary(["a"=>1, "b"=>2, "c"=>3, "a"=>4])
3-element Dictionary{String,Int64}
 "a" │ 1
 "b" │ 2
 "c" │ 3

julia> dictionary(zip(["a","b","c"], [1,2,3]))
3-element Dictionary{String,Int64}
 "a" │ 1
 "b" │ 2
 "c" │ 3
```
"""
function dictionary(iter)
    return _dictionary(first, last, Dictionary, iter)
end

# An auto-widening Dictionary constructor
function _dictionary(key, value, ::Type{Dictionary}, iter)
    tmp = iterate(iter)
    if tmp === nothing
        IT = Base.@default_eltype(iter)
        I = Core.Compiler.return_type(first, Tuple{IT})
        T = Core.Compiler.return_type(last, Tuple{IT})
        return Dictionary{I, T}()
    end
    (x, s) = tmp
    i = key(x)
    v = value(x)
    dict = Dictionary{typeof(i), typeof(v)}()
    insert!(dict, i, v)
    return __dictionary(key, value, dict, iter, s)
end

# An auto-widening AbstractDictionary constructor
function __dictionary(key, value, dict, iter, s)
    I = keytype(dict)
    T = eltype(dict)
    tmp = iterate(iter, s)
    while tmp !== nothing
        (x, s) = tmp
        i = key(x)
        v = value(x)
        if !(i isa I) && promote_type(typeof(i), I) != I
            new_inds = copy(keys(dict), promote_type(I, typeof(i)))
            new_dict = similar(new_inds, promote_type(T, typeof(v)))
            copyto!(new_dict.values, dict.values)
            (hadtoken, token) = gettoken!(new_dict, i)
            if !hadtoken
                @inbounds settokenvalue!(new_dict, token, v)
            end
            return __dictionary(key, value, new_dict, iter, s)
        elseif !(v isa T) && promote_type(typeof(v), T) != T
            new_dict = copy(dict, promote_type(T, typeof(v)))
            (hadtoken, token) = gettoken!(new_dict, i)
            if !hadtoken
                @inbounds settokenvalue!(new_dict, token, v)
            end
            return __dictionary(key, value, new_dict, iter, s)
        end
        (hadtoken, token) = gettoken!(dict, i)
        if !hadtoken
            @inbounds settokenvalue!(dict, token, v)
        end
        tmp = iterate(iter, s)
    end
    return dict
end

"""
    index(f, iter)

Return a dictionary associating the values `x` of iterable collection `iter` with the key
`f(x)`. If keys are repeated, only the first is kept. Somewhat similar to `unique(f, iter)`

See also the `dictionary` function.

# Examples

```julia
julia> index(first, ["Alice", "Bob", "Charlie"])
3-element Dictionary{Char,String}
 'A' │ "Alice"
 'B' │ "Bob"
 'C' │ "Charlie"

julia> index(iseven, 1:10)
2-element Dictionary{Bool,Int64}
 false │ 1
  true │ 2
```
"""
function index(f, iter)
    _dictionary(f, identity, Dictionary, iter)
end

# indices

Base.keys(dict::Dictionary) = getfield(dict, :indices)
_values(dict::Dictionary) = getfield(dict, :values)

# tokens

tokenized(dict::Dictionary) = _values(dict)

# values

function istokenassigned(dict::Dictionary, (_slot, index))
    return isassigned(_values(dict), index)
end

function istokenassigned(dict::Dictionary, index::Int)
    return isassigned(_values(dict), index)
end

@propagate_inbounds function gettokenvalue(dict::Dictionary, (_slot, index))
    return _values(dict)[index]
end

@propagate_inbounds function gettokenvalue(dict::Dictionary, index::Int)
    return _values(dict)[index]
end

issettable(::Dictionary) = true

@propagate_inbounds function settokenvalue!(dict::Dictionary{<:Any, T}, (_slot, index), value::T) where {T}
    _values(dict)[index] = value
    return dict
end

@propagate_inbounds function settokenvalue!(dict::Dictionary{<:Any, T}, index::Int, value::T) where {T}
    _values(dict)[index] = value
    return dict
end

# insertion

isinsertable(::Dictionary) = true

function gettoken!(dict::Dictionary{I}, i::I) where {I}
    (hadtoken, (slot, index)) = gettoken!(keys(dict), i, (_values(dict),))
    return (hadtoken, (slot, index))
end

function deletetoken!(dict::Dictionary{I, T}, (slot, index)) where {I, T}
    isbitstype(T) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), _values(dict), index-1)
    deletetoken!(keys(dict), (slot, index), (_values(dict),))
    return dict
end

function Base.empty!(dict::Dictionary{I, T}) where {I, T}
    empty!(_values(dict))
    empty!(keys(dict))

    return dict
end

function Base.filter!(pred, dict::Dictionary)
    indices = keys(dict)
    _filter!(token -> pred(@inbounds gettokenvalue(dict, token)), indices, (_values(dict),))
    return dict
end

function Base.filter!(pred, dict::PairDictionary{<:Any, <:Any, <:Dictionary})
    d = dict.d
    indices = keys(d)
    _filter!(token -> pred(@inbounds gettokenvalue(indices, token) => gettokenvalue(d, token)), indices, (_values(d),))
    return dict
end

# Factories

function Base.similar(indices::Indices{I}, ::Type{T}) where {I, T}
    return Dictionary{I, T}(indices, Vector{T}(undef, length(_values(indices))), nothing)
end

"""
    sort!(dict::AbstractDictionary; kwargs...)

Modify `dict` so that it is sorted by its values. The `kwargs` are the usual ordering
options supported by `sort`. Note that this only works on supported types.

See also `sort`, `sortkeys!` and `sortpairs!`.
"""
function Base.sort!(dict::Dictionary; kwargs...)
    inds = keys(dict)
    if inds.holes != 0
        rehash!(inds, length(inds.slots))
    end
    perm = sortperm(dict.values; kwargs...)
    iperm = invperm(perm)
    inds.values = @inbounds inds.values[perm]
    inds.hashes = @inbounds inds.hashes[perm]
    @inbounds for i in keys(inds.slots)
        s = inds.slots[i]
        if s > 0
            inds.slots[i] = iperm[s]
        end
    end
    permute!(dict.values, perm)
    return dict
end

"""
    sortkeys!(dict::AbstractDictionary; kwargs...)

Modify `dict` so that it is sorted by `keys(dict)`. The `kwargs` are the usual ordering
options supported by `sort`. Note that this only works on supported types.

See also `sortkeys`, `sort!` and `sortpairs!`.
"""
function sortkeys!(dict::Dictionary; kwargs...)
    inds = keys(dict)
    if inds.holes != 0
        rehash!(inds, length(inds.slots))
    end
    perm = sortperm(inds.values; kwargs...)
    iperm = invperm(perm)
    inds.values = @inbounds inds.values[perm]
    inds.hashes = @inbounds inds.hashes[perm]
    @inbounds for i in keys(inds.slots)
        s = inds.slots[i]
        if s > 0
            inds.slots[i] = iperm[s]
        end
    end
    permute!(dict.values, perm)
    return dict
end

"""
    sortpairs!(dict::AbstractDictionary; kwargs...)

Modify `dict` so that it is sorted by `pairs(dict)`. The `kwargs` are the usual ordering
options supported by `sort`.

See also `sortpairs`,`sort!` and `sortkeys!`.
"""
function sortpairs!(dict::Dictionary; by = identity, kwargs...)
    inds = keys(dict)
    if inds.holes != 0
        rehash!(inds, length(inds.slots))
    end
    vals = dict.values
    inds_vals = inds.values
    perm = sortperm(keys(dict.values); by = i -> by(@inbounds(inds_vals[i]) => @inbounds(vals[i])), kwargs...)
    iperm = invperm(perm)
    inds.values = @inbounds inds.values[perm]
    inds.hashes = @inbounds inds.hashes[perm]
    @inbounds for i in keys(inds.slots)
        s = inds.slots[i]
        if s > 0
            inds.slots[i] = iperm[s]
        end
    end
    permute!(dict.values, perm)
    return dict
end


text: struct JupiterType <: abstractPlanet end

ext: jl


text:
struct ScalarUniformPeriodicTecno1D{BalanceLaw, T, Fnum, Diss, Reconstruction,
                                    ComputeEntropyVar, EntropyVar, EdgeEntropyVar,
                                    Fluxes, Parallel} <: FiniteVolumeSemidiscretisation
    balance_law::BalanceLaw
    meshx::UniformPeriodicMesh1D{T}
    fnum::Fnum
    diss::Diss
    reconstruction::Reconstruction

    compute_entropy_var::ComputeEntropyVar
    entropy_var::EntropyVar
    edge_entropy_var::EdgeEntropyVar
    fluxes::Fluxes
    parallel::Parallel

    function ScalarUniformPeriodicTecno1D(
            balance_law::BalanceLaw, meshx::UniformPeriodicMesh1D{T},
            fnum::Fnum, diss::Diss, reconstruction::Reconstruction,
            compute_entropy_var::ComputeEntropyVar,
            entropy_var::EntropyVar, edge_entropy_var::EdgeEntropyVar,
            fluxes::Fluxes, parallel::Parallel) where {BalanceLaw, T, Fnum, Diss, Reconstruction, ComputeEntropyVar, EntropyVar, EdgeEntropyVar, Fluxes, Parallel}
        @assert typeof(fluxes) <: AbstractArray{variables(balance_law), 1}
        @assert typeof(entropy_var) <: AbstractArray{variables(balance_law), 1}
        @assert typeof(edge_entropy_var) <: AbstractArray{variables(balance_law), 2}
        @assert size(fluxes, 1) == size(entropy_var, 1) == size(edge_entropy_var, 2) == numcells(meshx)
        @assert size(edge_entropy_var, 1) == 2
        new{BalanceLaw, T, Fnum, Diss, Reconstruction, ComputeEntropyVar, EntropyVar, EdgeEntropyVar, Fluxes, Parallel}(
            balance_law, meshx, fnum, diss, reconstruction, compute_entropy_var, entropy_var, edge_entropy_var, fluxes, parallel)
    end
end

function ScalarUniformPeriodicTecno1D(balance_law, meshx, fnum, diss, reconstruction,
                                      compute_entropy_var, parallel=Val{:serial}())
    fluxes = Array{variables(balance_law)}(undef, numcells(meshx))
    entropy_var = Array{variables(balance_law)}(undef, numcells(meshx))
    edge_entropy_var = Array{variables(balance_law)}(undef, 2, numcells(meshx))
    ScalarUniformPeriodicTecno1D(balance_law, meshx, fnum, diss, reconstruction,
                                 compute_entropy_var, entropy_var, edge_entropy_var, fluxes, parallel)
end

function Base.show(io::IO, fv::ScalarUniformPeriodicTecno1D)
    print(io, "TeCNO finite volume method ",
            "\n  Balance law:    ", fv.balance_law,
            "\n  Mesh:           ", fv.meshx,
            "\n  Numerical flux: ", fv.fnum,
            "\n  Dissipation:    ", fv.diss,
            "\n  Reconstruction: ", fv.reconstruction)
end

function evaluate_coefficients(u, fv::ScalarUniformPeriodicTecno1D,
                               npoints=2*order(fv.reconstruction))
    evaluate_coefficients(u, fv.balance_law, fv.meshx, fv.reconstruction, npoints)
end

@noinline function (semidisc::ScalarUniformPeriodicTecno1D)(du, u, p, t)
    @boundscheck begin
        if size(u) != size(du)
            error("size(u) = $(size(u)) != $(size(du)) = size(du)")
        end
        @assert length(u) == numcells(semidisc.meshx)
        @assert length(u) >= stencil_width(semidisc.reconstruction) - 1
        if eltype(u) != variables(semidisc.balance_law)
            error("eltype(u) == $(eltype(u)) != $(variables(semidisc.balance_law)) == variables(semidisc.balance_law)")
        end
    end

    fill!(du, zero(eltype(du)))
    add_numerical_fluxes!(du, u, semidisc, t)

    nothing
end

function add_numerical_fluxes!(du, u, fv::ScalarUniformPeriodicTecno1D, t)
    @unpack balance_law, meshx, fnum, diss, reconstruction, compute_entropy_var, entropy_var, edge_entropy_var, fluxes, parallel = fv

    # L² entropy
    entropy_var .= compute_entropy_var.(u, balance_law)
    reconstruct!(edge_entropy_var, entropy_var, balance_law, meshx, reconstruction, stencil_width_val(reconstruction), parallel)
    compute_EC_fluxes!(fluxes, u, balance_law, meshx, fnum, stencil_width_val(reconstruction), parallel)
    compute_ES_fluxes!(fluxes, edge_entropy_var, balance_law, meshx, diss, parallel)
    compute_du!(du, fluxes, meshx, parallel)
end

# Entropy conservative fluxes
function compute_EC_fluxes!(fluxes, u, balance_law, meshx::UniformPeriodicMesh1D,
                            fnum, ::Val{1}, parallel)
    @inbounds fluxes[1] = fnum(u[end], u[1], balance_law)
    @inbounds for edge in 2:length(fluxes)
        fluxes[edge] = fnum(u[edge-1], u[edge], balance_law)
    end
    nothing
end

function compute_EC_fluxes!(fluxes, u, balance_law, meshx::UniformPeriodicMesh1D,
                            fnum, ::Val{3}, parallel)
    T = eltype(u)
    c1 = T(4) / 3
    c2 = T(-1) / 6

    @inbounds begin
        fluxes[1] = c1 * fnum(u[end], u[1], balance_law) +
                        c2 * fnum(u[end-1], u[1], balance_law) +
                        c2 * fnum(u[end], u[2], balance_law)
        fluxes[2] = c1 * fnum(u[1], u[2], balance_law) +
                        c2 * fnum(u[end], u[2], balance_law) +
                        c2 * fnum(u[1], u[3], balance_law)
        fluxes[end] = c1 * fnum(u[end-1], u[end], balance_law) +
                        c2 * fnum(u[end-2], u[end], balance_law) +
                        c2 * fnum(u[end-1], u[1], balance_law)
    end
    @inbounds for edge in 3:length(fluxes)-1
        fluxes[edge] = c1 * fnum(u[edge-1], u[edge], balance_law) +
                        c2 * fnum(u[edge-2], u[edge], balance_law) +
                        c2 * fnum(u[edge-1], u[edge+1], balance_law)
    end
    nothing
end

function compute_EC_fluxes!(fluxes, u, balance_law, meshx::UniformPeriodicMesh1D,
                            fnum, ::Val{5}, parallel)
    T = eltype(u)
    c1 = T(3) / 2
    c2 = T(-3) / 10
    c3 = T(1) / 30

    @inbounds begin
        fluxes[1] = c1 * fnum(u[end], u[1], balance_law) +
                        c2 * fnum(u[end-1], u[1], balance_law) +
                        c2 * fnum(u[end], u[2], balance_law) +
                        c3 * fnum(u[end-2], u[1], balance_law) +
                        c3 * fnum(u[end-1], u[2], balance_law) +
                        c3 * fnum(u[end], u[3], balance_law)
        fluxes[2] = c1 * fnum(u[1], u[2], balance_law) +
                        c2 * fnum(u[end], u[2], balance_law) +
                        c2 * fnum(u[1], u[3], balance_law) +
                        c3 * fnum(u[end-1], u[2], balance_law) +
                        c3 * fnum(u[end], u[3], balance_law) +
                        c3 * fnum(u[1], u[4], balance_law)
        fluxes[3] = c1 * fnum(u[2], u[3], balance_law) +
                        c2 * fnum(u[1], u[3], balance_law) +
                        c2 * fnum(u[2], u[4], balance_law) +
                        c3 * fnum(u[end], u[3], balance_law) +
                        c3 * fnum(u[1], u[4], balance_law) +
                        c3 * fnum(u[2], u[5], balance_law)
        fluxes[end-1] = c1 * fnum(u[end-2], u[end-1], balance_law) +
                        c2 * fnum(u[end-3], u[end-1], balance_law) +
                        c2 * fnum(u[end-2], u[end], balance_law) +
                        c3 * fnum(u[end-4], u[end-1], balance_law) +
                        c3 * fnum(u[end-3], u[end], balance_law) +
                        c3 * fnum(u[end-2], u[1], balance_law)
        fluxes[end] = c1 * fnum(u[end-1], u[end], balance_law) +
                        c2 * fnum(u[end-2], u[end], balance_law) +
                        c2 * fnum(u[end-1], u[1], balance_law) +
                        c3 * fnum(u[end-3], u[end], balance_law) +
                        c3 * fnum(u[end-2], u[1], balance_law) +
                        c3 * fnum(u[end-1], u[2], balance_law)
    end
    @inbounds for edge in 4:length(fluxes)-2
        fluxes[edge] = c1 * fnum(u[edge-1], u[edge], balance_law) +
                        c2 * fnum(u[edge-2], u[edge], balance_law) +
                        c2 * fnum(u[edge-1], u[edge+1], balance_law) +
                        c3 * fnum(u[edge-3], u[edge], balance_law) +
                        c3 * fnum(u[edge-2], u[edge+1], balance_law) +
                        c3 * fnum(u[edge-1], u[edge+2], balance_law)
    end
    nothing
end

# Add dissipation operators to get entropy stable fluxes
function compute_ES_fluxes!(fluxes, edge_entropy_var, balance_law, meshx::UniformPeriodicMesh1D,
                            diss, parallel)
    @inbounds fluxes[1] += diss(edge_entropy_var[2,end], edge_entropy_var[1,1], balance_law)
    @inbounds for edge in 2:length(fluxes)
        fluxes[edge] += diss(edge_entropy_var[2,edge-1], edge_entropy_var[1,edge], balance_law)
    end
    nothing
end

function compute_ES_fluxes!(fluxes, edge_entropy_var, balance_law, meshx::UniformPeriodicMesh1D,
                            diss, parallel::Val{:threads})
    @inbounds fluxes[1] += diss(edge_entropy_var[2,end], edge_entropy_var[1,1], balance_law)
    @inbounds Threads.@threads for edge in 2:length(fluxes)
        fluxes[edge] += diss(edge_entropy_var[2,edge-1], edge_entropy_var[1,edge], balance_law)
    end
    nothing
end

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

function semidiscretise(fv::ScalarUniformPeriodicTecno1D, u₀func, tspan)
    u₀ = compute_coefficients(u₀func, fv.meshx, order(fv.reconstruction))

    ODEProblem(fv, u₀, tspan)
end

"""
    max_dt(t, u, fv::ScalarUniformPeriodicTecno1D, cfl=0.5)

Compute the maximal time step `dt` satisfying the CFL condition `dt <= cfl * dx / 2λ`,
where `dx` is the length of a cell and `λ` the greatest absolute value of the
speed in this cell.
"""
function max_dt(t, u, fv::ScalarUniformPeriodicTecno1D, cfl=0.5)
    @unpack balance_law, meshx = fv
    dt = mapreduce(cell->volume(cell,meshx) / max_abs_speed(u[cell],balance_law),
                    min, cell_indices(meshx), init=typemax(t))
    if dt == Inf
        dt = mapreduce(cell->volume(cell,meshx), min, cell_indices(meshx), init=dt)
    end
    cfl * dt
end

function integrate(func, u, fv::ScalarUniformPeriodicTecno1D)
    integrate(func, u, fv.meshx)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct AHMCSampleID <: SampleID
    chainid::Int32
    chaincycle::Int32
    stepno::Int64
    sampletype::Int32
    hamiltonian_energy::Float64
    tree_depth::Int32
    divergent::Bool
    step_size::Float64
end

const AHMCSampleIDVector{TV<:AbstractVector{<:Int32},UV<:AbstractVector{<:Int64},
                         FV<:AbstractVector{<:Float64},BV<:AbstractVector{<:Bool}} = StructArray{
    AHMCSampleID,
    1,
    NamedTuple{
        (:chainid, :chaincycle, :stepno, :sampletype, :hamiltonian_energy, :tree_depth, :divergent, :step_size),
        Tuple{TV,TV,UV,TV,FV,TV,BV,FV}
    },
    Int
}

AHMCSampleIDVector(contents::NTuple{8,Any}) = StructArray{AHMCSampleID}(contents)

AHMCSampleIDVector(::UndefInitializer, len::Integer) = AHMCSampleIDVector((
    Vector{Int32}(undef, len), Vector{Int32}(undef, len),
    Vector{Int64}(undef, len), Vector{Int32}(undef, len),
    Vector{Float64}(undef, len), Vector{Int32}(undef, len),
    Vector{Bool}(undef, len), Vector{Float64}(undef, len),
))

AHMCSampleIDVector() = AHMCSampleIDVector(undef, 0)

_create_undef_vector(::Type{AHMCSampleID}, len::Integer) = AHMCSampleIDVector(undef, len)

# Specialize comparison, currently StructArray seems fall back to `(==)(A::AbstractArray, B::AbstractArray)`
import Base.==
function(==)(A::AHMCSampleIDVector, B::AHMCSampleIDVector)
    A.chainid == B.chainid &&
    A.chaincycle == B.chaincycle &&
    A.stepno == B.stepno &&
    A.sampletype == B.sampletype &&
    A.hamiltonian_energy == B.hamiltonian_energy &&
    A.tree_depth == B.tree_depth &&
    A.divergent == B.divergent &&
    A.step_size == B.step_size
end

function Base.merge!(X::AHMCSampleIDVector, Xs::AHMCSampleIDVector...)
    for Y in Xs
        append!(X, Y)
    end
    X
end

Base.merge(X::AHMCSampleIDVector, Xs::AHMCSampleIDVector...) = merge!(deepcopy(X), Xs...)


text:
module F16Stevens

using StaticArrays
using FlightMechanicsUtils
using ..FlightMechanicsSimulator

import ..FlightMechanicsSimulator: get_mass, get_inertia_tensor, get_xcg_mac
import ..FlightMechanicsSimulator: get_chord, get_surface, get_wing_span
import ..FlightMechanicsSimulator: calculate_prop_forces_moments
import ..FlightMechanicsSimulator: calculate_prop_gyro_effects
import ..FlightMechanicsSimulator: calculate_pdot
import ..FlightMechanicsSimulator: calculate_aero_forces_moments
import ..FlightMechanicsSimulator: tgear

export F16
include("aircraft.jl")
include("adc.jl")
include("engine.jl")
include("aero.jl")

end


text:
function sde_determine_initdt(u0::uType,t::tType,tdir,dtmax,abstol,reltol,internalnorm,prob,order,integrator) where {tType,uType}

  if integrator.P !== nothing
    # Don't have a good estimate with jumps
    return tdir*dtmax/1e6
  end

  f = prob.f
  g = prob.g
  p = prob.p
  d₀ = internalnorm(ArrayInterface.aos_to_soa(u0)./(abstol.+internalnorm.(u0,t).*reltol),t)
  dtmin = nextfloat(integrator.opts.dtmin)
  smalldt = tType(1//10^(6))
  if !isinplace(prob)
    f₀ = f(u0,p,t)
    if integrator.opts.verbose && any(x->any(isnan,x),f₀)
      @warn("First function call for f produced NaNs. Exiting.")
    end
    g₀ = 3g(u0,p,t)
    if integrator.opts.verbose && any(x->any(isnan,x),g₀)
      @warn("First function call for g produced NaNs. Exiting.")
    end
  else
    f₀ = zero(u0)
    if prob.noise_rate_prototype !== nothing
      g₀ = zero(prob.noise_rate_prototype)
    else
      g₀ = zero(u0)
    end
    f(f₀,u0,p,t)
    if integrator.opts.verbose && any(x->any(isnan,x),f₀)
      @warn("First function call for f produced NaNs. Exiting.")
    end
    g(g₀,u0,p,t); g₀.*=3
    if integrator.opts.verbose && any(x->any(isnan,x),g₀)
      @warn("First function call for g produced NaNs. Exiting.")
    end
  end

  d₁ = internalnorm(max.(internalnorm.(f₀.+g₀,t),internalnorm.(f₀.-g₀,t))./(abstol.+internalnorm.(u0,t).*reltol),t)
  T0 = typeof(d₀)
  T1 = typeof(d₁)
  if d₀ < T0(1//10^(5)) || d₁ < T1(1//10^(5))
    dt₀ = smalldt
  else
    dt₀ = tType(0.01*(d₀/d₁))
  end
  dt₀ = min(dt₀,tdir*dtmax)
  u₁ = u0 .+ tdir.*dt₀.*f₀
  if !isinplace(prob)
    f₁ = f(u₁,p,t+tdir*dt₀)
    g₁ = 3g(u₁,p,t+tdir*dt₀)
  else
    f₁ = zero(u0)
    if prob.noise_rate_prototype !== nothing
      g₁ = zero(prob.noise_rate_prototype)
    else
      g₁ = zero(u0)
    end
    f(f₁,u0,p,t)
    g(g₁,u0,p,t); g₁.*=3
  end
  ΔgMax = max.(internalnorm.(g₀.-g₁,t),internalnorm.(g₀.+g₁,t))
  d₂ = internalnorm(max.(internalnorm.(f₁.-f₀.+ΔgMax,t),internalnorm.(f₁.-f₀.-ΔgMax,t))./(abstol.+internalnorm.(u0,t).*reltol),t)./dt₀
  if max(d₁,d₂)<=T1(1//Int64(10)^(15))
    dt₁ = max(smalldt,dt₀*1//10^(3))
  else
    dt₁ = tType(10.0^(-(2+log10(max(d₁,d₂)))/(order+.5)))
  end
  dt = tdir*max(dtmin, min(100dt₀,dt₁,tdir*dtmax))
end


text:
Lar = LinearAlgebraicRepresentation

"""
    bbox(vertices::Points)

The axis aligned bounding box of the provided set of n-dim `vertices`.

The box is returned as the couple of `Points` of the two opposite corners of the box.
"""
function bbox(vertices::Points)
    minimum = mapslices(x->min(x...), vertices, dims=1)
    maximum = mapslices(x->max(x...), vertices, dims=1)
    minimum, maximum
end

"""
    bbox_contains(container, contained)

Check if the axis aligned bounding box `container` contains `contained`.

Each input box must be passed as the couple of `Points` standing on the opposite corners of the box.
"""
function bbox_contains(container, contained)
    b1_min, b1_max = container
    b2_min, b2_max = contained
    all(map((i,j,k,l)->i<=j<=k<=l, b1_min, b2_min, b2_max, b1_max))
end

"""
    face_area(V::Points, EV::Cells, face::Cell)

The area of `face` given a geometry `V` and an edge topology `EV`.
"""
function face_area(V::Points, EV::Cells, face::Cell)
    return face_area(V, build_copEV(EV), face)
end

function face_area(V::Points, EV::ChainOp, face::Cell)
    function triangle_area(triangle_points::Points)
        ret = ones(3,3)
        ret[:, 1:2] = triangle_points
        return .5*det(ret)
    end

    area = 0

    fv = buildFV(EV, face)

    verts_num = length(fv)
    v1 = fv[1]

    for i in 2:(verts_num-1)

        v2 = fv[i]
        v3 = fv[i+1]

        area += triangle_area(V[[v1, v2, v3], :])
    end

    return area
end

"""
    skel_merge(V1::Points, EV1::ChainOp, V2::Points, EV2::ChainOp)

Merge two **1-skeletons**
"""
function skel_merge(V1::Points, EV1::ChainOp, V2::Points, EV2::ChainOp)
    V = [V1; V2]
    EV = blockdiag(EV1,EV2)
    return V, EV
end

"""
    skel_merge(V1::Points, EV1::ChainOp, FE1::ChainOp, V2::Points, EV2::ChainOp, FE2::ChainOp)

Merge two **2-skeletons**
"""
function skel_merge(V1::Points, EV1::ChainOp, FE1::ChainOp, V2::Points, EV2::ChainOp, FE2::ChainOp)
    FE = blockdiag(FE1,FE2)
    V, EV = skel_merge(V1, EV1, V2, EV2)
    return V, EV, FE
end

"""
    delete_edges(todel, V::Points, EV::ChainOp)

Delete edges and remove unused vertices from a **2-skeleton**.

Loop over the `todel` edge index list and remove the marked edges from `EV`.
The vertices in `V` which remained unconnected after the edge deletion are deleted too.
"""
function delete_edges(todel, V::Points, EV::ChainOp)
    tokeep = setdiff(collect(1:EV.m), todel)
    EV = EV[tokeep, :]

    vertinds = 1:EV.n
    todel = Array{Int64, 1}()
    for i in vertinds
        if length(EV[:, i].nzind) == 0
            push!(todel, i)
        end
    end

    tokeep = setdiff(vertinds, todel)
    EV = EV[:, tokeep]
    V = V[tokeep, :]

    return V, EV
end

"""
    buildFV(EV::Cells, face::Cell)

The list of vertex indices that expresses the given `face`.

The returned list is made of the vertex indices ordered following the traversal order to keep a coherent face orientation.
The edges are need to understand the topology of the face.

In this method the input face must be expressed as a `Cell`(=`SparseVector{Int8, Int}`) and the edges as `Cells`.
"""
function buildFV(EV::Cells, face::Cell)
    return buildFV(build_copEV(EV), face)
end

"""
    buildFV(copEV::ChainOp, face::Cell)

The list of vertex indices that expresses the given `face`.

The returned list is made of the vertex indices ordered following the traversal order to keep a coherent face orientation.
The edges are need to understand the topology of the face.

In this method the input face must be expressed as a `Cell`(=`SparseVector{Int8, Int}`) and the edges as `ChainOp`.
"""
function buildFV(copEV::ChainOp, face::Cell)
    startv = -1
    nextv = 0
    edge = 0

    vs = Array{Int64, 1}()

    while startv != nextv
        if startv < 0
            edge = face.nzind[1]
            startv = copEV[edge,:].nzind[face[edge] < 0 ? 2 : 1]
            push!(vs, startv)
        else
            edge = setdiff(intersect(face.nzind, copEV[:, nextv].nzind), edge)[1]
        end
        nextv = copEV[edge,:].nzind[face[edge] < 0 ? 1 : 2]
        push!(vs, nextv)

    end

    return vs[1:end-1]
end

"""
    buildFV(copEV::ChainOp, face::Array{Int, 1})

The list of vertex indices that expresses the given `face`.

The returned list is made of the vertex indices ordered following the traversal order to keep a coherent face orientation.
The edges are need to understand the topology of the face.

In this method the input face must be expressed as a list of vertex indices and the edges as `ChainOp`.
"""
function buildFV(copEV::ChainOp, face::Array{Int, 1})
    startv = face[1]
    nextv = startv

    vs = []
    visited_edges = []

    while true
        curv = nextv
        push!(vs, curv)

        edge = 0

        for edgeEx in copEV[:, curv].nzind
            nextv = setdiff(copEV[edgeEx, :].nzind, curv)[1]
            if nextv in face && (nextv == startv || !(nextv in vs)) && !(edgeEx in visited_edges)
                edge = edgeEx
                break
            end
        end

        push!(visited_edges, edge)

        if nextv == startv
            break
        end
    end

    return vs
end

"""
    build_copFE(FV::Cells, EV::Cells)

The signed `ChainOp` from 1-cells (edges) to 2-cells (faces)
"""
function build_copFE(FV::Cells, EV::Cells)
    faces = []

    for face in FV
        f = []
        for (i,v) in enumerate(face)
            edge = [v, face[(i==length(face)) ? 1 : i+1]]
            ord_edge = sort(edge)

            edge_idx = findfirst(e->e==ord_edge, EV)

            push!(f, (edge_idx, sign(edge[2]-edge[1])))
        end

        push!(faces, f)
    end

    FE = spzeros(Int8, length(faces), length(EV))

    for (i,f) in enumerate(faces)
        for e in f
            FE[i, e[1]] = e[2]
        end
    end

    return FE
end

"""
    build_copEV(EV::Cells, signed=true)

The signed (or not) `ChainOp` from 0-cells (vertices) to 1-cells (edges)
"""
function build_copEV(EV::Cells, signed=true)
    setValue = [-1, 1]
    if signed == false
        setValue = [1, 1]
    end

    maxv = max(map(x->max(x...), EV)...)
    copEV = spzeros(Int8, length(EV), maxv)

    for (i,e) in enumerate(EV)
        e = sort(collect(e))
        copEV[i, e] = setValue
    end

    return copEV
end

"""
    build_cops(edges::Cells, faces::Cells)

The vertices-edges and edges-faces chain operators (`copEV::ChainOp`, `copFE::ChainOp`)
"""
function build_cops(edges::Cells, faces::Cells)
    copEV = build_copEV(edges)
    FV = Cells(map(x->buildFV(copEV,x), faces))
    copFE = build_copFE(FV, edges)

    return [copEV, copFE]
end

"""
    vin(vertex, vertices_set)

Checks if `vertex` is one of the vertices inside `vertices_set`
"""
function vin(vertex, vertices_set)
    for v in vertices_set
        if vequals(vertex, v)
            return true
        end
    end
    return false
end

"""
    vequals(v1, v2)

Check the equality between vertex `v1` and vertex `v2`
"""
function vequals(v1, v2)
    err = 10e-8
    return length(v1) == length(v2) && all(map((x1, x2)->-err < x1-x2 < err, v1, v2))
end

"""
    triangulate(model::LARmodel)

Full constrained Delaunnay triangulation of the given 3-dimensional `LARmodel`
"""
function triangulate(model::LARmodel)
    V, topology = model
    cc = build_cops(topology...)
    return triangulate(V, cc)
end

"""
    triangulate(V::Points, cc::ChainComplex)

Full constrained Delaunnay triangulation of the given 3-dimensional model (given with topology as a `ChainComplex`)
"""
function triangulate(V::Points, cc::ChainComplex)
    copEV, copFE = cc

    triangulated_faces = Array{Any, 1}(undef, copFE.m)

	function vcycle( copEV::Lar.ChainOp, copFE::Lar.ChainOp, f::Int64 )
		edges,signs = findnz(copFE[f,:])
		vpairs = [s>0 ? findnz(copEV[e,:])[1] : reverse(findnz(copEV[e,:])[1])
					for (e,s) in zip(edges,signs)]
		vdict = Dict((v1,v2) for (v1,v2) in  vpairs)

		v0 = collect(vdict)[1][1]
		chain_0 = Int64[v0]
		v = vdict[v0]
		while v ≠ v0
			push!(chain_0,v)
			v = vdict[v]
		end
		return chain_0
	end

    for f in 1:copFE.m
        if f % 10 == 0
            print(".")
        end

        edges_idxs = copFE[f, :].nzind
        edge_num = length(edges_idxs)
        edges = zeros(Int64, edge_num, 2)

        #fv = Lar.buildFV(copEV, copFE[f, :])
        fv = vcycle(copEV, copFE, f)

        vs = V[fv, :]

        v1 = normalize(vs[2, :] - vs[1, :])
        v2 = [0 0 0]
        v3 = [0 0 0]
        err = 1e-8
        i = 3
        while -err < norm(v3) < err
            v2 = normalize(vs[i, :] - vs[1, :])
            v3 = cross(v1, v2)
            i = i + 1
        end
        M = reshape([v1; v2; v3], 3, 3)

        vs = (vs*M)[:, 1:2]

        for i in 1:length(fv)
            edges[i, 1] = fv[i]
            edges[i, 2] = i == length(fv) ? fv[1] : fv[i+1]
        end

        triangulated_faces[f] =
        	Triangle.constrained_triangulation(vs, fv, edges, fill(true, edge_num))

        tV = (V*M)[:, 1:2]

        area = face_area(tV, copEV, copFE[f, :])
        if area < 0
            for i in 1:length(triangulated_faces[f])
                triangulated_faces[f][i] = triangulated_faces[f][i][end:-1:1]
            end
        end
    end

    return triangulated_faces
end

"""
    point_in_face(point, V::Points, copEV::ChainOp)

Check if `point` is inside the area of the face bounded by the edges in `copEV`
"""
function point_in_face(point, V::Points, copEV::ChainOp)

    function pointInPolygonClassification(V,EV)

        function crossingTest(new, old, status, count)
        if status == 0
            status = new
            return status, (count + 0.5)
        else
            if status == old
                return 0, (count + 0.5)
            else
                return 0, (count - 0.5)
            end
        end
        end

        function setTile(box)
        tiles = [[9,1,5],[8,0,4],[10,2,6]]
        b1,b2,b3,b4 = box
        function tileCode(point)
            x,y = point
            code = 0
            if y>b1 code=code|1 end
            if y<b2 code=code|2 end
            if x>b3 code=code|4 end
            if x<b4 code=code|8 end
            return code
        end
        return tileCode
        end

        function pointInPolygonClassification0(pnt)
            x,y = pnt
            xmin,xmax,ymin,ymax = x,x,y,y
            tilecode = setTile([ymax,ymin,xmax,xmin])
            count,status = 0,0

            for k in 1:EV.m
                edge = EV[k,:]
                p1, p2 = V[edge.nzind[1], :], V[edge.nzind[2], :]
                (x1,y1),(x2,y2) = p1,p2
                c1,c2 = tilecode(p1),tilecode(p2)
                c_edge, c_un, c_int = xor(c1, c2), c1|c2, c1&c2

                if (c_edge == 0) & (c_un == 0) return "p_on"
                elseif (c_edge == 12) & (c_un == c_edge) return "p_on"
                elseif c_edge == 3
                    if c_int == 0 return "p_on"
                    elseif c_int == 4 count += 1 end
                elseif c_edge == 15
                    x_int = ((y-y2)*(x1-x2)/(y1-y2))+x2
                    if x_int > x count += 1
                    elseif x_int == x return "p_on" end
                elseif (c_edge == 13) & ((c1==4) | (c2==4))
                        status, count = crossingTest(1,2,status,count)
                elseif (c_edge == 14) & ((c1==4) | (c2==4))
                        status, count = crossingTest(2,1,status,count)
                elseif c_edge == 7 count += 1
                elseif c_edge == 11 count = count
                elseif c_edge == 1
                    if c_int == 0 return "p_on"
                    elseif c_int == 4
                        status, count = crossingTest(1,2,status,count)
                    end
                elseif c_edge == 2
                    if c_int == 0 return "p_on"
                    elseif c_int == 4
                        status, count = crossingTest(2,1,status,count)
                    end
                elseif (c_edge == 4) & (c_un == c_edge) return "p_on"
                elseif (c_edge == 8) & (c_un == c_edge) return "p_on"
                elseif c_edge == 5
                    if (c1==0) | (c2==0) return "p_on"
                    else
                        status, count = crossingTest(1,2,status,count)
                    end
                elseif c_edge == 6
                    if (c1==0) | (c2==0) return "p_on"
                    else
                        status, count = crossingTest(2,1,status,count)
                    end
                elseif (c_edge == 9) & ((c1==0) | (c2==0)) return "p_on"
                elseif (c_edge == 10) & ((c1==0) | (c2==0)) return "p_on"
                end
            end

            if (round(count)%2)==1
                return "p_in"
            else
                return "p_out"
            end
        end
        return pointInPolygonClassification0
    end

    return pointInPolygonClassification(V, copEV)(point) == "p_in"
end

"""
    lar2obj(V::Points, cc::ChainComplex)

Triangulated OBJ string representation of the model passed as input.

Use this function to export LAR models into OBJ

# Example

```julia
	julia> cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
	[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
	[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )

	julia> cube_2 = Lar.Struct([Lar.t(0,0,0.5), Lar.r(0,0,pi/3), cube_1])

	julia> V, FV, EV = Lar.struct2lar(Lar.Struct([ cube_1, cube_2 ]))

	julia> V, bases, coboundaries = Lar.chaincomplex(V,FV,EV)

	julia> (EV, FV, CV), (copEV, copFE, copCF) = bases, coboundaries

	julia> FV # bases[2]
	18-element Array{Array{Int64,1},1}:
	 [1, 3, 4, 6]
	 [2, 3, 5, 6]
	 [7, 8, 9, 10]
	 [1, 2, 3, 7, 8]
	 [4, 6, 9, 10, 11, 12]
	 [5, 6, 11, 12]
	 [1, 4, 7, 9]
	 [2, 5, 11, 13]
	 [2, 8, 10, 11, 13]
	 [2, 3, 14, 15, 16]
	 [11, 12, 13, 17]
	 [11, 12, 13, 18, 19, 20]
	 [2, 3, 13, 17]
	 [2, 13, 14, 18]
	 [15, 16, 19, 20]
	 [3, 6, 12, 15, 19]
	 [3, 6, 12, 17]
	 [14, 16, 18, 20]

	julia> CV # bases[3]
	3-element Array{Array{Int64,1},1}:
	 [2, 3, 5, 6, 11, 12, 13, 14, 15, 16, 18, 19, 20]
	 [2, 3, 5, 6, 11, 12, 13, 17]
	 [1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 17]

	julia> copEV # coboundaries[1]
	34×20 SparseMatrixCSC{Int8,Int64} with 68 stored entries: ...

	julia> copFE # coboundaries[2]
	18×34 SparseMatrixCSC{Int8,Int64} with 80 stored entries: ...

	julia> copCF # coboundaries[3]
	4×18 SparseMatrixCSC{Int8,Int64} with 36 stored entries: ...

	objs = Lar.lar2obj(V'::Lar.Points, [coboundaries...])

	open("./two_cubes.obj", "w") do f
    	write(f, objs)
	end

```
"""
function lar2obj(V::Points, cc::ChainComplex)
    copEV, copFE, copCF = cc

    obj = ""
    for v in 1:size(V, 1)
        obj = string(obj, "v ",
        	round(V[v, 1], digits=6), " ",
        	round(V[v, 2], digits=6), " ",
        	round(V[v, 3], digits=6), "\n")
    end

    print("Triangulating")
    triangulated_faces = triangulate(V, cc[1:2])
    println("DONE")

    for c in 1:copCF.m
        obj = string(obj, "\ng cell", c, "\n")
        for f in copCF[c, :].nzind
            triangles = triangulated_faces[f]
            for tri in triangles
                t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
                obj = string(obj, "f ", t[1], " ", t[2], " ", t[3], "\n")
            end
        end
    end

    return obj
end

"""
    obj2lar(path)

Read OBJ file at `path` and create a 2-skeleton as `Tuple{Points, ChainComplex}` from it.

This function does not care about eventual internal grouping inside the OBJ file.
"""
function obj2lar(path)
    vs = Array{Float64, 2}(undef, 0, 3)
    edges = Array{Array{Int, 1}, 1}()
    faces = Array{Array{Int, 1}, 1}()

    open(path, "r") do fd
        for line in eachline(fd)
            elems = split(line)
            if length(elems) > 0
                if elems[1] == "v"

                    x = parse(Float64, elems[2])
                    y = parse(Float64, elems[3])
                    z = parse(Float64, elems[4])
                    vs = [vs; x y z]

                elseif elems[1] == "f"
                    # Ignore the vertex tangents and normals
                    v1 = parse(Int, split(elems[2], "/")[1])
                    v2 = parse(Int, split(elems[3], "/")[1])
                    v3 = parse(Int, split(elems[4], "/")[1])

                    e1 = sort([v1, v2])
                    e2 = sort([v2, v3])
                    e3 = sort([v1, v3])

                    if !(e1 in edges)
                        push!(edges, e1)
                    end
                    if !(e2 in edges)
                        push!(edges, e2)
                    end
                    if !(e3 in edges)
                        push!(edges, e3)
                    end

                    push!(faces, sort([v1, v2, v3]))
                end
            end
        end
    end

    return vs, build_cops(edges, faces)
end

"""
    binaryRange(n)

Generate the first `n` binary numbers in string padded for max `2^n` length
"""
function binaryRange(n)
    return string.(range(0, length=2^n), base=2, pad=n)
end

function space_arrangement(V::Points, EV::ChainOp, FE::ChainOp, multiproc::Bool=false)

    fs_num = size(FE, 1)
    sp_idx = Lar.Arrangement.spatial_index(V, EV, FE)

    global rV = Lar.Points(undef, 0,3)
    global rEV = SparseArrays.spzeros(Int8,0,0)
    global rFE = SparseArrays.spzeros(Int8,0,0)

    if (multiproc == true)
        in_chan = Distributed.RemoteChannel(()->Channel{Int64}(0))
        out_chan = Distributed.RemoteChannel(()->Channel{Tuple}(0))

        @async begin
            for sigma in 1:fs_num
                put!(in_chan, sigma)
            end
            for p in Distributed.workers()
                put!(in_chan, -1)
            end
        end

        for p in Distributed.workers()
            @async Base.remote_do(
                frag_face_channel, p, in_chan, out_chan, V, EV, FE, sp_idx)
        end

        for sigma in 1:fs_num
            rV, rEV, rFE = skel_merge(rV, rEV, rFE, take!(out_chan)...)
        end

    else

       for sigma in 1:fs_num
           # print(sigma, "/", fs_num, "\r")
           nV, nEV, nFE = Lar.Arrangement.frag_face(
           	V, EV, FE, sp_idx, sigma)
           a,b,c = Lar.skel_merge(
           	rV, rEV, rFE, nV, nEV, nFE)
           global rV=a; global rEV=b; global rFE=c
       end

#		depot_V = Array{Array{Float64,2},1}(undef,fs_num)
#		depot_EV = Array{ChainOp,1}(undef,fs_num)
#		depot_FE = Array{ChainOp,1}(undef,fs_num)
#        for sigma in 1:fs_num
#            print(sigma, "/", fs_num, "\r")
#            nV, nEV, nFE = Arrangement.frag_face( V, EV, FE, sp_idx, sigma)
#            depot_V[sigma] = nV
#            depot_EV[sigma] = nEV
#            depot_FE[sigma] = nFE
#        end
#		rV = vcat(depot_V...)
#		rEV = SparseArrays.blockdiag(depot_EV...)
#		rFE = SparseArrays.blockdiag(depot_FE...)

    end

    rV, rEV, rFE = Lar.Arrangement.merge_vertices(rV, rEV, rFE)

    rCF = Arrangement.minimal_3cycles(rV, rEV, rFE)

    return rV, rEV, rFE, rCF
end

###  2D triangulation
Lar = LinearAlgebraicRepresentation
"""
	obj2lar2D(path::AbstractString)::Lar.LARmodel

Read a *triangulation* from file, given its `path`. Return a `LARmodel` object
"""
function obj2lar2D(path::AbstractString)::Lar.LARmodel
    vs = Array{Float64, 2}(undef, 0, 3)
    edges = Array{Array{Int, 1}, 1}()
    faces = Array{Array{Int, 1}, 1}()

    open(path, "r") do fd
		for line in eachline(fd)
			elems = split(line)
			if length(elems) > 0
				if elems[1] == "v"
					x = parse(Float64, elems[2])
					y = parse(Float64, elems[3])
					z = parse(Float64, elems[4])
					vs = [vs; x y z]
				elseif elems[1] == "f"
					# Ignore the vertex tangents and normals
					v1 = parse(Int, elems[2])
					v2 = parse(Int, elems[3])
					v3 = parse(Int, elems[4])
					append!(edges, map(sort,[[v1,v2],[v2,v3],[v3,v1]]))
					push!(faces, [v1, v2, v3])
				end
				edges = collect(Set(edges))
			end
		end
	end
    return (vs, [edges,faces])::Lar.LARmodel
end

"""
	lar2obj2D(V::Lar.Points,
			cc::Lar.ChainComplex)::String

Produce a *triangulation* from a `LARmodel`. Return a `String` object
"""
function lar2obj2D(V::Lar.Points, cc::Lar.ChainComplex)::String
    @assert length(cc) == 2
    copEV, copFE = cc
    V = [V zeros(size(V, 1))]

    global obj = ""
    for v in 1:size(V, 1)
        	global obj = string(obj, "v ",
        	round(V[v, 1]; digits=6), " ",
        	round(V[v, 2]; digits=6), " ",
        	round(V[v, 3]; digits=6), "\n")
    end

    #triangulated_faces = triangulate2D(V, cc)
    triangulated_faces = triangulate(V, cc)

	obj = string(obj, "\n")
	for f in 1:copFE.m
		triangles = triangulated_faces[f]
		for tri in triangles
			t = tri
			#t = copCF[c, f] > 0 ? tri : tri[end:-1:1]
			obj = string(obj, "f ", t[1], " ", t[2], " ", t[3], "\n")
		end
	end

    return obj
end

#TODO: finish by using a string as an IObuffer
#"""
#	lar2tria2lar(V::Lar.Points, cc::Lar.ChainComplex)::Lar.LARmodel
#
#Return a triangulated `LARmodel` starting from a stadard LARmodel.
#Useful for colour drawing a complex of non-convex cells.
#
#"""
#function lar2tria2lar(V::Lar.Points, cc::Lar.ChainComplex)::Lar.LARmodel
#	obj = Lar.lar2obj2D(V::Lar.Points, cc::Lar.ChainComplex)
#	vs, (edges,faces) = Lar.obj2lar2D(obj::AbstractString)::Lar.LARmodel
#	return (vs, [edges,faces])::Lar.LARmodel
#end

"""
	triangulate2D(V::Lar.Points,
			cc::Lar.ChainComplex)::Array{Any, 1}

Compute a *CDT* for each face of a `ChainComplex`. Return an `Array` of triangles.
"""
function triangulate2D(V::Lar.Points, cc::Lar.ChainComplex)::Array{Any, 1}
    copEV, copFE = cc
    triangulated_faces = Array{Any, 1}(undef, copFE.m)
    if size(V,2)==2
		V = [V zeros(size(V,1),1)]
	end

    for f in 1:copFE.m
        edges_idxs = copFE[f, :].nzind
        edge_num = length(edges_idxs)
        edges = Array{Int64,1}[] #zeros(Int64, edge_num, 2)

        fv = Lar.buildFV(copEV, copFE[f, :])
        vs = V[fv, :]

        for i in 1:length(fv)
        	edge = Int64[0,0]
            edge[1] = fv[i]
            edge[2] = i == length(fv) ? fv[1] : fv[i+1]
            push!(edges,edge::Array{Int64,1})
        end
        edges = hcat(edges...)'
        edges = convert(Array{Int64,2}, edges)

        triangulated_faces[f] = Triangle.constrained_triangulation(
        vs, fv, edges, fill(true, edge_num))
        tV = V[:, 1:2]

        area = Lar.face_area(tV, copEV, copFE[f, :])
        if area < 0
            for i in 1:length(triangulated_faces[f])
                triangulated_faces[f][i] = triangulated_faces[f][i][end:-1:1]
            end
        end
    end

    return triangulated_faces
end

"""
	lar2cop(CV::Lar.Cells)::Lar.ChainOp

Convert an array of array of integer indices to vertices into a sparse matrix.

# Examples

For a single 3D unit cube we get:

```
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> Matrix(Lar.lar2cop(EV))
12×8 Array{Int8,2}:
 1  1  0  0  0  0  0  0
 0  0  1  1  0  0  0  0
 0  0  0  0  1  1  0  0
 0  0  0  0  0  0  1  1
 1  0  1  0  0  0  0  0
 0  1  0  1  0  0  0  0
 0  0  0  0  1  0  1  0
 0  0  0  0  0  1  0  1
 1  0  0  0  1  0  0  0
 0  1  0  0  0  1  0  0
 0  0  1  0  0  0  1  0
 0  0  0  1  0  0  0  1

julia> Matrix(Lar.lar2cop(FV))
6×8 Array{Int8,2}:
 1  1  1  1  0  0  0  0
 0  0  0  0  1  1  1  1
 1  1  0  0  1  1  0  0
 0  0  1  1  0  0  1  1
 1  0  1  0  1  0  1  0
 0  1  0  1  0  1  0  1

julia> Matrix(Lar.lar2cop(CV))
1×8 Array{Int8,2}:
 1  1  1  1  1  1  1  1
```
"""
function lar2cop(CV::Lar.Cells)::Lar.ChainOp
	I = Int64[]; J = Int64[]; Value = Int8[];
	for k=1:size(CV,1)
		n = length(CV[k])
		append!(I, k * ones(Int64, n))
		append!(J, CV[k])
		append!(Value, ones(Int64, n))
	end
	return SparseArrays.sparse(I,J,Value)
end

"""
	cop2lar(cop::Lar.ChainOp)::Lar.Cells

Convert a sparse array of type `ChainOp` into an array of array of type `Cells`.

Notice that `cop2lar` is the inverse function of `lar2cop`. their composition is the identity function.

# Example

```
julia> V,(VV,EV,FV,CV) = Lar.cuboid([1,1,1],true);

julia> Lar.cop2lar(Lar.lar2cop(EV))
12-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
   ...
 [2, 6]
 [3, 7]
 [4, 8]

julia> Lar.cop2lar(Lar.lar2cop(FV))
6-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4]
 [5, 6, 7, 8]
 [1, 2, 5, 6]
 [3, 4, 7, 8]
 [1, 3, 5, 7]
 [2, 4, 6, 8]

julia> Lar.cop2lar(Lar.lar2cop(CV))
1-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4, 5, 6, 7, 8]
```
"""
function cop2lar(cop::Lar.ChainOp)::Lar.Cells
	[findnz(cop[k,:])[1] for k=1:size(cop,1)]
end

function FV2EVs(copEV::Lar.ChainOp, copFE::Lar.ChainOp)
	EV = [findnz(copEV[k,:])[1] for k=1:size(copEV,1)]
	FE = [findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
	EVs = [[EV[e] for e in fe] for fe in FE]
	return EVs
end

"""
	randomcuboids(n,scale=1.0

Generate the `LAR` model of a collection of `n` random cuboids.
Position, orientation and measure of sides are all random.
"""
function randomcuboids(n,scale=1.0)
	assembly = []
	for k=1:n
		corner = rand(Float64, 2)
		sizes = rand(Float64, 2)
		V,(_,EV,_) = Lar.cuboid(corner,true,corner+sizes)
		center = (corner + corner+sizes)/2
		angle = rand(Float64)*2*pi
		obj = Lar.Struct([ Lar.t(center...), Lar.r(angle),
				Lar.s(scale,scale), Lar.t(-center...), (V,EV) ])
		push!(assembly, obj)
	end
	Lar.struct2lar(Lar.Struct(assembly))
end


text:
using Aeroelasticity
using LinearAlgebra
using GXBeam
using ForwardDiff
using Test

include("testutil.jl")

include(joinpath("models", "aerodynamics", "steady.jl"))
include(joinpath("models", "aerodynamics", "quasisteady.jl"))
include(joinpath("models", "aerodynamics", "wagner.jl"))
include(joinpath("models", "aerodynamics", "peters.jl"))
include(joinpath("models", "aerodynamics", "liftingline.jl"))

include(joinpath("models", "structures", "section.jl"))
include(joinpath("models", "structures", "gxbeam.jl"))

include(joinpath("models", "couplings", "steady-section.jl"))
include(joinpath("models", "couplings", "quasisteady-section.jl"))
include(joinpath("models", "couplings", "wagner-section.jl"))
include(joinpath("models", "couplings", "peters-section.jl"))


text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: March, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Uncertainty quantification for AVP for Sleipner via Langevin dynamics

using LinearAlgebra, Statistics, PyPlot, Images, JLD2
push!(LOAD_PATH, "./src/"); using WaveOpsAVP
push!(LOAD_PATH, "../LangevinDynamics/src/"); using LangevinDynamics
using Random; Random.seed!(1)

## Load model, synthetics, and MAP result

@load "./data/Sleipner/Sleipner_data.jld"
@load "./data/Sleipner/results/MAP_results.jld"

## Setting log-posterior loss

# Pre- and post-processing for log-posterior
nz = size(Mtrue.m, 1)
mmin = R(1/3000^2)
mmax = Inf
function proj_bounds(m::Array{R}, mmin, mmax)
    m_ = copy(m)
    m_[m_ .< mmin] .= mmin
    m_[m_ .> mmax] .= mmax
    return m_
end
Mbg_ = Mbg
# Mbg_ = Mmap
# Mbg_ = Mtrue

# Log-likelihood
mloglike_fun(m::Model; compute_grad::Bool = false) = mLogLikelihood_gauss(dat_n, m, σ2_n, F; compute_grad = compute_grad)

# Log-prior
D = derivative_linop(Mbg.geom)
# D = LinFunctionMod2Mod(m->m, m->m)
M0 = Model(Mbg.geom, zeros(R, nz, 1))
# M0 = Mbg
Dx = D.eval(Mtrue-M0)
# (gauss)
# σ2_pr = R(2e-8)^2
σ2_pr = R(1e-8)^2
mlogpr_fun(m::Model; compute_grad::Bool = false) = mLogPrior_gauss(m, σ2_pr; m0 = M0, D = D, compute_grad = compute_grad)
# # (L1-smooth)
# ϵ2 = R(0.001)*norm(Dx.m, Inf)^2
# σ2_pr = R(1e-3)*sum(sqrt.(Dx.m.^R(2).+ϵ2))[1]
# mlogpr_fun(m::Model; compute_grad::Bool = false) = mLogPrior_L1smooth(m, σ2_pr, ϵ2; m0 = M0, D = D, compute_grad = compute_grad)

# Log-posterior
function mlogpost_fun(m::Array{R, 2})
    f, g = mLogPost(Model(Mbg.geom, proj_bounds(m, mmin, mmax)), mloglike_fun, mlogpr_fun; compute_grad = true)
    return f[1], g.m
end

## Langevin sampling

# SGLD parameters
maxit = 2^10
ϵ0 = R(1e-17)
ϵ1 = ϵ0*R(1e-2)
gamma = R(0.55)
# gamma = R(0)
burnin = 100
m0 = Mmap.m

# Run SGLD
time = @elapsed fval, μ, σ2, ϵ, minv = langevinSolver(mlogpost_fun, maxit, m0; ϵ0 = ϵ0, ϵ1 = ϵ1, gamma = gamma, burnin = burnin, save = true)

# Saving result
@save "./data/Sleipner/results/UQ_langevin_results.jld" dat_n σ2_pr σ2_n fval minv time μ σ2 ϵ