text:
module Traits
@doc """This package provides an implementation of traits, aka interfaces or type-classes.
     It is based on the premises that traits are:

     - contracts on one type or between several types.  The contract can
        contain required methods but also other assertions and just
        belonging to a group (i.e. the trait).
     - they are structural types: i.e. they needn't be declared explicitly
        """ -> current_module()

export istrait, istraittype, issubtrait, check_return_types,
       traitgetsuper, traitgetpara, traitmethods,
       @traitdef, @traitimpl, @traitfn, TraitException, All

if !(VERSION>v"0.4-")
    error("Traits.jl needs Julia version 0.4.-")
end

## patches for bugs in base
include("base_fixes.jl")

## common helper functions
include("helpers.jl")

#######
# Flags
#######
type Flags
    # A type to hold flags
    check_return_types::Bool
end
getflag(f::Flags, n::Symbol) = getfield(f, n)
setflag!(f::Flags, n::Symbol, val) = setfield!(f, n, val)

"""
Flags for various Traits behavior:

- first arg: check return types?
"""
const flags = Flags(true)

# By setting them in Main before using, they can be turned on or off.
# TODO: update to use functions.
if isdefined(Main, :Traits_check_return_types)
    println("Traits.jl: not using return types of @traitdef functions")
    setflag!(flags, :check_return_types, Main.Traits_check_return_types)
end

@doc "Toggles return type checking.  Will issue warning because of const declaration, ignore:"->
function check_return_types(flg::Bool)
    setflag!(flags, :check_return_types, flg)
    nothing
end

#######
# Types
#######
@doc """`abstract Trait{SUPER}`

     All traits are direct decedents of abstract type Trait.  The type parameter
     SUPER of Trait is needed to specify super-traits (a tuple).""" ->
abstract Trait{SUPER<:Tuple}

# Type of methods field of concrete traits:
typealias FDict Dict{Union{Function,DataType},Function}

# A concrete trait type has the form
## Tr{X,Y,Z} <: Trait{ Tuple{ST1{X,Y},ST2{Z}} }
#
# immutable Tr1{X1} <: Traits.Trait{ Tuple{} }
#     methods::FDict
#     constraints::Vector{Bool}
#     assoctyps::Vector{Any}
#     Tr1() = new(FDict(methods_defined), Bool[], [])
# end
#
# where methods field holds the function signatures, like so:
# Dict{Function,Any} with 3 entries:
#   start => _start(Int64) = Tuple{Vararg{Any}}
#   next  => _next(Int64,Any) = Tuple{Vararg{Any}}
#   done  => _done(Int64,Any) = Bool

# used to dispatch to helper methods
immutable _TraitDispatch end
immutable _TraitStorage end

# @doc """Type All is to denote that any type goes in type signatures in
#      @traitdef.  This is a bit awkward:

#      - method_exists(f, s) returns true if there is a method of f with
#        signature sig such that s<:sig.  Thus All<->Union{}
#      - Base.return_types works the other way around, there All<->Any

#      See also https://github.com/JuliaLang/julia/issues/8974"""->
# abstract All

# General trait exception
type TraitException <: Exception
    msg::AbstractString
end

# Helper dummy types used in istrait below
abstract _TestType{T}
immutable _TestTvar{T}<:_TestType{T} end # used for TypeVar testing
#Base.show{T<:_TestType}(io::IO, x::Type{T}) = print(io, string(x.parameters[1])*"_")

#########
# istrait, one of the core functions
#########

# Update after PR #10380
@doc """Tests whether a DataType is a trait.  (But only istrait checks
        whether it's actually full-filled)""" ->
istraittype(x) = false
istraittype{T<:Trait}(x::Type{T}) = true
function istraittype{T<:Tuple}(xs::Type{T})
    for x in xs
        istraittype(x) || return false
    end
    return true
end

@doc """Tests whether a set of types fulfill a trait.
     A Trait Tr is defined for some parameters if:

     - all the functions of a trait are defined for them
     - all the trait constraints are fulfilled

     Example:

     `istrait(Tr{Int, Float64})`

     or with a tuple of traits:

     `istrait( (Tr1{Int, Float64}, Tr2{Int}) )`
     """ ->
function istrait{T<:Trait}(Tr::Type{T}; verbose=false)
    if verbose
        println_verb(x) = println("**** Checking $(deparameterize_type(Tr)): " * x)
    else
        println_verb = x->x
    end

    if !hasparameters(Tr)
        throw(TraitException("Trait $Tr has no type parameters."))
    end
    # check supertraits
    !istrait(traitgetsuper(Tr); verbose=verbose) && return false

    # check instantiating
    tr = nothing
    try
        tr = Tr()
    catch err
        println_verb("""Could not instantiate instance for type encoding the trait $Tr.
                     This usually indicates that something is amiss with the @traitdef
                     or that one of the generic functions is not defined.
                     The error was: $err""")
        return false
    end

    # check constraints
    if !reduce(&, tr.constraints)
        println_verb("Not all constraints are satisfied for $T")
        return false
    end

    # Check call signature of all methods:
    for (gf,_gf) in tr.methods
        println_verb("*** Checking function $gf")
        # Loop over all methods defined for each function in traitdef
        for tmm in methods(_gf)
            tm = FakeMethod(tmm, ret=true)
            println_verb("** Checking method $tmm")
            checks = false
            # Only loop over methods which have the right number of arguments:
            for fm in methods(gf, NTuple{length(tm.sig),Any})
                if isfitting(tmm, fm, verbose=verbose)
                    checks = true
                    break
                end
            end
            if !checks # if check==false no fitting method was found
                println_verb("""No method of the generic function/call-overloaded `$gf` matched the
                             trait specification: `$tm`""")
                return false
            end
        end
    end

    # check return-type.  Specifed return type tret and return-type of
    # the methods frets should fret<:tret.  This is backwards to
    # argument types checking above.
    if flags.check_return_types
        for (gf,_gf) in tr.methods
            println_verb("*** Checking return types of function $gf")
            for tmm in methods(_gf) # loop over all methods defined for each function in traitdef
                tm = FakeMethod(tmm, ret=true)
                tm = replace_concrete_tvars(tm)
                println_verb("** Checking return types of method $tmm")
                tret_typ = tm.ret
                fret_typ = Base.return_types(gf, tm.sig)
                # at least one of the return types need to be a subtype of tret_typ
                checks = false
                for fr in fret_typ
                    if fr<:tret_typ
                        checks = true
                    end
                end
                if !checks
                    println_verb("""For function $gf: no return types found which are subtypes of the specified return type:
                                 $tret_typ
                                 List of found return types:
                                 $fret_typ
                                 Returning false.
                                 """)
                    return false
                end
            end
        end
    end
    return true
end
# check a tuple of traits against a signature
function istrait{T<:Tuple}(Trs::Type{T}; verbose=false)
    for Tr in Trs
        istrait(Tr; verbose=verbose) || return false
    end
    return true
end

## Helpers for istrait

# FakeMethod is a datatype similar to Method.
immutable FakeMethod
    sig::Type
    tvars::SimpleVector
    va::Bool
    ret::Union{Type,TypeVar} # return type
    FakeMethod(sig,tv,va) = new(sig,tv,va)
    FakeMethod(sig,tv,va,ret) = new(sig,tv,va,ret)
end
function FakeMethod(m::Method; ret=false)
    # make sure the tvars are always in a SimpleVec
    tvs = isa(m.tvars, SimpleVector) ? m.tvars : Base.svec(m.tvars)
    if ret
        # last
        sig = Tuple{m.sig.parameters[2:end]...}
        ret = m.sig.parameters[1]
        return FakeMethod(sig, tvs, m.va, ret)
    else
        return FakeMethod(m.sig, tvs, m.va)
    end
end
@doc """isfitting checks whether the signature of a method `tm`
     specified in the trait definition is fulfilled by one method `fm`
     of the corresponding generic function.  This is the core function
     which is called by istraits.

     Checks that tm.sig<:fm.sig and that the parametric constraints on
     fm and tm are equal where applicable.  Lets call this relation tm<<:fm.

     So, summarizing, for a trait-signature to be satisfied (fitting)
     the following condition need to hold:

     A) `tsig<:sig` for just the types themselves (sans parametric
         constraints)

     B) The parametric constraints parameters on `sig` and `tsig` need
        to feature in the same argument positions.  Except when the
        corresponding function parameter is constraint by a concrete
        type: then make sure that all the occurrences are the same
        concrete type.

     Examples, left trait-method, right implementation-method:
     {T<:Real, S}(a::T, b::Array{T,1}, c::S, d::S) <<: {T<:Number, S}(a::T, b::AbstractArray{T,1}, c::S, d::S)
     -> true

     {T<:Integer}(T, T, Integer) <<: {T<:Integer}(T, T, T)
     -> false as parametric constraints are not equal
     """ ->
function isfitting(tmm::Method, fmm::Method; verbose=false) # tm=trait-method, fm=function-method
    println_verb = verbose ? println : x->x # TODO maybe move this function out
    tm = FakeMethod(tmm, ret=true)
    fm = FakeMethod(fmm)

    # Replace type parameters which are constraint by a concrete type
    # (because Vector{TypeVar(:V, Int)}<:Vector{Int}==false but we need ==true)
    tm = replace_concrete_tvars(tm)
    fm = replace_concrete_tvars(fm)

    # Special casing call-overloading:
    if fmm.func.code.name==:call && tmm.func.code.name!=:call # true if only fm is call-overloaded
        # prepend ::Type{...} to signature
        tm = FakeMethod(Tuple{fm.sig[1], tm.sig...}, tm.tvars, tm.va)
        # check whether there are method parameters too:
        for ftv in fm.tvars
            flocs = find_tvar(fm.sig, ftv)
            if flocs[1]==1 # yep, has a constraint like call{T}(::Type{Array{T}},...)
                if length(flocs)==1
                    tm = FakeMethod(tm.sig, Base.svec(ftv, tm.tvars...) , tm.va)
                else
                    println_verb("This check is not implemented, returning false.")
                    # Note that none of the 1000 methods of call in
                    # Base end up here.
                    return false
                end
            end
        end
        # # There is a strange bug which is prevented by this never
        # # executing @show.  I'll try and investigate this in branch
        # # m3/heisenbug
        # if length(tm.sig)==-10
        #     @show tm
        #     error("This is not possible")
        # end
    end

    ## Check condition A:
    # If there are no function parameters then just compare the
    # signatures.
    if tm.tvars==Base.svec() && fm.tvars==Base.svec()
        out = tm.sig<:fm.sig
        println_verb("Reason fail/pass: no tvars in trait-method only checking signature. Result: $out")
        return out
    end

    # If !(tm.sig<:fm.sig) then tm<<:fm is false
    # but the converse is not true:
    if !(tm.sig<:fm.sig)
        println_verb("""Reason fail: !(tm.sig<:fm.sig)
                     tm.sig = $(tm.sig)
                     fm.sig = $(fm.sig)""")
        return false
    end

    # False if there are not the same number of arguments: (I don't
    # think this test is necessary as it is tested above.)
    if length(tm.sig)!=length(fm.sig)!
        println_verb("Reason fail: not same argument length.")
        return false
    end
    # Getting to here means that that condition (A) is fulfilled.

    ## Check condition B:
    # If there is only one argument then we're done as parametric
    # constraints play no role:
    if length(tm.sig)==1
        println_verb("Reason pass: length(tm.sig)==1")
        return true
    end

    # First special case if tm.tvars==Base.svec() && !(fm.tvars==Base.svec())
    if tm.tvars==Base.svec()
        fm.tvars==Base.svec() && error("Execution shouldn't get here as this should have been checked above!")
        for (i,ftv) in enumerate(fm.tvars)
            # If all the types in tm.sig, which correspond to a
            # parameter constraint argument of fm.sig, are the same then pass.
            typs = Tuple{tm.sig.parameters[find_tvar(fm.sig, ftv)]...}
            if length(typs)==0
                println_verb("Reason fail: this method $fmm is not callable because the static parameter does not occur in signature.")
                return false
            elseif length(typs)==1 # Necessarily the same
                continue
            else # length(typs)>1
                if !reduce(&, map(isleaftype, typs)) # note isleaftype can have some issues with inner constructors
                    println_verb("Reason fail: not all parametric-constraints in function-method $fmm are on leaftypes in traitmethod $tmm.")
                    return false
                else
                    # Now check that all of the tm.sig-types have the same type at the parametric-constraint sites.
                    if !allequal(find_correponding_type(tm.sig, fm.sig, ftv))
                        println_verb("Reason fail: not all parametric-constraints in function-method $fmm correspond to the same type in traitmethod $tmm.")
                        return false
                    end
                end
            end
        end

        println_verb("""Reason pass: All occurrences of the parametric-constraint in $fmm correspond to the
                        same type in trait-method $tmm.""")
        return true
    end

    # Strategy: go through constraints on trait-method and check
    # whether they are fulfilled in function-method.
    for tv in tm.tvars
        # find all occurrences in the signature
        locs = find_tvar(tm.sig, tv)
        if length(locs)==0
            throw(TraitException("""The parametric-constraint of trait-method $tmm has to feature
                                 in at least one argument of the signature."""))
        end
        # Find the tvar in fm which corresponds to tv.
        ftvs = Any[]
        for ftv in fm.tvars
            flocs = find_tvar(fm.sig, ftv)
            if flocs==locs
                push!(ftvs,ftv)
            end
        end

        if length(ftvs)==0
            ## This should pass, because the trait-parameter is a leaftype:
            # @traitdef Tr01{X} begin
            #     g01{T<:X}(T, T) -> T
            # end
            # g01(::Int, ::Int) = Int
            # @assert istrait(Tr01{Int}, verbose=true)
            if isleaftype(tv.ub) # note isleaftype can have some issues with inner constructors
                # Check if the method definition of fm has the same
                # leaftypes in the same location.
                if mapreduce(x -> x==tv.ub, &, true, fm.sig[locs])
                    println_verb("Reason pass: parametric constraints only on leaftypes.")
                    return true
                end
            end
            println_verb("Reason fail: parametric constraints on function method not as severe as on trait-method.")
            return false
        end

        # Check that they constrain the same thing in each argument.
        # E.g. this should fail: {K,V}(::Dict{K,V}, T) <<: {T}(::Dict{V,K}, T).
        # Do this by substituting a concrete type into the respective
        # TypeVars and check that arg(tv')<:arg(ftv')
        checks = false
        for ft in ftvs
            for i in find(locs)
                targ = subs_tvar(tv, tm.sig[i], _TestTvar{i})
                farg = subs_tvar(ft, fm.sig[i], _TestTvar{i})
                checks = checks || (targ<:farg)
            end
        end
        if !checks
            println_verb("Reason fail: parametric constraints on args $(tm.sig[i]) and $(fm.sig[i]) on different TypeVar locations!")
            return false
        end
    end

    println_verb("Reason pass: all checks passed")
    return true
end

# helpers for isfitting
function subs_tvar(tv::TypeVar, arg::DataType, TestT::DataType)
    # Substitute `TestT` for a particular TypeVar `tv` in an argument `arg`.
    #
    # Example:
    # Array{I<:Int64,N} -> Array{_TestTvar{23},N}
    if isleaftype(arg) || length(arg.parameters)==0 # concrete type or abstract type with no parameters
        return arg
    else # It's a parameterized type: do substitution on all parameters:
        pa = [ subs_tvar(tv, arg.parameters[i], TestT) for i=1:length(arg.parameters) ]
        typ = deparameterize_type(arg)
        return typ{pa...}
    end
end
subs_tvar(tv::TypeVar, arg::TypeVar, TestT::DataType) = tv===arg ? TestT : arg  # note === this it essential!
subs_tvar(tv::TypeVar, arg, TestT::DataType) = arg # for anything else

function replace_concrete_tvars(m::FakeMethod)
    # Example:
    # FakeMethod(Tuple{T<:Int64,Array{T<:Int64,1},Integer},svec(T<:Int64,),false, T<:Int64)
    # ->
    # FakeMethod(Tuple{Int64,   Array{Int64,1},   Integer},svec()         ,false, Int64)
    newtv = []
    newsig = Any[m.sig.parameters...] # WTF: need the .parameters... here even though
                                      #      I added the start, next & done methods to Type{Tuple}
    if isdefined(m, :ret)
        ret_tuple = isa(m.ret, Tuple)
        newret = ret_tuple ? Any[m.ret.parameters...] : Any[m.ret]
    end
    for tv in m.tvars
        if !isleaftype(tv.ub)
            push!(newtv, tv)
        else
            replaced = true
            newsig = Any[subs_tvar(tv, arg, tv.ub) for arg in newsig]
            if isdefined(m, :ret)
                newret = Any[subs_tvar(tv, arg, tv.ub) for arg in newret]
            end
        end
    end
    if isdefined(m, :ret)
        if ret_tuple
            return FakeMethod(Tuple{newsig...}, Base.svec(newtv...), m.va, Tuple{newret...})
        else
            return FakeMethod(Tuple{newsig...}, Base.svec(newtv...), m.va, newret[1])
        end
    else
        return FakeMethod(Tuple{newsig...}, Base.svec(newtv...), m.va)
    end
end

# Finds the types in tmsig which correspond to TypeVar ftv in fmsig
function find_correponding_type(tmsig::Tuple, fmsig::Tuple, ftv::TypeVar)
    out = Any[]
    for (ta,fa) in zip(tmsig,fmsig)
        if isa(fa, TypeVar)
            fa===ftv && push!(out, ta)
        elseif isa(fa, DataType) || isa(fa, Tuple)
            append!(out, find_correponding_type(ta,fa,ftv))
        else
            println("ta=$ta\nfa=$fa")
            error("Not implemented")
        end
    end
    return out
end
function find_correponding_type(ta::DataType, fa::DataType, ftv::TypeVar)
    # gets here if fa is not a TypeVar
    out = Any[]
    if !( deparameterize_type(ta)<:deparameterize_type(fa)) # ||
        # length(ta.parameters)!=length(fa.parameters)  # don't check for length.  If not the same length, assume that the first parameters are corresponding...
        push!(out, _TestType{:no_match})  # this will lead to a no-match in isfitting
        return out
    end
    for (tp,fp) in zip(ta.parameters,fa.parameters)
        if isa(fp, TypeVar)
            fp===ftv && push!(out, tp)
        elseif isa(fp, DataType) || isa(fa, Tuple)
            append!(out, find_correponding_type(tp,fp,ftv))
        end
    end
    return out
end

# find_tvar finds index of arguments in a function signature `sig` where a
# particular TypeVar `tv` features. Example:
#
# find_tvar( Tuple{T, Int, Array{T}} -> [1,3]
function find_tvar{T<:Tuple}(sig::Type{T}, tv)
    ns = length(sig)
    out = Int[]
    for i = 1:ns
        if length(find_tvar(sig[i], tv))>0
            push!(out,i)
        end
    end
    return out
end
find_tvar(sig::TypeVar, tv) = sig===tv ? [1] : Int[]  # note ===, this is essential!
function find_tvar(arg::DataType, tv)
    ns = length(arg.parameters)
    for i=1:ns
        if length(find_tvar(arg.parameters[i], tv))>0
            return [1]
        end
    end
    return Int[]
end
find_tvar(arg, tv) = Int[]

######################
# Sub and supertraits:
######################
@doc """Returns the super traits""" ->
traitgetsuper{T<:Trait}(t::Type{T}) =  t.super.parameters[1]
traitgetpara{T<:Trait}(t::Type{T}) =  Tuple{t.parameters...}

@doc """Checks whether a trait, or a tuple of them, is a subtrait of
     the second argument.""" ->
function issubtrait{T1<:Trait,T2<:Trait}(t1::Type{T1}, t2::Type{T2})
    if t1==t2
        return true
    end
    if t2 in traitgetsuper(t1)
        return true
    end
    for t in traitgetsuper(t1)
        issubtrait(t, t2) && return true
    end
    return false
end

# TODO: think about how to handle tuple traits and empty traits
function issubtrait{T1<:Trait, Tu<:Tuple}(t1::Type{T1}, t2::Type{Tu})
    if t2==()
        # the empty trait is the super-trait of all traits
        true
    else
        throw(TraitException(""))
    end
end

# traits in a tuple have no order, really, this should reflected.
# Maybe use a set instead?  Subtrait if it is a subset?
function issubtrait{T1<:Tuple, T2<:Tuple}(t1::Type{T1}, t2::Type{T2})
    if length(t1)!=length(t2)
        return false
    end
    checks = true
    for (p1,p2) in zip(t1, t2)
        checks = checks && issubtrait(p1,p2)
    end
    return checks
end

## Trait definition
include("traitdef.jl")

# Trait implementation
include("traitimpl.jl")

# Trait functions
include("traitfns.jl")

## Common traits
include("commontraits.jl")

end # module


text:
#=

using Static
using Base: RefValue

_elcheck(::Type{T1}, ::Type{T2}) where {T1,T2<:Ref} = _elcheck(T, eltype(T2))
_elcheck(::Type{T1}, ::Type{T2}) where {T1,T2} = __elcheck(is_static(T2), T, T2)
__elcheck(::True, ::Type{T1}, ::Type{T2}) where {T1,T2} = _elcheck(T1, eltype(T2))
_elcheck(::Type{T}, ::Type{T}) where {T} = True()
_elcheck(::Type{T1}, ::Type{T2}) where {T1,T2} = False()

_check_elconvert(::Type{T1}, x::T2) where {T1,T2} = __check_elconvert(_elcheck(T1, T2), T1, x)
__check_elconvert(::True, ::Type{T}, x) where {T} = x
__check_elconvert(::False, ::Type{T}, x) where {T} = _elconvert(T, x)
_elconvert(::Type{T}, x::Ref) where {T} = Ref{T}()
_elconvert(::Type{T1}, x::T2) where {T1,T2} = __elconvert(is_static(T2), T1, x)
__elconvert(::True, ::Type{T1}, x::T2) where {T1,T2} = static(convert(T1, x))
__elconvert(::False ::Type{T1}, x::T2) where {T1,T2} = convert(T1, x)

dynamic(x::Base.RefValue) = x[]

function _unit_range_last!(start, stop::RefValue)
    stop[] = _unit_range_last!(dynamic(start), dynamic(stop))
    return stop
end
_unit_range_last!(start, stop) = Base.unitrange_last(dynamic(start), stop)

const TRef{T} = Union{T,RefValue{T}}

struct URange{T,B<:TRef{T},E<:TRef{T}} <: AbstractUnitRange{T}
    start::B  # B for begin
    stop::E   # E for end

    global _URange(start::B, stop::E) where {B,E} = new{eltype(start),B,E}(start, stop)
end

URange{Int,One,Int}(start::One, stop::Int) = _URange(start, max(0, stop))
URange{Int,One,Int}(start::One, stop::Int) = _URange(start, max(0, stop))
function URange{T,B,E}(start::B, stop::E) where {T,B,E}
    @assert known(_elcheck(T, B))
    @assert known(_elcheck(T, E))
    return _URange(start, _unit_range_last!(start, stop))
end
function URange{T}(start::B, stop::E) where {T,B,E}
    start2 = _check_elconvert(T, start2)
    return _URange(start2, _unit_range_last!(start2, _check_elconvert(T, stop)))
end
function URange(start::B, stop::E) where {B,E}
    return URange{typejoin(eltype(B),eltype(E))}(start, stop)
end

const Length{L} = URange{Int,One,L}
Length(stop::Int) = _URange(One(), max(0, stop))
Length(stop::StaticInt) = _URange(One(), max(static(0), stop))
Length(stop::Ref{Int}) = _URange(One(), max(static(0), stop))

const StaticLength{L} = Length{StaticInt{L}}
const MutableLength = Length{RefValue{Int}}
const DynamicLength = Length{Int}

can_grow_tail(::Type{URange{T,B,E}}) where {T,B,E} = ismutable(E)
can_grow_front(::Type{URange{T,B,E}}) where {T,B,E} = ismutable(B)
function ArrayInterface.can_grow_tail(::Type{T}) wehre {T<:URange}
    return can_grow_tail(T) || can_grow_front(T)
end

can_grow_tail(::Type{T}) where {T<:}

function grow_tail!()
end

# LinRange
struct LengthRange{T,R<:URange{T},L} <: AbstractRange{T}
    span::R
    length::Length{L}
end

struct DynamicStepRange{T,R<:URange{T},S} <: AbstractRange{T}
    span::R
    step::S
end

=#

ext: jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    # slider
    sl = slider(1:3)
    sl[] = (1:5, 3)
    precompile(sl)
    destroy(sl)

    # checkbox
    cb = checkbox(true)
    cb[] = false
    precompile(cb)
    destroy(cb)

    # togglebutton
    tb = togglebutton(true)
    tb[] = false
    precompile(tb)
    destroy(tb)

    # button
    btn = button("Push me")
    btn[] = nothing
    precompile(btn)
    destroy(btn)

    # textbox
    tb1, tb2 = textbox("Edit me"), textbox(3; range=1:5)
    tb1[] = "done"
    tb2[] = 4
    precompile(tb1)
    precompile(tb2)
    destroy(tb1); destroy(tb2)

    # textarea
    ta = textarea("Lorem ipsum...")
    ta[] = "...muspi meroL"
    precompile(ta)
    destroy(ta)

    # dropdown
    dd = dropdown(["one", "two", "three"])
    precompile(dd)
    destroy(dd)
    @assert precompile(Tuple{typeof(dropdown),Tuple{Vararg{String, 100}}})   # time: 0.042523906
    @assert precompile(Tuple{typeof(dropdown),Vector{Pair{String, Function}}})   # time: 0.01169991

    # label
    lbl = label("Info")
    lbl[] = "other info"
    precompile(lbl)
    destroy(lbl)
    @assert precompile(Tuple{typeof(show),IOBuffer,Label})   # time: 0.07879969

    # spinbutton
    sb = spinbutton(1:3)
    sb[] = (1:4, 4)
    precompile(sb)
    destroy(sb)

    # cyclicspinbutton
    csb = cyclicspinbutton(1:3, Observable(true))
    csb[] = 4
    precompile(csb)
    destroy(csb)

    # progressbar
    pb = progressbar(1:3)
    pb[] = 2
    precompile(pb)
    destroy(pb)
    pb = progressbar(1.0 .. 3.0)
    pb[] = 2.2
    precompile(pb)
    destroy(pb)

    # player
    p = player(1:3)
    precompile(p)
    p[] = 2
    destroy(p)

    # # timewidget
    # tw = timewidget(Dates.Time(1,1,1))
    # tw[] = Dates.Time(2,2,2)
    # precompile(tw)
    # destroy(tw)

    # # datetimewidget
    # dtw = datetimewidget(DateTime(1,1,1,1,1,1))
    # dtw[] = DateTime(2,2,2,2,2,2)
    # destroy(dtw)

    # canvas
    @assert precompile(Tuple{typeof(canvas),Int,Int})   # time: 0.08071636
    @assert precompile(Tuple{typeof(canvas)})   # time: 0.00260139
    for U in (UserUnit, DeviceUnit)
        @assert precompile(Tuple{typeof(canvas),Type{U}})   # time: 0.020233927
        @assert precompile(Tuple{typeof(canvas),Type{U},Int,Int})   # time: 0.020233927
        @assert precompile(Tuple{Type{GtkAspectFrame},Canvas{U},String,Float64,Float64,Float64})   # time: 0.021899309
        @assert precompile(Tuple{typeof(gc_preserve),GtkWindowLeaf,Canvas{U}})

        # @assert precompile(Tuple{typeof(draw),Function,Canvas{U}})   # time: 0.00281273
        # @assert precompile(Tuple{typeof(draw),Function,Canvas{U},Observable{Any}})   # time: 0.006168176
        # @assert precompile(Tuple{typeof(draw),Function,Canvas{U},Observable{Any},Vararg{Observable{Any}, 100}})   # time: 0.007357308

        @assert precompile(Tuple{typeof(fill!),Canvas{U},RGB{N0f8}})   # time: 0.002759036
        @assert precompile(Tuple{typeof(fill!),Canvas{U},RGBA{N0f8}})

        @assert precompile(Tuple{Type{MouseButton{U}}})   # time: 0.001963396
        @assert precompile(Tuple{Type{MouseHandler{U}},GtkCanvas})   # time: 0.017885247
        @assert precompile(Tuple{typeof(mousedown_cb),Ptr{GObject},Ptr{Gtk.GdkEventButton},MouseHandler{U}})
        @assert precompile(Tuple{typeof(mouseup_cb),Ptr{GObject},Ptr{Gtk.GdkEventButton},MouseHandler{U}})
        @assert precompile(Tuple{typeof(mousescroll_cb),Ptr{GObject},Ptr{Gtk.GdkEventScroll},MouseHandler{U}})   # time: 0.018706173
        @assert precompile(Tuple{typeof(mousemove_cb),Ptr{GObject},Ptr{Gtk.GdkEventMotion},MouseHandler{U}})   # time: 0.016197048
        @assert precompile(Tuple{typeof(rubberband_move),Canvas{U},RubberBand{U},MouseButton{U},Cairo.CairoContext})   # time: 0.001334047

        @assert precompile(Tuple{typeof(set_coordinates),Canvas{U},BoundingBox})   # time: 0.003914427
        @assert precompile(Tuple{typeof(set_coordinates),Canvas{U},ZoomRegion{RInt}})   # time: 0.002466032
    end
    @assert precompile(Tuple{typeof(convertunits),Type{DeviceUnit},Canvas{UserUnit},UserUnit,UserUnit})   # time: 0.002624784
    @assert precompile(Tuple{typeof(convertunits),Type{UserUnit},Canvas{UserUnit},DeviceUnit,DeviceUnit})   # time: 0.001084311

    # image_surface & copy!
    for T in (Gray{N0f8}, GrayA{N0f8}, Gray{Float32}, Gray{Float64},
              RGB{N0f8}, RGBA{N0f8}, RGB{Float32}, RGBA{Float32}, RGB{Float64}, RGBA{Float64})
        @assert precompile(Tuple{typeof(image_surface),Matrix{T}})   # time: 0.046262104 etc.
        for U in (UserUnit, DeviceUnit)
            @assert precompile(Tuple{typeof(copy!),Canvas{U},Matrix{T}})   # time: 0.09909209
        end
    end

    # zoom & friends
    @assert precompile(Tuple{typeof(zoom),ZoomRegion{RInt},Float64})   # time: 0.018138554
    for U in (UserUnit, DeviceUnit)
        @assert precompile(Tuple{Type{XY{U}},Float64,Float64})   # time: 0.010058553
        @assert precompile(Tuple{typeof(zoom),ZoomRegion{RInt},Float64,XY{U}})   # time: 0.002209134
    end
    @assert precompile(Tuple{typeof(zoom),ZoomRegion{RInt},Int})   # time: 0.001822723
    @assert precompile(Tuple{Type{ZoomRegion},Tuple{UnitRange{Int}, UnitRange{Int}}})   # time: 0.009671802
    @assert precompile(Tuple{Type{ZoomRegion},Tuple{Base.OneTo{Int}, Base.OneTo{Int}}})   # time: 0.006405025
    @assert precompile(Tuple{Type{ZoomRegion},Matrix{RGB{N0f8}}})   # time: 0.001308999
    @assert precompile(Tuple{Type{ZoomRegion},Tuple{UnitRange{Int}, UnitRange{Int}},Tuple{UnitRange{Int}, UnitRange{Int}}})   # time: 0.001152972
    @assert precompile(Tuple{Type{ZoomRegion},XY{ClosedInterval{RInt}},BoundingBox})   # time: 0.001061675

    @assert precompile(Tuple{typeof(setindex!),Observable{ZoomRegion{RInt}},Tuple{UnitRange{Int}, UnitRange{Int}}})   # time: 0.003122286
    @assert precompile(Tuple{typeof(setindex!),Observable{ZoomRegion{RInt}},XY{ClosedInterval{Int}}})   # time: 0.001308547
    @assert precompile(Tuple{Type{BoundingBox},XY{ClosedInterval{Int}}})   # time: 0.001497522

    # pan & friends
    @assert precompile(Tuple{typeof(pan),ClosedInterval{Int},Float64,ClosedInterval{Int}})
    @assert precompile(Tuple{typeof(pan_x),ZoomRegion{RInt},Float64})   # time: 0.001364175
    @assert precompile(Tuple{typeof(pan_y),ZoomRegion{RInt},Float64})   # time: 0.001364175

    # advanced interaction
    for U in (UserUnit, DeviceUnit), T in (RInt, Float64)
        @assert precompile(Tuple{typeof(init_zoom_rubberband),Canvas{UserUnit},Observable{ZoomRegion{RInt}}})   # time: 0.110777095
        @assert precompile(Tuple{typeof(init_zoom_rubberband),Canvas{UserUnit},Observable{ZoomRegion{RInt}},typeof(GtkObservables.zrb_init_default), typeof(GtkObservables.zrb_reset_default), Int})   # time: 0.110777095
        @assert precompile(Tuple{typeof(init_zoom_scroll),Canvas{UserUnit},Observable{ZoomRegion{RInt}}})   # time: 0.01061211
        @assert precompile(Tuple{typeof(init_pan_drag),Canvas{UserUnit},Observable{ZoomRegion{RInt}}})   # time: 0.004753301
        @assert precompile(Tuple{typeof(init_pan_scroll),Canvas{UserUnit},Observable{ZoomRegion{RInt}}})   # time: 0.00264137
        # Special handling of the interactives to precompile the callback functions
        zr = T == RInt ? ZoomRegion(rand(3,3)) : (xy = XY(T(1) .. T(3), T(1) .. T(3)); ZoomRegion(xy, xy))
        guidict = init_pan_scroll(canvas(U), Observable(zr))
        @assert precompile(guidict["pan"])
        for f in (init_pan_drag, init_zoom_rubberband)
            guidict = f(canvas(U), Observable(zr))
            @assert precompile(guidict["init"])
            @assert precompile(guidict["drag"])
            @assert precompile(guidict["finish"])
        end
        guidict = init_zoom_scroll(canvas(U), Observable(zr))
        @assert precompile(guidict["zoom"])
    end

    # misc
    @assert precompile(Tuple{typeof(axes),ZoomRegion{RInt}})   # time: 0.012924676
    @assert precompile(Tuple{typeof(+),XY{Float64},XY{Float64}})   # time: 0.011136752
    @assert precompile(Tuple{typeof(show),IOBuffer,XY{Float64}})   # time: 0.001530783
    @assert precompile(Tuple{typeof(show),IOBuffer,XY{Int}})   # time: 0.001503061
    @assert precompile(Tuple{typeof(convert),Type{XY{RInt}},XY{Float64}})   # time: 0.001046295
end


text:
const VERSION_KEY = "__versions__"

const VERSIONS = Dict("Julia" => string(VERSION), "BenchmarkTools" => string(BENCHMARKTOOLS_VERSION))

type ParametersPreV006
    seconds::Float64
    samples::Int
    evals::Int
    overhead::Int
    gctrial::Bool
    gcsample::Bool
    time_tolerance::Float64
    memory_tolerance::Float64
end

type TrialPreV006
    params::Parameters
    times::Vector{Int}
    gctimes::Vector{Int}
    memory::Int
    allocs::Int
end

function JLD.readas(p::ParametersPreV006)
    return Parameters(p.seconds, p.samples, p.evals, Float64(p.overhead), p.gctrial,
                      p.gcsample, p.time_tolerance, p.memory_tolerance)
end

function JLD.readas(t::TrialPreV006)
    new_times = convert(Vector{Float64}, t.times)
    new_gctimes = convert(Vector{Float64}, t.gctimes)
    return Trial(t.params, new_times, new_gctimes, t.memory, t.allocs)
end

function save(filename, args...)
    JLD.save(filename, VERSION_KEY, VERSIONS, args...)
    return nothing
end

@inline function load(filename, args...)
    # no version-based rules are needed for now, we just need
    # to check that version information exists in the file.
    if JLD.jldopen(file -> JLD.exists(file, VERSION_KEY), filename, "r")
        result = JLD.load(filename, args...)
    else
        JLD.translate("BenchmarkTools.Parameters", "BenchmarkTools.ParametersPreV006")
        JLD.translate("BenchmarkTools.Trial", "BenchmarkTools.TrialPreV006")
        result = JLD.load(filename, args...)
        JLD.translate("BenchmarkTools.Parameters", "BenchmarkTools.Parameters")
        JLD.translate("BenchmarkTools.Trial", "BenchmarkTools.Trial")
    end
    return result
end


text:
using Test, SpectralDistances, DSP, ControlSystems, Optim
import SpectralDistances: softmax
ip(x,y) = x'y/norm(x)/norm(y)
# const Continuous = SpectralDistances.Continuous

@testset "ISA" begin
    @info "Testing ISA"

    ## Extreme case 1- each measure is like a point, the bc should be close to the Euclidean mean of the center of each measure
    d = 2
    k = 4
    X0 = [1 1 2 2; 1 2 1 2]
    X = [X0 .+ 110rand(d) for _ in 1:4]
    S = ISA(X, iters=100, printerval=100)
    X̂ = [X[i][:,S[i]] for i in eachindex(X)]
    bc = mean(X̂)
    @test mean(bc) ≈ mean(mean, X)
    @test mean(bc, dims=2) ≈ mean(mean.(X, dims=2))
    if isinteractive()
        scatter(eachrow(reduce(hcat,X))...)
        scatter!(eachrow(bc)...)
    end

    ## Extreme case 2- variability between measures is small in comparison to the intra-measure variability. The bc should be close to the coordinate-wise Euclidean mean
    d = 2
    k = 4
    X0 = [1 1 2 2; 1 2 1 2]
    X = [X0 .+ 0.1rand(d) for _ in 1:4]
    S = ISA(X, iters=100, printerval=100)
    X̂ = [X[i][:,S[i]] for i in eachindex(X)]
    bc = mean(X̂)
    @test mean(bc) ≈ mean(mean, X)
    @test mean(bc, dims=2) ≈ mean(mean.(X, dims=2))
    if isinteractive()
        scatter(eachrow(reduce(hcat,X))...)
        scatter!(eachrow(bc)...)
    end

    # Extreme case 3, with weights where one weight is extremely large
    d = 2
    k = 4
    X0 = [1 1 2 2; 1 2 1 2]
    X = [X0 .+ 110rand(d) for _ in 1:4]
    w = s1([10000;ones(length(X)-1)])
    S = ISA(X, w, iters=100, printerval=100)
    X̂ = [w[i]*X[i][:,S[i]] for i in eachindex(X)]
    # bc = mean(X̂ .* w)
    bc = sum(X̂)
    @test mean(bc) ≈ mean(X[1]) rtol=0.1
    if isinteractive()
        scatter(eachrow(reduce(hcat,X))...)
        scatter!(eachrow(bc)..., alpha=0.2)
    end

end

@testset "ERD barycenters" begin
    @info "Testing ERD barycenters"

    d = EuclideanRootDistance(domain=Continuous(), weight=unitweight)
    models = [AR(ContinuousRoots([-1])), AR(ContinuousRoots([-3]))]
    bc = @inferred barycenter(d,models)
    @test bc isa AR
    @test roots.(Continuous(), bc)[1] ≈ -2

    models = [AR(ContinuousRoots([-1.0im,+im])), AR(ContinuousRoots([-3.0im,+3im]))]
    bc = @inferred barycenter(d,models)
    @test bc isa AR
    @test roots.(Continuous(), bc) ≈ [-2.0im,+2im]

    d = @inferred EuclideanRootDistance(domain=Continuous())
    models = [AR(ContinuousRoots([-1.0])), AR(ContinuousRoots([-3.0]))]
    bc = @inferred barycenter(d,models)
    @test bc isa AR
    @test roots.(Continuous(), bc)[1] ≈ -2

    models = [AR(ContinuousRoots([-im,+1.0im])), AR(ContinuousRoots([-3.0im,+3im]))]
    bc = barycenter(d,models)
    @test bc isa AR
    @test roots.(Continuous(), bc) ≈ [-2.0im,+2im]

    models = [AR(ContinuousRoots([-im,+im])), AR(ContinuousRoots([-3.0im,+2im]))]
    bc = barycenter(d,models)
    @test bc isa AR
    @test roots.(Continuous(), bc) ≈ [-2im,+1.5im]

    d = EuclideanRootDistance(domain=Continuous(), weight=residueweight)
    models = [AR(ContinuousRoots([-5.0-im,-1+im])), AR(ContinuousRoots([-1.0-im,-1+im]))]
    bc = barycenter(d,models)
    @test bc isa AR
    @test roots.(Continuous(), bc) ≈ [-1.1538-im,-1+im] rtol=0.01

end

##

@testset "alg1 alg2" begin
    @info "Testing alg1 alg2"

    Y = [[1. 1], [2. 2], [3. 3]] .+ Ref([0 0.001])
    X = [1.1 1.01]
    a = ones(2) ./ 2
    b = [ones(2) ./2 for _ in eachindex(Y)]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3, θ=0.5, printerval=190, γ=0.01, iters=50000, inneriters=50000, uniform=true) # TODO: there are type instabilities in this function but Cthulhu is not working at the time of writing.
    @test Xo ≈ [2 2] rtol=1e-2

    X = [0. 0.1]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3,θ=0.5, γ=0.01, inneriters=50000, uniform=true)
    @test Xo ≈ [2 2] rtol=1e-2

    X = [1. 3.]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/10.0, θ=0.5, γ=0.01, inneriters=50000, uniform=true)
    @test Xo ≈ [2 2] rtol=1e-2

    X = [0. 4.]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3, θ=0.5, γ=0.01, inneriters=50000, uniform=true)
    @test Xo ≈ [2 2] rtol=1e-2

    Y = [[1. 2], [2. 3], [3. 4]]
    X = [2.0 3]
    a = ones(2) |> s1
    b = [ones(2) ./2 for _ in eachindex(Y)]
    a1 = SpectralDistances.alg1(X,Y,a,b;β=1/3, tol=1e-5)
    @test a1[1] ≈ a1[2]

    X = [2.1 3]
    a1 = SpectralDistances.alg1(X,Y,a,b;β=1/3, tol=1e-3, printerval=1)
    @test a1[1] < a1[2]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3, θ=0.5, printerval=100, iters=500, tol=1e-3, γ=0.01, innertol=1e-3, inneriters=50000, uniform=true)
    @test Xo ≈ [2 3] rtol=5e-1
    # @test ao ≈ b[1] rtol=0.1

    X = [0. 0.1]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3, θ=0.5, printerval=100, innertol=1e-3, inneriters=50000, tol=1e-3, γ=0.1, uniform=true)
    @test Xo ≈ [2 3] rtol=5e-1

    X = [1. 3.]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3, θ=0.5, printerval=100, innertol=1e-4, inneriters=100, tol=1e-5, uniform=true)
    @test Xo ≈ [2 3] rtol=5e-1

    X = [0. 4.]
    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/2, θ=0.5, uniform=true)
    @test Xo ≈ [2 3] rtol=5e-1

    Y = [[1. 1], [2. 2], [3. 3]]
    X = [1.1 1.01]
    a = ones(2) |> s1
    b = [[0.2, 0.8] for _ in eachindex(Y)]
    @test SpectralDistances.alg1(Y[2],Y,b[1],b;β=1/3, printerval=1) ≈ b[1]

    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/3, θ=0.5, printerval=100)
    @test Xo ≈ [2 2] rtol=5e-1
    @test_broken ao ≈ b[1] rtol=1e-2

    ##
    Random.seed!(0)
    d = 2
    k = 4
    Y0 = 0.1*[1 1 2 2; 1 2 1 2]
    Y = [Y0[:,randperm(k)] .+ 1rand(d) for _ in 1:4]
    X = mean(Y) .+ Y0
    a = ones(k) |> s1
    b = [ones(k) |> s1 for _ in eachindex(Y)]
    a1 = SpectralDistances.alg1(X,Y,a,b;β=1/100.0, printerval=100)
    # @test a1 ≈ a rtol=0.01

    Xo,ao = SpectralDistances.alg2(X,Y,a,b;β=1/2, innertol=1e-5, tol=1e-6, printerval=20, inneriters=100000, solver=IPOT, γ=0.01, uniform=true)
    @test Xo ≈ mean(Y) rtol=1e-1
    @test ao ≈ a rtol = 0.1

    if isinteractive()
        scatter(eachrow(reduce(hcat,Y))...)
        scatter!(eachrow(X)..., alpha=0.8)
        scatter!(eachrow(Xo)..., alpha=0.4)
    end
    ##

    using JuMP, GLPK

    r6(x) = x#round.(x, digits=6)

    M = @inferred SpectralDistances.distmat_euclidean(X,Y[1])
    M2 = similar(M)
    @inferred SpectralDistances.distmat_euclidean!(M2,X,Y[1])
    @test M == M2

    g1,a1,b1 = ot_jump(M,a,b[1])
    g2,a2,b2 = @inferred sinkhorn_log!(M,a,b[1], β=0.0001, iters=1000000, tol=1e-7)
    g3,a3,b3 = @inferred IPOT(M,a,b[1], β=0.5, iters=100000)
    g4,a4,b4 = ot_convex(M,a,b[1])
    @test isapprox(ip(a1,a2), 1, atol=1e-1) || isapprox(ip(a1,a3), 1, atol=1e-1)
    @test isapprox(ip(b1,b2), 1, atol=1e-1) || isapprox(ip(b1,b3), 1, atol=1e-1)
    @test isapprox(ip(a1,a4), 1, atol=1e-1) || isapprox(ip(a2,a4), 1, atol=1e-1)
    @test isapprox(ip(b1,b4), 1, atol=1e-1) || isapprox(ip(b2,b4), 1, atol=1e-1)

    a = ones(k) |> s1
    b = [[1.0,2,3,4] |> s1 for _ in eachindex(Y)]
    M = SpectralDistances.distmat_euclidean(X,Y[1])
    g1,a1,b1 = ot_jump(M,a,b[1]) .|> r6
    g2,a2,b2 = sinkhorn_log(M,a,b[1], β=0.001, iters=50000, printerval=5000, tol=1e-9) .|> r6
    g3,a3,b3 = IPOT(M,a,b[1], β=0.5, iters=10000, printerval=5000, tol=1e-9) .|> r6
    @test isapprox(ip(a1,a2), 1, atol=1e-1) || isapprox(ip(a1,a3), 1, atol=1e-1)
    @test isapprox(ip(b1,b2), 1, atol=1e-1) || isapprox(ip(b1,b3), 1, atol=1e-1)

end

##

@testset "Barycentric coordinates" begin
    @info "Testing Barycentric coordinates"
    Random.seed!(0)
    d = 2
    k = 4
    S = 3
    X0 = [1 1 2 2; 1 2 1 2]

    res = map(1:10) do _
        X = [X0[:,randperm(k)] .+ 0.01.*randn.() .+ 10rand(d) for _ in 1:S]
        λ0 = randn(S) |> SpectralDistances.softmax
        p = [ones(k) |> s1 for _ in 1:S]
        q = ones(k) |> s1

        β = 0.1
        ql = barycenter(X, λ0, inneriters=200000, tol=1e-9, innertol=1e-8, β=β, solver=sinkhorn_log!)

        λh = barycentric_coordinates(X,ql,p,q, β=β, solver=sinkhorn_log!, robust=false, tol=1e-10)
        if isinteractive()
            scatter(eachrow(reduce(hcat,X))..., lab="X")
            scatter!(eachrow(ql)..., lab="initial bc")
        end
        costfun = λ -> SpectralDistances.sinkhorn_cost(X, ql, p, q, λ;
            solver = sinkhorn_log!,
            tol    = 1e-10,
            β      = β)

        isinteractive() && @show norm(λ0-λh)
    end
    @test median(res) < 0.01
    # mean(norm(s1(rand(4)) - s1(rand(4))) for _ in 1:10000)

    res = map(1:5) do _
        X = [X0[:,randperm(k)] .+ 10rand(d) for _ in 1:S]
        p = [randn(k) |> softmax for _ in 1:S]

        ql = X[1] .+ 0.01 .* randn.()
        q = p[1]
        λ0 = [1,0,0]

        β = 1/10.0
        λh = barycentric_coordinates(X,ql,p,q, β=β, solver=sinkhorn_log!, robust=false)
        if isinteractive()
            scatter(eachrow(reduce(hcat,X))..., lab="X")
            scatter!(eachrow(ql)..., lab="ql")
        end

        isinteractive() && @show norm(λ0-λh)
    end
    @test median(res) < 0.01

end

# using DoubleFloats
# df(v) = reduce(hcat, v)
# v = [[Double64(1) for _ in 1:2] for _ in 1:2]
# Test.@inferred df(v)

# v = randn(ComplexF64, 3)
# tv(v) = @views [real(v[1:2]); imag(v[1:2])]
# @inferred tv(v)

@testset "barycentric coordinates with models" begin
    @info "Testing barycentric coordinates with models"
    # The example below is likely to mix up the two lightly damped poles with euclidean root distance, making the bc poles end up inbetween the two clusters. The SRD should fix it
    Random.seed!(0)
    d = @inferred EuclideanRootDistance(domain=SpectralDistances.Continuous(),p=2)
    models = @inferred examplemodels(8)
    Xe = @inferred barycenter(d, models)
    # models = change_precision.(Float64, models)
    # Xe = change_precision(Float64, Xe)
    λ = barycentric_coordinates(d,models, Xe)
    @test λ ≈ s1(ones(length(λ))) atol=0.02

    λ = s1(rand(8))
    Xe = @inferred barycenter(d, models, λ)
    λ2 = barycentric_coordinates(d,models, Xe)
    @test ip(λ2,λ) > 0.85

    models = @inferred examplemodels(5) # When the number of models is smaller than the number of poles it's very exact
    λ = s1(rand(5))
    Xe = @inferred barycenter(d, models, λ)
    λ2 = barycentric_coordinates(d,models, Xe)
    @test λ2 ≈ λ atol=1e-6

    # Change weighting function
    d = @inferred EuclideanRootDistance(domain=SpectralDistances.Continuous(),p=2, weight=simplex_residueweight)
    models = @inferred examplemodels(8)
    Xe = barycenter(d, models)
    λ = barycentric_coordinates(d,models, Xe)
    @test λ ≈ s1(ones(length(λ))) atol=0.1

    λ = s1(rand(8))
    Xe = barycenter(d, models, λ)
    λ2 = barycentric_coordinates(d,models, Xe)
    @test ip(λ2,λ) > 0.85

    models = @inferred examplemodels(5) # When the number of models is smaller than the number of poles it's very exact
    λ = s1(rand(5))
    Xe = barycenter(d, models, λ)
    λ2 = barycentric_coordinates(d,models, Xe)
    @test λ2 ≈ λ atol=0.01

    models = examplemodels(8)
    λ = barycentric_coordinates(d,models, models[1], verbose=false, iters=5000, α0=20)
    isinteractive() && bar(λ)
    @test λ[1] > 0.5
    @test sum(λ) ≈ 1

    G = tf.(models)
    if isinteractive()
        plot()
        pzmap!.(G)
        pzmap!(tf(Xe), m=:c, title="Barycenter EuclideanRootDistance")
    end
    ##
    d = @inferred OptimalTransportRootDistance(domain=SpectralDistances.Continuous(),p=2, weight=residueweight, β=0.02)
    Xe = barycenter(d, models, solver=sinkhorn_log!)

    G = tf.(models)
    if isinteractive()
        plot()
        pzmap!.(G)
        pzmap!(tf(Xe), m=:c, title="Barycenter OptimalTransportRootDistance", lab="BC")
    end

    using Optim
    options = Optim.Options(
        store_trace = true,
        show_trace = false,
        show_every = 1,
        iterations = 50,
        allow_f_increases = true,
        time_limit = 100,
        x_tol = 1e-7,
        f_tol = 1e-7,
        g_tol = 1e-7,
        f_calls_limit = 0,
        g_calls_limit = 0,
    )

    method = LBFGS()
    method = ParticleSwarm()
    # d = OptimalTransportRootDistance(domain=SpectralDistances.Continuous(),p=2, weight=unitweight, β=0.1)
    Random.seed!(0)
    successrate = mean(1:5) do _
        λ = @inferred barycentric_coordinates(d,models,Xe, method, options=options, solver=sinkhorn_log!, robust=true, uniform=true, tol=1e-6)
        isinteractive() && bar(λ)
        median(λ) > 0.02
    end
    @test successrate >= 3/5

    G = tf.(models)
    if isinteractive()
        plot()
        pzmap!.(G, lab="")
        pzmap!(tf(Xe), m=:c, title="Barycenter OptimalTransportRootDistance", lab="BC")
        pzmap!(G[argmax(λ)], m=:c, lab="Largest bc coord", legend=true)
        # It's okay if the green dot do not match the blue exactly, there are limited models to choose from.
    end
end

@testset "kbarycenters" begin
    @info "Testing kbarycenters"

    using Clustering

    d = 2
    k = 4
    X0 = 0.1*[1 1 2 2; 1 2 1 2]
    X1 = 0.1*[1 2 3 4; 1 2 1 2]
    X00 = [X0[:,randperm(k)] .+ 0.01 .*randn.() .+ 1rand(d) for _ in 1:4]
    X11 = [X1[:,randperm(k)] .+ 0.01 .*randn.() .+ 1rand(d) .+ 2 for _ in 1:4]
    X = [X00; X11]

    p = [ones(k) |> s1 for _ in eachindex(X)]

    ##
    Q = SpectralDistances.kbarycenters(X,p,2, iters=7, solver=sinkhorn_log!, uniform=true, seed=:rand, tol=1e-3, innertol=1e-3, β=0.2, inneriters=50000, verbose=true)

    if isinteractive()
        scatter(eachrow(reduce(hcat,X))..., markerstrokewidth=false)
        scatter!(eachrow(reduce(hcat,Q.barycenters))..., markerstrokewidth=false, legend=false)
    end

    ##

    C = @inferred SpectralDistances.distmat_euclidean(X[1], X[5])
    c1 = @inferred SpectralDistances.kwcostfun(C,X[1], X[5],p[1],p[1],solver=IPOT, tol=1e-5, β=0.5, printerval=1, iters=200)
    c2 = @inferred SpectralDistances.kwcostfun(C,X[1], X[5],p[1],p[1],solver=sinkhorn_log, β=0.01)
    @test c1 ≈ c2 rtol=0.01

end

##

# d = 2
# k = 4
# S = 4
# X0 = [1 1 2 2; 1 2 1 2]
#
# using JuMP, GLPK
# ##
# X = [X0[:,randperm(k)] .+ 10rand(d) for _ in 1:S]
#
# λ0 = [0.9; 0.1ones(S-1)] |> s1
# ql = barycenter(X, λ0)
#
# ql2 = SpectralDistances.barycenter2(copy(X), λ0, printerval=10, γ=1, θ=0.5, iters=400, tol=1e-6, solver=SpectralDistances.sinkhorn_log)
#
# scatter(eachrow(reduce(hcat,X[2:end]))..., lab="X")
# scatter!(eachrow(X[1])..., lab="X1")
# # scatter!(eachrow(ql)..., lab="ql")
# scatter!(eachrow(ql2)..., lab="ql2")
## Measures with nonuniform weights, the bc should be pulled to atom 1 and 4 in the first measure. The trick seems to be to run ISA without weights and then use the weights to form the estimate

# d = 2
# k = 4
# X0 = [1 1 2 2; 1 2 1 2]
# X = [X0 .+ 0.3rand(d,k) .+ 0.020rand(d) for _ in 1:6]
# w = [ones(1,k) for _ in 1:length(X)]
# w[1][1] = 100
# w[1][4] = 100
# for i = 2:length(X)
#     w[i][1] = 0.01
#     w[i][4] = 0.01
# end
# w = s1.(w)
# W = reduce(vcat,w)
# W ./= sum(W,dims=1)
# w2 = [W[i,:]' for i in 1:length(X)]
# S = ISA(X, iters=1000, printerval=10)
# X̂ = [w2[i].*X[i][:,S[i]] for i in eachindex(X)]
# bc = sum(X̂)
# # @test mean(bc) ≈ mean(X[1]) rtol=0.1
# scatter(eachrow(reduce(hcat,X))...)
# scatter!([X[1][1,:]],[X[1][2,:]])
# scatter!(eachrow(bc)..., m=:square, legend=false, alpha=0.4)
# =============================================================================
## Tests for sinkhorn_diff, it does not seem to produce the correct gradient ==============
# C = [[mean(abs2, x1-x2) for x1 in eachcol(Xi), x2 in eachcol(ql)] for Xi in X]
# λl = SpectralDistances.softmax(1e-3randn(S))
#
# c,Pl,w = sinkhorn_diff(X,ql,p,q,C,λl, β=β, L=32, solver=sinkhorn_log)
# display(w)
#
# using ForwardDiff
# g = ForwardDiff.gradient(λl->SpectralDistances.sinkhorn_cost(X,ql,p,q,C,λl, β=β, L=32, solver=sinkhorn_log), λl)
# display(g)
# @test w ≈ g
## =================================================

# @test mean(getindex.(res,2)) >= 0.9
# ##
# X0 = [1. 1 2 2 3 3; 1 2 3 1 2 3]
# # X0 = [X0 X0]
# d,k = size(X0)
# S = 7
# X = [X0[:,randperm(k)] .+ 0.050randn(d) .+ 0.0001 .* randn(d,k) for _ in 1:S]
# w = ones(S) |> s1
# sw = ISA(X, w, iters=100, printerval=1)
# X̂ = SpectralDistances.barycentric_weighting(X,w,sw)
# # X̂2 = SpectralDistances.barycentric_weighting2(X,w,sw)
# scatter(eachrow(reduce(hcat,X))...)
# scatter!(eachrow(X̂)..., alpha=0.5)
# # scatter!(eachrow(X̂2)..., alpha=0.5)
# # In the plot above, the bc should have the same shape as the acnhors
#
#
# error("keep track of the cost function noted in the paper and ensure that it decreases")
# ##
#
# # Now testing for models
# ζ = [0.1, 0.3, 0.7]
# models = map(1:50) do _
#     pol = [1]
#     for i = eachindex(ζ)
#         pol = SpectralDistances.polyconv(pol, [1,2ζ[i] + 0.1randn(),1+0.1randn()])
#     end
#     AR(Continuous(),pol)
# end
#
# # λ0 = rand(length(models)) |> s1
# # qmodel = barycenter(distance, models, λ0)
#
# distance = OptimalTransportRootDistance(domain=Continuous(),p=2, weight=unitweight)
#
# # We choose the point to be projected equal to one of the anchor points. In this case, the barycentric coordinates is a one-hot vector
# qmodel = models[1]
# q_proj, λ = barycentric_coordinates(distance, models, qmodel, γ = 0.2, robust=true)
# @test λ[1] ≈ 1 atol=0.2
#
# qmodel = models[end]
# q_proj, λ = barycentric_coordinates(distance, models, qmodel, γ=0.2, robust=false)
# @test λ[end] ≈ 1 atol=0.2
#
#
# qbc = barycenter(distance, models, λ) # TODO: the problem might be in this method of barycenter, not sure if it's actually any good
#
#
#
# plot()
# pzmap!.(tf.(models), lab="")
# pzmap!(tf(qmodel), m=:c, lab="q")
# pzmap!(tf(q_proj), m=:c, lab="q_proj", legend=true)


text:
facts("Rendering.") do

    output = IOBuffer()

    context("Query output.") do

        writemime(output, MIME"text/plain"(), run(@query("")))
        writemime(output, MIME"text/plain"(), run(@query("",1)))
        writemime(output, MIME"text/plain"(), run(@query("",2)))
        writemime(output, MIME"text/plain"(), run(@query("",3)))
        writemime(output, MIME"text/plain"(), run(@query("",4)))
        writemime(output, MIME"text/plain"(), run(@query("",5)))

        writemime(output, MIME"text/plain"(), run(@query(@document())))
        writemime(output, MIME"text/plain"(), run(@query(Docile.@document())))

        writemime(output, MIME"text/plain"(), run(@query(Docile.meta())))

        writemime(output, MIME"text/plain"(), run(@query(Docile.meta)))

        writemime(output, MIME"text/plain"(), run(@query(Docile.Interface.Entry, 2)))

        writemime(output, MIME"text/plain"(), run(@query(Docile)))

    end

    context("Doctests output.") do

        writemime(output, MIME"text/plain"(), doctest(Docile))
        writemime(output, MIME"text/plain"(), doctest(Docile.Interface))
        writemime(output, MIME"text/plain"(), doctest(Lexicon))

        writemime(output, MIME"text/plain"(), passed(doctest(Docile)))
        writemime(output, MIME"text/plain"(), failed(doctest(Docile.Interface)))
        writemime(output, MIME"text/plain"(), skipped(doctest(Lexicon)))

    end

    context("Saving static content.") do
        index  = Index()
        config = Config(include_internal = true)
        for modname in (Lexicon, Docile, Docile.Interface), ft in ("md", "html")
            dir = joinpath(tempdir(), randstring())
            f = joinpath(dir, "$(modname).$(ft)")
            # TODO: update / save index is not implemented for html
            if ft == "md"
                update!(index, save(f, modname, config; mathjax = true))
                save(joinpath(dir, "index_$(modname).$(ft)"), index, config)
            else
                save(f, modname, config; mathjax = true)
            end
            rm(dir, recursive = true)
        end
    end

    if VERSION < v"0.4.0-dev+4393"
        context("Testing relpath.") do
            sep = Base.path_separator
            filepaths = [
                "$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "$(sep)home$(sep)user$(sep)dir_withendsep$(sep)",
                "$(sep)home$(sep)dir2_withendsep$(sep)",
                "$(sep)home$(sep)test.md",
                "$(sep)home",
                # Special cases
                "$(sep)",
                "$(sep)home$(sep)$(sep)$(sep)"
                ]

            startpaths = [
                "$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)genindex.md",
                "$(sep)multi_docs$(sep)genindex.md",
                "$(sep)home$(sep)user$(sep)dir_withendsep$(sep)",
                "$(sep)home$(sep)dir2_withendsep$(sep)",
                "$(sep)home$(sep)test.md",
                "$(sep)home",
                # Special cases
                "$(sep)",
                "$(sep)home$(sep)$(sep)$(sep)"
                ]

            # generated with python's relpath
            relpath_expected_results = [
                "..$(sep)Lexicon.md",
                "..$(sep)..$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "..$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "..$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "..$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)Lexicon.md",
                "..$(sep)lib$(sep)file1.md",
                "..$(sep)..$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "home$(sep)user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "user$(sep).julia$(sep)v0.4$(sep)Lexicon$(sep)docs$(sep)api$(sep)lib$(sep)file1.md",
                "..$(sep)..$(sep)..$(sep)..$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "..$(sep)..$(sep)home$(sep)user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "..$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "..$(sep)user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "..$(sep)user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "home$(sep)user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "user$(sep).julia$(sep)v0.4$(sep)Docile$(sep)docs$(sep)api$(sep)Docile.md",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)dir_withendsep",
                "..$(sep)..$(sep)home$(sep)user$(sep)dir_withendsep", ".", "..$(sep)user$(sep)dir_withendsep",
                "..$(sep)user$(sep)dir_withendsep", "user$(sep)dir_withendsep",
                "home$(sep)user$(sep)dir_withendsep", "user$(sep)dir_withendsep",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)dir2_withendsep",
                "..$(sep)..$(sep)home$(sep)dir2_withendsep", "..$(sep)..$(sep)dir2_withendsep", ".",
                "..$(sep)dir2_withendsep", "dir2_withendsep", "home$(sep)dir2_withendsep", "dir2_withendsep",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)test.md",
                "..$(sep)..$(sep)home$(sep)test.md", "..$(sep)..$(sep)test.md", "..$(sep)test.md", ".",
                "test.md", "home$(sep)test.md", "test.md", "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..",
                "..$(sep)..$(sep)home", "..$(sep)..", "..", "..", ".", "home", ".",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..", "..$(sep)..",
                "..$(sep)..$(sep)..", "..$(sep)..", "..$(sep)..", "..", ".", "..",
                "..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..$(sep)..", "..$(sep)..$(sep)home",
                "..$(sep)..", "..", "..", ".", "home", "."
                ]

            idx = 0
            for filep in filepaths
                for startp in startpaths
                    res = Lexicon.relpath(filep, startp)
                    idx += 1
                    @fact res --> relpath_expected_results[idx] "Excpected: $(relpath_expected_results[idx])"
                end
            end

            # Additional cases
            @fact_throws ArgumentError Lexicon.relpath("$(sep)home$(sep)user$(sep)dir_withendsep$(sep)", "")
            @fact_throws ArgumentError Lexicon.relpath("", "$(sep)home$(sep)user$(sep)dir_withendsep$(sep)")

        end
    end

end


text:
using VoronoiCells
using VoronoiDelaunay
using Test

@testset "Clipping" begin
    rect = VoronoiCells.BoundingBox

    @testset "Lines wholly within the rectangle are not modified" begin
        A = VoronoiDelaunay.Point2D(1.3, 1.5)
        B = VoronoiDelaunay.Point2D(1.7, 1.5)
        l = VoronoiDelaunay.Line(A, B)

        @test l == VoronoiCells.clip(l, rect)
    end

    @testset "Clip line with one endpoint outside rectangle" begin
        A = VoronoiDelaunay.Point2D(1.3, 1.5)
        B = VoronoiDelaunay.Point2D(VoronoiCells.right(rect) + 1, 1.5)
        l = VoronoiDelaunay.Line(A, B)

        expected_B = VoronoiDelaunay.Point2D(VoronoiCells.right(rect), gety(B))
        expected_line = VoronoiDelaunay.Line(A, expected_B)

        @test VoronoiCells.clip(l, rect) == expected_line
    end

    @testset "Clip line with both endpoints outside rectangle" begin
        A = VoronoiDelaunay.Point2D(VoronoiCells.left(rect) - 1, 1.5)
        B = VoronoiDelaunay.Point2D(VoronoiCells.right(rect) + 1, 1.5)
        l = VoronoiDelaunay.Line(A, B)

        expected_A = VoronoiDelaunay.Point2D(VoronoiCells.left(rect), gety(A))
        expected_B = VoronoiDelaunay.Point2D(VoronoiCells.right(rect), gety(B))
        expected_line = VoronoiDelaunay.Line(expected_A, expected_B)

        @test VoronoiCells.clip(l, rect) == expected_line
    end

    @testset "Line does not intersect rectangle" begin
        l = VoronoiDelaunay.Line(
            VoronoiDelaunay.Point2D(1.3, VoronoiCells.lower(rect) - 1),
            VoronoiDelaunay.Point2D(1.7, VoronoiCells.lower(rect) - 1)
        )

        @test isnothing(VoronoiCells.clip(l, rect))
    end

    @testset "Clip line whose intersection with rectangle is one point on boundary" begin
        A = VoronoiDelaunay.Point2D(VoronoiCells.right(rect), 1.5)
        B = VoronoiDelaunay.Point2D(VoronoiCells.right(rect) + 1, 1.5)
        l = VoronoiDelaunay.Line(A, B)

        expected_B = VoronoiDelaunay.Point2D(VoronoiCells.right(rect), gety(B))
        expected_line = VoronoiDelaunay.Line(A, expected_B)

        @test VoronoiCells.clip(l, rect) == expected_line
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    API implementation for two body orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

get_epoch(orbp::OrbitPropagatorTwoBody) = orbp.tbd.epoch

function get_mean_elements(orbp::OrbitPropagatorTwoBody)
    tbd = orbp.tbd
    orb = KeplerianElements(tbd.epoch + tbd.Δt/86400,
                            tbd.a_0,
                            tbd.e_0,
                            tbd.i_0,
                            tbd.Ω_0,
                            tbd.ω_0,
                            tbd.f_k)

    return orb
end

"""
    init_orbit_propagator(::Val{:twobody}, epoch::Number, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number; μ::T = m0) where T

Initialize the two body orbit propagator.

# Args

* `epoch`: Initial orbit epoch [Julian Day].
* `a_0`: Initial mean semi-major axis [m].
* `e_0`: Initial mean eccentricity.
* `i_0`: Initial mean inclination [rad].
* `Ω_0`: Initial mean right ascension of the ascending node [rad].
* `ω_0`: Initial mean argument of perigee [rad].
* `f_0`: Initial mean true anomaly [rad].
* `orb_0`: Instance of the structure `KeplerianElements` with the initial mean
           orbital elements [SI].

## Keywords

* `μ`: (OPTIONAL) Standard gravitational parameter of the central body
       \\[m^3/s^2] (**Default** = `m0`).

"""
function init_orbit_propagator(::Val{:twobody}, epoch::Number,
                               a_0::Number, e_0::Number, i_0::Number,
                               Ω_0::Number, ω_0::Number, f_0::Number;
                               μ::T = m0) where T

    # Create the new Two Body propagator structure.
    tbd = twobody_init(epoch, a_0, e_0, i_0, Ω_0, ω_0, f_0; μ = μ)

    # Create and return the orbit propagator structure.
    return OrbitPropagatorTwoBody(tbd)
end

function init_orbit_propagator(::Val{:twobody}, orb_0::Orbit; μ::T = m0) where T
    # Convert the orbit representation to Keplerian elements.
    k_0 = convert(KeplerianElements, orb_0)

    return init_orbit_propagator(Val(:twobody), k_0.t, k_0.a, k_0.e, k_0.i,
                                 k_0.Ω, k_0.ω, k_0.f; μ = μ)
end

function propagate!(orbp::OrbitPropagatorTwoBody{T}, t::Number) where T
    # Auxiliary variables.
    tbd = orbp.tbd

    # Propagate the orbit.
    r_i, v_i = twobody!(tbd, t)

    return r_i, v_i
end

function step!(orbp::OrbitPropagatorTwoBody, Δt::Number)
    # Auxiliary variables.
    tbd = orbp.tbd

    # Propagate the orbit.
    r_i, v_i = twobody!(tbd, tbd.Δt + Δt)

    # Return the information about the step.
    return r_i, v_i
end


text:
pages = Any[
        "index.md",
        "type_hierarchy.md",
        "Usage" => [
            "numerical_integration.md",
            "quadrature_rules.md",
            "orthogonal_polynomials_canonical.md",
            "gaussian_mixture_model.md",
            "multiple_discretization.md",
            "scalar_products.md",
            "Polynomial Chaos" => [ "Basic Usage" => "pce_tutorial.md",
                                    "Chi Squared, One DOF" => "chi_squared_k1.md",
                                    "Chi Squared, Several DOFs" => "chi_squared_k_greater1.md",
                                    "Random ODE" => "random_ode.md"
                                  ],
            "Optimal Power Flow" => "DCsOPF.md"
            ],
    "math.md",
    "functions.md"
    ]

max_stars_repo_name: timueh/PolyChaos

max_issues_repo_name: timueh/PolyChaos

max_forks_repo_name: timueh/PolyChaos


text:
"""
    AbstractPopulation

An abstract type for population-based uncertain values.
"""
abstract type AbstractPopulation <: AbstractUncertainValue end

export AbstractPopulation


text:
using Base.Test
using openBF

function runTest(test_folder)

  cd(test_folder)
  try
    run(`cp ../../main.jl main.jl`)
    run(`$JULIA_HOME/julia main.jl $test_folder`)
    cd("..")
    return true

  catch
    cd("..")
    return false
  end
end

# Run tests

tic()
println("\nTest 1 - Single artery")
@time @test runTest("single-artery")

println("\nTest 2 - Bifurcation")
@time @test runTest("bifurcation")

println("\nTest 3 - External pressure")
@time @test runTest("external-pressure")

println("\nTest 4 - Anastomosis")
@time @test runTest("anastomosis")

println("\nTest 5 - Conjunction")
@time @test runTest("conjunction")
toc()


text:
module ExploreCMEO

import BinDeps #To download database from internet.
import HDF5
import HDF5: HDF5Group, HDF5Dataset

import Gtk
const _Gtk = Gtk.ShortNames

const URL_PUBLISHEDDATA = "https://raw.githubusercontent.com/alainman-krh/ExploreCMEODonnees.jl/master/ExploreCMEO.h5"
const PATH_DB = Ref(abspath("./ExploreCMEO.h5")) #In cwd

include("base.jl")
include("db_access.jl")
include("dbaccess_fieldview.jl")
include("gtk_ext.jl")
include("gtk_base.jl")
include("gtk_editdlg.jl")
include("gtk_top.jl")

#==Initialization
===============================================================================#
function __init__()
#	PATH_DB[] = abspath(joinpath(dirname(pathof(ExploreCMEODonnees)), "..", "ExploreCMEO.h5"))
	dbpath = PATH_DB[]
	if !isfile(dbpath)
		download_database(dbpath)
	end
end

end # module


text:
immutable Relation{S,T} <: AbstractVector{T}
    src::S
end

# map from tuples of attributes (as Symbols) to respective Attributed leaf type
let type_cache = Dict{Tuple{Vararg{Symbol}}, DataType}()
    global Relation
    function Relation{T<:Attributed}(::Type{T}, cols::AbstractVector...)::Relation
        A = type_cache[tuple(fieldnames(T)...)] # This key has to exist
        S = A{map(typeof, cols)...}
        return Relation{S,T}(S(cols...))
    end

    function Relation(; kwargs...)
        n = length(kwargs)
        attrs, cols = Vector{Symbol}(n), Vector{Any}(n)
        for (j,kwarg) in enumerate(kwargs)
            attrs[j]=kwarg[1]; cols[j]=kwarg[2];
        end
        return Relation(tuple(attrs...), cols...)
    end

    # function Relation(cols::Union{Vector{AbstractVector}, Tuple{Vararg{AbstractVector}}},
    function Relation(attrs::Tuple{Vararg{Symbol}}, cols::AbstractVector...)::Relation
        length(attrs) == length(cols) || throw(ArgumentError("Attributes, columns must be same length."))
        # We use @get! so that we only run eval if type_cache has no key attrs
        A = Base.@get! type_cache attrs begin
            _T = gensym()
            type_name_ex = Expr(:curly, _T)
            type_decl_ex = Expr(:(<:), type_name_ex, :Attributed)
            type_fields_ex = Expr(:block)
            for attr in attrs
                param = gensym()
                push!(type_name_ex.args, param)
                push!(type_fields_ex.args, Expr(:(::), attr, param))
            end
            typedef_ex = Expr(:type, false, type_decl_ex, type_fields_ex)
            # NOTE: uses of eval!
            eval(Relations, typedef_ex)
            eval(Relations, _T)
        end
        S, T = A{map(typeof, cols)...}, A{map(eltype, cols)...}
        return Relation{S,T}(S(cols...))
    end
end

Relation{T}(a::Attributed, ::Type{T}) = Relation{typeof(a),T}(a)

# as an iterator
Base.eltype{S,T}(r::Relation{S,T}) = T
Base.ndims{A<:Attributed}(r::Relation{A}) = 1

# as a container
Base.empty!(r::Relation) = (foreach(empty!, r.src); r)
Base.similar{S<:Attributed,T}(r::Relation{S,T}) = empty!(Relation(T, map(similar, r.src)...))
Base.similar{S<:Attributed,T}(r::Relation{S,T}, n::Integer) = Relation(T, map(c->similar(c,n), r.src)...)
Base.copy{S<:Attributed,T}(r::Relation{S,T}) = Relation{S,T}(S(map(copy, r.src)...))

# as an array
Base.size{A<:Attributed}(r::Relation{A}) = (length(r),)

# as a vector
Base.getindex{A<:Attributed,T}(r::Relation{A,T}, i) = T(ith_all(i, r.src)...)
Base.getindex{S<:Attributed,T}(r::Relation{S,T}, I::Vector) =
    Relation(T, map(c->getindex(c,I), r.src)...)
Base.push!{A<:Attributed}(r::Relation{A}, v) = (foreach(push!, r.src, v))
Base.linearindexing{A<:Attributed}(r::Relation{A}) = Base.LinearFast()
# NOTE: not all sources will support length
Base.length{A<:Attributed}(r::Relation{A}) = length(first(r.src))

# as a bundle of vectors (NOTE: not all sources will support this)
columns{S<:Attributed}(r::Relation{S}) = r.src

# as a relation
function attributes end
function index end
attributes{S,T}(r::Relation{S,T}) = fieldnames(T)
function index{S,T}(r::Relation{S,T})::Dict{Symbol, Int}
    index = Dict{Symbol, Int}()
    for (j, attr) in enumerate(fieldnames(T))
        index[attr] = j
    end
    return index
end
project{S<:Attributed}(r::Relation{S}, attrs::Symbol...) =
    Relation(attrs, [ getfield(r.src, attr) for attr in attrs ]...)

function together{S1<:Attributed,S2<:Attributed}(r1::Relation{S1}, r2::Relation{S2})
    attrs1, attrs2 = attributes(r1), attributes(r2)
    attrs = vcat(attrs1, attrs2)
    cols = vcat([ getfield(r1.src, attr) for attr in attrs1 ],
                [ getfield(r2.src, attr) for attr in attrs2 ])
    return Relation(tuple(attrs...), cols...)
end

# as queryable
SQ.prepare{S<:Attributed}(r::Relation{S}) = r
SQ.prepare{S<:Relation}(gr::Grouped{S}) = gr
SQ._with{S<:Attributed}(q::SQ.Node, ::Tuple{Relation{S}}) = collect(q)
SQ._with{S<:Attributed, T<:Attributed}(q::SQ.Node, ::Tuple{Relation{S}, Relation{T}}) = collect(q)

SQ.as{S<:Attributed}(r::Relation{S}, ::Type{Relation}) = r
SQ.as{S<:Attributed}(g::Grouped{Relation{S}}, ::Type{Relation}) = g


text:
__precompile__()
# ----------------------------------------------------------------------------
module FAD # Forward Automatic Differentiation demo
#=
    see also:
        - the walkthrough at http://vladium.com/tutorials/study_julia_with_me/multiple_dispatch/

    working with this module from REPL:

    julia>include("<path>/FAD.jl")
    julia>using Main.FAD
    ... REPL can now use derivative(), ∂(), root_solve(), etc ...
=#
# ............................................................................

import Base; # needed in order to add base ops overloads

export derivative, ∂, myerf, root_solve; # for convenience of 'using Main.FAD' in REPL

# ............................................................................

struct Context <: Number
    v   ::Number
    ∂   ::Number
end
# ............................................................................
# binary ops:

# note: ':()' is necessary to reference certain functions with single-character names

Base.:(+)(lhs ::Context, rhs ::Context) = Context(lhs.v + rhs.v, lhs.∂ + rhs.∂)
Base.:(-)(lhs ::Context, rhs ::Context) = Context(lhs.v - rhs.v, lhs.∂ - rhs.∂)

Base.:(*)(lhs ::Context, rhs ::Context) = Context(lhs.v * rhs.v, lhs.v * rhs.∂ + lhs.∂ * rhs.v)
Base.:(/)(lhs ::Context, rhs ::Context) = Context(lhs.v / rhs.v, (lhs.∂ * rhs.v - lhs.v * rhs.∂) / rhs.v^2)

# unary ops:

Base.:(+)(x ::Context) = x
Base.:(-)(x ::Context) = Context(- x.v, - x.∂)

# math ops:

Base.sin(d ::Context) = Context(sin(d.v),   cos(d.v) * d.∂)
Base.cos(d ::Context) = Context(cos(d.v), - sin(d.v) * d.∂)
# ... add more as needed for your own functions ...

# ............................................................................
# conversion/promotion:

Base.promote_rule(::Type{Context}, ::Type{<: Number}) = Context
Base.convert(::Type{Context}, x ::Real) = Context(x, 0.0)

# ............................................................................
"""
    Given a function 'f' of one variable, return its derivative with respect to that variable.
"""
function derivative(f ::Function)
    return x ::Number -> f(Context(x, 1.0)).∂ # discard value, return derivative
end

∂ = derivative # add a nice-looking alias

# ............................................................................

# to get "official" erf(x), do
# pkg> add SpecialFunctions
# julia> using SpecialFunctions

function myerf(x ::Number)
    Σ = 0.0
    x² = x * x
    for k in 0 : 20 # hardcoding the number of summation terms for simplicity
        Σ += x / (factorial(k) * (2k + 1))
        x *= -x²
    end
    return 2.0 / √π * Σ
end
# ............................................................................
"""
    Newton's method for finding root(s) of 'f'
"""
function root_solve(f ::Function, x₀ ::Number; ϵ = 1e-8)
    i = 1
    while true
        ctx = f(Context(x₀, 1.0))
        println("[$i]: f($x₀)\t= $(ctx.v)")
        abs(ctx.v) < ϵ && break
        x₀ -= ctx.v / ctx.∂ # ignoring issues with division by small floats, etc for simplicity
        i += 1
    end
    return x₀
end

end # of module
# ----------------------------------------------------------------------------


text:
import PlanningVsCompensatoryCaching: load_data, sort_semantically,
datatotextable, plotallindividualdata, logcategoricalposterior, plotcomparison,
pgfsave, compute_posterior

@isdefined(RESULT_DIR) || (RESULT_DIR = joinpath(@__DIR__, "results/"))

@info "Analysing Experiment 1..."

data1, data1_sorted, cached1 = load_data(joinpath(@__DIR__, "data", "experiment1.yaml"))

# Export raw and sorted data.
write(joinpath(RESULT_DIR, "datatable.tex"), datatotextable(data1, data1_sorted))

# Define qualitatively some hypotheses and plot them together with the data.
data_sorted_with_hypothesis = deepcopy(data1_sorted)
data_sorted_with_hypothesis[:FPH1] = (comp_order = [:A, :B, :C],
                                      prefed_order = [:P, :M],
                                      cached = [5 1; 1 1; 1 1])
data_sorted_with_hypothesis[:FPH2] = (comp_order = [:A, :B, :C],
                                      prefed_order = [:P, :M],
                                      cached = [3 1; 1 3; 3 2])
data_sorted_with_hypothesis[:CCH] = (comp_order = [:A, :B, :C],
                                     prefed_order = [:P, :M],
                                     cached = [2 3; 4 1; 1 3])
plotallindividualdata(data_sorted_with_hypothesis, dir = RESULT_DIR,
                      xlabelat = "FPH1", titles = Dict(:FPH1 => "FPH 1",
                                                       :FPH2 => "FPH 2"));

models = [[(bird_indep = b, compartement_indep = c, food_indep = f, hyp = nothing)
           for b in (true, false), f in (true, false), c in (true, false) if !(b && c && f)]...;
          [(bird_indep = b, compartement_indep = false, food_indep = false, hyp = h)
           for b in (true, false), h in (:cch, :fph1, :fph2)]...]

# Compute Posterior for each condition.
posteriors = compute_posterior(cached1, models)
@show posteriors

# Generate comparison plots.
p1 = plotcomparison(posteriors, models, title = "\$P(M|D)\$", ymax = 5);
pgfsave(joinpath(RESULT_DIR, "posteriors.tex"), p1, include_preamble = false)

"""
Conclusion

1 bird clear FPH 1 (Lisbon)
2 birds unclear (Jerusalem, Quito)
3 birds rather CCH (Caracas, Washington, Wellington)

BUT: I don't want to overinterpret these results since the data is quite well
explained with a model that distributes uniformly among compartments
"""

#####
##### Data Analysis: Experiment 2
#####

import PlanningVsCompensatoryCaching: load_data, sort_semantically,
datatotextable2, plotallindividualdata, logcategoricalposterior, plotcomparison,
pgfsave, compute_posterior

@info "Analysing Experiment 2..."

data2, data2_sorted, cached2 = load_data(joinpath(@__DIR__, "data", "experiment2.yaml"))

# Export raw and sorted data.
write(joinpath(RESULT_DIR, "datatable2.tex"), datatotextable2(data2, data2_sorted))

# Define qualitatively some hypotheses and plot them together with the data.
data_sorted_with_hypothesis = deepcopy(data2_sorted)
data_sorted_with_hypothesis[:FPH1] = (comp_order = [:A, :B, :C],
                                      prefed_order = [:F, :N],
                                      cached = [4; 1; 1])
data_sorted_with_hypothesis[:FPH2] = (comp_order = [:A, :B, :C],
                                      prefed_order = [:F, :N],
                                      cached = [2; 1; 2])
data_sorted_with_hypothesis[:CCH] = (comp_order = [:A, :B, :C],
                                     prefed_order = [:F, :N],
                                     cached = [1; 2; 1])
data_sorted_with_hypothesis[:FPH1NF] = (comp_order = [:A, :B, :C],
                                        prefed_order = [:N, :F],
                                        cached = [1; 3; 1])
data_sorted_with_hypothesis[:FPH2NF] = (comp_order = [:A, :B, :C],
                                        prefed_order = [:N, :F],
                                        cached = [1; 3; 1])
data_sorted_with_hypothesis[:CCHNF] = (comp_order = [:A, :B, :C],
                                       prefed_order = [:N, :F],
                                       cached = [2; 1; 2])
plotallindividualdata(data_sorted_with_hypothesis, exp = 2, xlabelat = "FPH1",
                      y_axis_line_style = "{->}", dir = RESULT_DIR, ylabel = "cached",
                      titles = Dict(:FPH1 => "FPH 1", :FPH2 => "FPH 2",
                                    :CCHNF => "CCH", :FPH1NF => "FPH 1",
                                    :FPH2NF => "FPH 2"));

models = [[(bird_indep = b, compartement_indep = c, food_indep = true, hyp = nothing)
           for b in (true, false), c in (true, false) if !(b && c)]...;
          [(bird_indep = b, compartement_indep = false, food_indep = true, hyp = h)
           for b in (true, false), h in (:cch, :fph1, :fph2)]...]

# Compute Posterior for each condition and 40 caching trials.
idx_fn = findall(x -> x ∈ (:Caracas, :Washington, :Quito), collect(keys(data2)))
idx_nf = findall(x -> x ∈ (:Wellington, :Rome, :Lisbon), collect(keys(data2)))
lp_fn = compute_posterior(cached2[idx_fn], models, group = :FN, unnormalizedlog = true)
lp_nf = compute_posterior(cached2[idx_nf], models, group = :NF, unnormalizedlog = true)
posteriors = (x -> x ./ sum(x))(exp.(lp_fn .+ lp_nf))

@show posteriors

p1 = plotcomparison(posteriors, models, title = "\$P(M|D)\$", altpattern = "north east lines", ymax = 5, exp2 = true, legendtoname = "modelcomparisonlegend2", legendentries = ["bird independent", "bird dependent"]);
pgfsave(joinpath(RESULT_DIR, "posteriors2.tex"), p1, include_preamble = false)

"""
Conclusion

3 birds rather FPH 1 (Caracas, Washington, Quito)
3 birds rather CCH (Rome, Lisbon Wellington)

BUT: I don't want to overinterpret these results since the data is quite well
explained with a model that distributes uniformly among compartments
"""


text:
# load ImpvolEquilibrium first so that methods are accessible
@everywhere include("../../../equilibrium.jl")
@everywhere using ImpvolEquilibrium, Logging
@everywhere Logging.configure(level=INFO)

@everywhere using FileIO, JLD2
@everywhere parameters = load("../common_parameters.jld2")["parameters"]

# parameters that govern counterfactual
@everywhere include("change_parameters.jl")

@time results = pmap(t -> (t, ImpvolEquilibrium.period_wrapper(parameters, t)), 1:parameters[:T])
include("../../../save.jl")


text:
include("src/ssp/windprofile.jl")
using .windprofile
using Plots
#using BenchmarkTools

z0 = 0.03; α = 0.21; vref = 2.7; href = 10
h = range(0,1000,step = 1.0)

# log Linear - law model
vh1 = similar(h) # pre-allocation
for i=2:length(h)
    vh1[i] = windprofile_logLL(h[i],z0,vref,href)
end
plot(vh1,h)

# logarithmic law model
vh2 = similar(h) # pre-allocation
for i=2:length(h)
    vh2[i] = windprofile_logL(h[i],z0,vref,href)
end
plot!(vh2,h)

# DH log model
vh3 = similar(h) # pre-allocation
for i=2:length(h)
    vh3[i] = windprofile_DH(h[i],z0,vref,href)
end
plot!(vh3,h)

# power law model
vh4 = similar(h) # pre-allocation
for i=2:length(h)
    vh4[i] = windprofile_PL_α(h[i],1/7,vref,href)
end
plot!(vh4,h)

# power law model mean anual
vh5 = similar(h) # pre-allocation
for i=2:length(h)
    vh5[i] = windprofile_PL_α(h[i],α,vref,href)
end
plot!(vh5,h)

# power law model Counihan model
vh6 = similar(h) # pre-allocation
αC = 0.24 + 0.096*log(z0) + 0.016*(log(z0))^2
for i=2:length(h)
    vh6[i] = windprofile_PL_α(h[i],αC,vref,href)
end
plot!(vh6,h)

# power law SH model
vh7 = similar(h) # pre-allocation
αSH = 1.03 + 0.31*log(z0) + 0.03*(log(z0))^2
for i=2:length(h)
    vh7[i] = windprofile_PL_α(h[i],αSH,vref,href)
end
plot!(vh7,h)

# power law αP model
vh8 = similar(h) # pre-allocation
for i=2:length(h)
    vh8[i] = windprofile_PL_α_P(h[i],z0,vref,href)
end
plot!(vh8,h)


text:
## MFCC.jl
## (c) 2013-2015 David A. van Leeuwen, (c) 2005--2012 Dan Ellis

## Recoded from / inspired by melfcc from Dan Ellis's rastamat package.

module MFCC

export powspec, audspec, postaud, lifter
export mfcc, deltas, warp, sdc, znorm, znorm!, stmvn
export feasave, feaload, feasize
export feacalc

using DSP
using HDF5
using WAV
using SpecialFunctions ## erfinv
using Statistics

include("rasta.jl")
include("mfccs.jl")
include("feacalc.jl")
include("io.jl")

end


text:
using Test, BitBasis

@testset "bit literal" begin
    @test bit"1011" == 0b1011
    @test bit"100_111" == 0b100_111
    @test bit"10_100_11" == bit"1010011"
    @test bit"10_100_11" != bit"01010011"
    @test bit"1011" === bit_literal(1, 1, 0, 1)
    @test BitBasis.parse_bit(BigInt, "10101010101010101010101010101010010101010101"^10) isa LongBitStr
    @test_throws ErrorException BitBasis.parse_bit(
        Int64,
        "10101010101010101010101010101010010101010101"^10,
    )
end

@testset "operations" begin
    @test (bit"1010" + bit"0101") == bit"1111" == 15
    @test (bit"1111" - bit"0101") == bit"1010" == 10
    @test (bit"01010" * 2) == (2 * bit"01010") == bit"10100" == 20
    @test (bit"1010" ÷ 2) == bit"0101" == 5
    @test bcat(bit"101", bit"100", bit"111") == bit"101100111"
    @test bcat(bit"101" for k in 1:3) == bit"101101101"
    @test (bit"00101" << 2) == bit"10100"
    @test (bit"1101" >> 2) == bit"0011"
    @test bit"10011" == bit"10011"

    @test collect(bit"1101101") == Int64[1, 0, 1, 1, 0, 1, 1]
    @test length(bit"000":bit"111") == 8
    v = zeros(8)
    v[buffer(bit"101")+1] = 1
    @test onehot(bit"101") == v
    @test onehot(bit"101", 2) == [v v]

    @test repeat(bit"101", 3) == bit"101101101"
end

@testset "conversions" begin
    for T in [Int8, Int16, Int32, Int64, Int128, BigInt]
        res = T(bit"10")
        @test typeof(res) == T
        @test res == 2
        res = convert(T, bit"10")
        @test res == 2
        @test typeof(res) == T

        @test convert(BitStr{2,BigInt}, convert(T, bit"10")) == lbit"10"
        @test convert(BitStr{2,Int64}, convert(T, bit"10")) === bit"10"
    end
end

@testset "indexing" begin
    A = rand(16)
    @test_throws ErrorException A[bit"1101"]
    @test bit"1101"[1] == 1
    @test bit"1101"[2] == 0
    @test bit"1101"[3] == 1
    @test bit"1101"[4] == 1

    mask = Bool[1, 0, 0, 0, 1, 1]
    @test bit"110101"[mask] == [1, 1, 1]
end

@testset "bitstr arithmetic" begin
    x = BitStr64{8}(99)
    y = BitStr64{8}(7)
    @test length(x) == 8
    for op in [+, -, *, ÷, |, ⊻, &, %, mod, mod1]
        z = op(x, y)
        @test z isa BitStr
        @test z == op(Int64(x), Int64(y))
        z = op(x, 7)
        @test z isa BitStr
        @test z == op(Int64(x), 7)
        z = op(7, y)
        @test z isa BitStr
        @test z == op(7, Int64(y))
    end
    @test -x === BitStr64{8}(-99)

    for op in [==, ≈]
        z = op(x, y)
        zz = op(Int64(x), Int64(y))
        @test typeof(z) == typeof(zz)
        @test z == zz
        z = op(x, 7)
        zz = op(Int64(x), 7)
        @test typeof(z) == typeof(zz)
        @test z == zz
        z = op(7, y)
        zz = op(7, Int64(y))
        @test typeof(z) == typeof(zz)
        @test z == zz
    end

    for op in [>, <, >=, <=]
        @test op(x, y) == op(Int64(x), Int64(y))
    end

    for op in [leading_ones, leading_zeros, count_ones, count_zeros, trailing_ones, trailing_zeros]
        op(x) == op(Int64(x))
    end
end

@testset "bitstr constructors" begin
    @test BitStr{8}(99) === BitStr{8,Int}(99)
    x = BitStr64{8}(99)
    y = BitStr64{8}(UInt64(7))
    @test zero(BitStr64{8}) == zero(BitStr64{8}(5)) == BitStr64{8}(0)
    @test reinterpret(BitStr{8,UInt64}, 99) === BitStr{8,UInt64}(99)
    @test reinterpret(Int64, BitStr64{8}(99)) === Int64(99)
    @test all(reinterpret(BitStr64{8}, [Int64(99)]) .=== [BitStr64{8}(99)])
    @test reinterpret(BitStr{8,Int}, bit"1100") === bit"00001100"
    @test x == 99
    @test y == 7
    @test y != BitStr64{4}(7)
    @test bit"1001" == BitStr64{4}(9)
    @test lbit"1001" == bit"1001"
    @test typemax(x) == bit"11111111"
    @test typemax(typeof(x)) == bit"11111111"
    @test typemin(x) == bit"00000000"
    @test typemin(typeof(x)) == bit"00000000"
    @test typemin(BitStr64{4}) == bit"0000"
    @test typemax(BitStr64{4}) == bit"1111"
    @test typemin(x) == bit"00000000"
    @test typemin(typeof(x)) == bit"00000000"
    @test_throws ErrorException x + BitStr64{10}(x)
    @test BitStr64{10}(x) != x
    @test !(BitStr64{10}(x) ≈ x)
end

@testset "bitstr binaryop" begin
    x = bit"001110"
    @test baddrs(x) == [2, 3, 4]
    @test breflect(x) === bit"011100"
    @test breflect(x, [bit"000011", bit"101000"]) === bit"100101"
    msk = bmask(BitStr64{6}, 1, 2, 3)
    @test flip(x, msk) === bit"001001"
    @test setbit(x, msk) === bit"001111"
    @test readbit(x, 3) === bit"000001"
    @test readbit(x, 3, 4) === bit"000011"
    @test anyone(x, msk) === true
    @test allone(x, msk) === false
    @test ismatch(x, msk, bit"000110") === true
    @test ismatch(x, msk, bit"100110") === false
    @test length(basis(x)) == 2^6
    @test length(basis(typeof(x))) == 2^6
    @test basis(x)[1] === bit"000000"

    @test neg(x) === bit"110001"
    @test swapbits(x, bit"101000") === bit"100110"

    @test basis(bit"000") == bit"000":bit"111"
    @test basis(BitStr{3,Int64}) == bit"000":bit"111"
    @test eltype(x) == Int64 == eltype(x[1])
end

@testset "bitstr indexing" begin
    x = bit"001110"
    @test x[1] === Int64(0)
    @test x[2] === Int64(1)
    @test_throws BoundsError x[10]
    @test x[[true, false, true, true, true, false]] == [0, 1, 1, 0]
    @test x[1:2] == [0, 1]
    @test x[[1, 4]] == [0, 1]
    @test lastindex(x) == 6
    @test Base.IteratorSize(x) == Base.HasLength()
end

@testset "bitstr readout" begin
    x = bit"001110"
    @test bint(x) === Int64(14)
    @test bint_r(x) === Int64(28)
    @test bfloat(x) === 28 / 2^6
    @test bfloat_r(x) === 14 / 2^6
    @test btruncate(x, 3) === bit"000110"
end


text:
module GeophysicalFlows

using
  FourierFlows,
  Statistics,
  SpecialFunctions

using FFTW: irfft

include("utils.jl")
include("twodturb.jl")
include("barotropicqg.jl")
include("barotropicqgql.jl")
include("multilayerqg.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Git"
version = v"2.31.0"

# Collection of sources required to build Git
sources = [
    ArchiveSource("https://mirrors.edge.kernel.org/pub/software/scm/git/git-$(version).tar.xz",
                  "e8f162cbdb3283e13cd7388d864ed23485f1b046a19e969f12ed2685fb789a40"),
    ArchiveSource("https://github.com/git-for-windows/git/releases/download/v$(version).windows.1/Git-$(version)-32-bit.tar.bz2",
                  "faff066b2bb4d7640c6a873c87208514218741c201b67b96eacf11be59d9ae1c"; unpack_target = "i686-w64-mingw32"),
    ArchiveSource("https://github.com/git-for-windows/git/releases/download/v$(version).windows.1/Git-$(version)-64-bit.tar.bz2",
                  "a5723d9307e044ab9f20f3d43fd014bc6506378da869c98370d3a98b3de7017f"; unpack_target = "x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""
install_license ${WORKSPACE}/srcdir/git-*/COPYING

if [[ "${target}" == *-ming* ]]; then
    # Fast path for Windows: just copy the content of the tarball to the prefix
    cp -r ${WORKSPACE}/srcdir/${target}/mingw${nbits}/* ${prefix}
    exit
fi

cd $WORKSPACE/srcdir/git-*/

# We need a native "msgfmt" to cross-compile
apk add tcl gettext

# Git doesn't want to be cross-compiled, but we swear that we're not going to
# cross-compile
sed -i 's/cross_compiling=yes/cross_compiling=no/' configure

if [[ "${target}" == *-apple-* ]]; then
    LDFLAGS="-lgettextlib"
elif [[ ${target} == *-freebsd* ]]; then
    LDFLAGS="-L${prefix}/lib -lcharset"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-curl \
    --with-expat \
    --with-openssl \
    --with-iconv=${prefix} \
    --with-libpcre2 \
    --with-zlib=${prefix} \
    CPPFLAGS="-I${prefix}/include" \
    LDFLAGS="${LDFLAGS}"
make -j${nproc}
make install INSTALL_SYMLINKS="yes, please"
"""

platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("git", :git),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("LibCURL_jll"),
    Dependency("Expat_jll", v"2.2.7"; compat="~2.2.7"),
    Dependency("OpenSSL_jll"),
    Dependency("Gettext_jll"),
    Dependency("Libiconv_jll"),
    Dependency("PCRE2_jll", v"10.35.0"),
    Dependency("Zlib_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
package javax.swing;

import java.beans.PropertyChangeEvent;
import java.beans.PropertyChangeListener;
import java.io.IOException;
import java.io.InputStream;
import java.io.HighInputStream;
import java.net.URL;

import javax.swing.event.HyperlinkListener;
import javax.swing.text.JTextComponent;

public class JEditorPane extends JTextComponent implements Scrollable {
    private URL page;
    private InputStream content;
    private HyperlinkListener hlListener;
    private PropertyChangeListener pcListener;

    public URL getPage() {
	return page;
    }

    protected InputStream getStream(URL page) throws IOException {
	InputStream retVal = new HighInputStream();
	return retVal;
    }

    public void setPage(URL newPage) throws IOException {
	URL oldPage = page;
	page = newPage;
	content = getStream(newPage);
	pcListener.propertyChange(new PropertyChangeEvent(this,"page",oldPage,newPage));
    }

    public void addHyperlinkListener(HyperlinkListener l) {
	hlListener = l;
    }

    public void addPropertyChangeListener(String prop, PropertyChangeListener l) {
	pcListener = l;
    }

    public void setContentType(String type) {}
    public void setEditable(boolean editable) {}
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
# Copyright 2015 Conrad Scott

function readUCIBagOfWords(filename :: AbstractString, vocab :: Dictionary)
  header = readtable( filename,
                      header = false, separator = ' ',
                      skipblanks = false, ignorepadding = false,
                      names = [:Header], eltypes = [Int64],
                      nrows = 3 )

  (lines, columns) = size(header)

  if lines != 3
    error("Expected three rows in the header of the UCI bag-of-words file; \"$filename\" only has $lines line(s)")
  end

  if columns != 1
    error("The header of the UCI bag-of-words file should only have one column of data; \"$filename\" has $columns columns")
  end

  (documents, words, triples) = header[:Header]

  if documents < 0
    error("The number of documents in the header should be non-negative; \"$filename\" specifies $documents document(s)")
  end

  if words < 0
    error("The number of words in the header should be non-negative; \"$filename\" specifies $words words(s)")
  end

  if words > length(vocab)
    error("The number of words in the header should not exceed the provided vocabulary; \"$filename\" specifies $words words(s) but the dictionary only has $(length(vocab)) words")
  end

  if triples < 0
    error("The number of triples in the header should be non-negative; \"$filename\" specifies $triples triple(s)")
  end

  body = readtable( filename,
                    header = false, separator = ' ',
                    skipblanks = false, ignorepadding = false,
                    names = [:DocId, :WordId, :Count],
                    eltypes = [Int64, Int64, Int64],
                    skipstart = 3 )

  (lines, columns) = size(body)

  if lines != triples
    error("The header claims that the file contains $triples lines in the body; \"$filename\" has $lines line(s)")
  end

  if columns != 3
    error("Expected three columns in the body of the UCI bag-of-words format file; \"$filename\" has $columns column(s)")
  end

  if minimum(body[:DocId]) < 1
    error("The document ids should be positive integers; \"$filename\" contains a document id of $(minimum(body[:DocId]))")
  end

  if maximum(body[:DocId]) > documents
    error("The header claims that there are $documents documents; \"$filename\" contains a document id of $(maximum(body[:DocId]))")
  end

  if minimum(body[:WordId]) < 1
    error("The word ids should be positive integers; \"$filename\" contains a word id of $(minimum(body[:WordId]))")
  end

  if maximum(body[:WordId]) > words
    error("The header claims that there are $words words; \"$filename\" contains a word id of $(maximum(body[:WordId]))")
  end

  if minimum(body[:Count]) < 0
    error("The document word counts should be non-negative integers; \"$filename\" contains a count of $(minimum(body[:Count])))")
  end

  if maximum(body[:WordId]) > length(vocab)
    error("The word ids should be in the provided dictionary; \"$filename\" contains a word id of $(maximum(body[:WordId])) but the dictionary only has $(length(vocab)) words")
  end

  return Corpus(vocab, sparse(body[:WordId], body[:DocId], body[:Count], words, documents))
end