text:
@noinline not_planar_err() = throw(ArgumentError("not a planar diagram expression"))
@noinline not_planar_err(ex) = throw(ArgumentError("not a planar diagram expression: $ex"))

# Preprocessors used by `@planar` and `@plansor`
function _conj_to_adjoint(ex::Expr)
    if ex.head == :call && ex.args[1] == :conj && TO.istensor(ex.args[2])
        obj, leftind, rightind = TO.decomposetensor(ex.args[2])
        return Expr(:typed_vcat, Expr(TO.prime, obj),
                        Expr(:tuple, rightind...), Expr(:tuple, leftind...))
    else
        return Expr(ex.head, [_conj_to_adjoint(a) for a in ex.args]...)
    end
end
_conj_to_adjoint(ex) = ex

# replacement of TensorOperations functionality:
# adds checks for matching number of domain and codomain indices
# special cases adjoints so that t and t' are considered the same object
# ignore braiding tensors
function _extract_tensormap_objects(ex)
    inputtensors = collect(filter(!=(:τ), _remove_adjoint.(TO.getinputtensorobjects(ex))))
    outputtensors = _remove_adjoint.(TO.getoutputtensorobjects(ex))
    newtensors = TO.getnewtensorobjects(ex)
    (any(==(:τ), newtensors) || any(==(:τ), outputtensors)) &&
        throw(ArgumentError("The name τ is reserved for the braiding, and should not be assigned to."))
    @assert !any(_is_adjoint, newtensors)
    existingtensors = unique!(vcat(inputtensors, outputtensors))
    alltensors = unique!(vcat(existingtensors, newtensors))
    tensordict = Dict{Any,Any}(a => gensym() for a in alltensors)
    pre = Expr(:block, [Expr(:(=), tensordict[a], a) for a in existingtensors]...)
    pre2 = Expr(:block)
    ex = TO.replacetensorobjects(ex) do obj, leftind, rightind
        _is_adj = _is_adjoint(obj)
        if _is_adj
            leftind, rightind = rightind, leftind
            obj = _remove_adjoint(obj)
        end
        newobj = get(tensordict, obj, obj)
        if (obj in existingtensors)
            nl = length(leftind)
            nr = length(rightind)
            nlsym = gensym()
            nrsym = gensym()
            objstr = string(obj)
            errorstr1 = "incorrect number of input-output indices: ($nl, $nr) instead of "
            errorstr2 = " for $objstr."
            checksize = quote
                $nlsym, $nrsym = numout($newobj), numin($newobj)
                ($nlsym == $nl && $nrsym == $nr) ||
                    throw(IndexError($errorstr1 * string(($nlsym, $nrsym)) * $errorstr2))
            end
            push!(pre2.args, checksize)
        end
        return _is_adj ? _add_adjoint(newobj) : newobj
    end
    post = Expr(:block, [Expr(:(=), a, tensordict[a]) for a in newtensors]...)
    pre = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre)
    pre2 = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre2)
    post = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), post)
    return Expr(:block, pre, pre2, ex, post)
end
_is_adjoint(ex) = isa(ex, Expr) && ex.head == TO.prime
_remove_adjoint(ex) = _is_adjoint(ex) ? ex.args[1] : ex
_add_adjoint(ex) = Expr(TO.prime, ex)

# used by `@planar`: realize explicit braiding tensors
function _construct_braidingtensors(ex::Expr)
    if TO.isdefinition(ex) || TO.isassignment(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            list = TO.gettensors(_conj_to_adjoint(rhs))
            if TO.isassignment(ex) && TO.istensor(lhs)
                obj, l, r = TO.decomposetensor(lhs)
                lhs_as_rhs = Expr(:typed_vcat, _add_adjoint(obj),
                                    Expr(:tuple, r...), Expr(:tuple, l...))
                push!(list, lhs_as_rhs)
            end
        else
            return ex
        end
    elseif TO.istensorexpr(ex)
        list = TO.gettensors(_conj_to_adjoint(ex))
    else
        return Expr(ex.head, map(_construct_braidingtensors, ex.args)...)
    end

    i = 1
    translatebraidings = Dict{Any,Any}()
    while i <= length(list)
        t = list[i]
        if _remove_adjoint(TO.gettensorobject(t)) == :τ
            translatebraidings[t] = Expr(:call, GlobalRef(TensorKit, :BraidingTensor))
            deleteat!(list, i)
        else
            i += 1
        end
    end

    unresolved = Any[] # list of indices that we couldn't yet figure out
    indexmap = Dict{Any,Any}()
    # indexmap[i] contains the expression to resolve the space for index i
    for (t, construct_expr) in translatebraidings
        obj, leftind, rightind = TO.decomposetensor(t)
        length(leftind) == length(rightind) == 2 ||
            throw(ArgumentError("The name τ is reserved for the braiding, and should have two input and two output indices."))
        if _is_adjoint(obj)
            i1b, i2b, = leftind
            i2a, i1a, = rightind
        else
            i2b, i1b, = leftind
            i1a, i2a, = rightind
        end

        obj_and_pos1a = _findindex(i1a, list)
        obj_and_pos2a = _findindex(i2a, list)
        obj_and_pos1b = _findindex(i1b, list)
        obj_and_pos2b = _findindex(i2b, list)

        if !isnothing(obj_and_pos1a)
            indexmap[i1b] = Expr(:call, :space, obj_and_pos1a...)
            indexmap[i1a] = Expr(:call, :space, obj_and_pos1a...)
        elseif !isnothing(obj_and_pos1b)
            indexmap[i1b] = Expr(TO.prime, Expr(:call, :space, obj_and_pos1b...))
            indexmap[i1a] = Expr(TO.prime, Expr(:call, :space, obj_and_pos1b...))
        else
            push!(unresolved,(i1a,i1b))
        end

        if !isnothing(obj_and_pos2a)
            indexmap[i2b] = Expr(:call, :space, obj_and_pos2a...)
            indexmap[i2a] = Expr(:call, :space, obj_and_pos2a...)
        elseif !isnothing(obj_and_pos2b)
            indexmap[i2b] = Expr(TO.prime, Expr(:call, :space, obj_and_pos2b...))
            indexmap[i2a] = Expr(TO.prime, Expr(:call, :space, obj_and_pos2b...))
        else
            push!(unresolved,(i2a,i2b))
        end
    end
    # simple loop that tries to resolve as many indices as possible
    changed = true
    while changed == true
        changed = false
        i = 1
        while i <= length(unresolved)
            (i1, i2) = unresolved[i]
            if i1 in keys(indexmap)
                changed = true
                indexmap[i2] = indexmap[i1]
                deleteat!(unresolved, i)
            elseif i2 in keys(indexmap)
                changed = true
                indexmap[i1] = indexmap[i2]
                deleteat!(unresolved, i)
            else
                i += 1
            end
        end
    end
    !isempty(unresolved) &&
        throw(ArgumentError("cannot determine the spaces of indices " *
                                string(tuple(unresolved...)) *
                                    "for the braiding tensors in $(ex)"))

    pre = Expr(:block)
    for (t, construct_expr) in translatebraidings
        obj, leftind, rightind = TO.decomposetensor(t)
        if _is_adjoint(obj)
            i1b, i2b, = leftind
            i2a, i1a, = rightind
        else
            i2b, i1b, = leftind
            i1a, i2a, = rightind
        end
        push!(construct_expr.args, indexmap[i1b])
        push!(construct_expr.args, indexmap[i2b])
        s = gensym()
        push!(pre.args, Expr(:(=), s, construct_expr))
        ex = TO.replacetensorobjects(ex) do o, l, r
            if o == obj && l == leftind && r == rightind
                return obj  == :τ ? s : Expr(TO.prime, s)
            else
                return o
            end
        end
    end
    return Expr(:block, pre, ex)
end
_construct_braidingtensors(x) = x

# used by non-planar parser of `@plansor`: remove explicit braiding tensors
function _remove_braidingtensors(ex::Expr)
    outgoing = []

    if TO.isdefinition(ex) || TO.isassignment(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            list = TO.gettensors(_conj_to_adjoint(rhs))
            if TO.istensor(lhs)
                obj, l, r = TO.decomposetensor(lhs)
                outgoing = [l; r]
            end
        else
            return ex
        end
    elseif TO.istensorexpr(ex)
        list = TO.gettensors(_conj_to_adjoint(ex))
    else
        return Expr(ex.head, map(_remove_braidingtensors, ex.args)...)
    end

    τs = Any[]
    i = 1
    while i <= length(list)
        t = list[i]
        if _remove_adjoint(TO.gettensorobject(t)) == :τ
            push!(τs, t)
            deleteat!(list, i)
        else
            i += 1
        end
    end

    indexmap = Dict{Any,Any}()
    # to remove the braidingtensors, we need to map certain indices to other indices
    for t in τs
        obj, leftind, rightind = TO.decomposetensor(t)
        length(leftind) == length(rightind) == 2 ||
            throw(ArgumentError("The name τ is reserved for the braiding, and should have two input and two output indices."))
        if _is_adjoint(obj)
            i1b, i2b, = leftind
            i2a, i1a, = rightind
        else
            i2b, i1b, = leftind
            i1a, i2a, = rightind
        end

        i1a = get(indexmap, i1a, i1a)
        i1b = get(indexmap, i1b, i1b)
        i2a = get(indexmap, i2a, i2a)
        i2b = get(indexmap, i2b, i2b)

        obj_and_pos1a = _findindex(i1a, list)
        obj_and_pos2a = _findindex(i2a, list)
        obj_and_pos1b = _findindex(i1b, list)
        obj_and_pos2b = _findindex(i2b, list)

        if i1a in outgoing
            indexmap[i1a] = i1a
            indexmap[i1b] = i1a
        elseif i1b in outgoing
            indexmap[i1a] = i1b
            indexmap[i1b] = i1b
        else
            if i1a isa Int && i1b isa Int
                indexmap[i1a] = max(i1a, i1b)
                indexmap[i1b] = max(i1a, i1b)
            else
                indexmap[i1a] = i1a
                indexmap[i1b] = i1a
            end
        end

        if i2a in outgoing
            indexmap[i2a] = i2a
            indexmap[i2b] = i2a
        elseif i2b in outgoing
            indexmap[i2a] = i2b
            indexmap[i2b] = i2b
        else
            if i2a isa Int && i2b isa Int
                indexmap[i2a] = max(i2a, i2b)
                indexmap[i2b] = max(i2a, i2b)
            else
                indexmap[i2a] = i2a
                indexmap[i2b] = i2a
            end
        end
    end

    # simple loop that tries to simplify the indicemaps (a=>b,b=>c -> a=>c,b=>c)
    changed = true
    while changed == true
        changed = false
        i = 1
        for (k,v) in indexmap
            if v in keys(indexmap) && indexmap[v] != v
                changed = true
                indexmap[k] = indexmap[v]
            end
        end
    end

    ex = TO.replaceindices(i -> get(indexmap, i, i), ex)
    return _purge_braidingtensors(ex)
end
_remove_braidingtensors(x) = x

function _purge_braidingtensors(ex::Expr) # actually remove the braidingtensors
    args = collect(filter(ex.args) do a
        if a isa Expr && a.head == :call && a.args[1] == :conj
            a = a.args[2]
        end
        if a isa Expr && TO.istensor(a) && _remove_adjoint(TO.gettensorobject(a)) == :τ
            _, leftind, rightind = TO.decomposetensor(a)
            (leftind[1] == rightind[2] && leftind[2] == rightind[1]) ||
                throw(ArgumentError("unable to remove braiding tensor $a"))
            return false
        end
        return true
    end)

    # multiplication with only a single argument is (rightfully) seen as invalid syntax
    if ex.head == :call && args[1] == :* && length(args) == 2
        return _purge_braidingtensors(args[2])
    else
        return Expr(ex.head, map(_purge_braidingtensors, args)...)
    end
end
_purge_braidingtensors(x) = x

function _check_planarity(ex::Expr)
    if ex.head == :macrocall && ex.args[1] == Symbol("@notensor")
    elseif TO.isassignment(ex) || TO.isdefinition(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            if TO.istensorexpr(lhs)
                @assert TO.istensor(lhs)
                indlhs = only(get_possible_planar_indices(lhs)) # should have only one element
            else
                indlhs = Any[]
            end
            indsrhs = get_possible_planar_indices(rhs)
            isempty(indsrhs) && not_planar_err(rhs)
            i = findfirst(ind -> iscyclicpermutation(ind, indlhs), indsrhs)
            i === nothing && not_planar_err(ex)
        end
    else
        foreach(ex.args) do a
            _check_planarity(a)
        end
    end
    return ex
end
_check_planarity(ex) = ex

# decompose contraction trees in order to fix index order of temporaries
# to ensure that planarity is guaranteed
_decompose_planar_contractions(ex, temporaries) = ex
function _decompose_planar_contractions(ex::Expr, temporaries)
    if ex.head == :macrocall && ex.args[1] == Symbol("@notensor")
        return ex
    end
    if TO.isassignment(ex) || TO.isdefinition(ex)
        lhs, rhs = TO.getlhs(ex), TO.getrhs(ex)
        if TO.istensorexpr(rhs)
            pre = Vector{Any}()
            rhs = _extract_contraction_pairs(rhs, lhs, pre, temporaries)
            return Expr(:block, pre..., Expr(ex.head, lhs, rhs))
        else
            return ex
        end
    end
    if TO.istensorexpr(ex)
        pre = Vector{Any}()
        rhs = _extract_contraction_pairs(ex, (Any[], Any[]), pre, temporaries)
        return Expr(:block, pre..., rhs)
    end
    if ex.head == :block
        return Expr(ex.head,
                    [_decompose_planar_contractions(a, temporaries) for a in ex.args]...)
    end
    if ex.head == :for || ex.head == :function
        return Expr(ex.head, ex.args[1],
                        _decompose_planar_contractions(ex.args[2], temporaries))
    end
    return ex
end

# decompose a contraction into elementary binary contractions of tensors without inner traces
# if lhs is an expression, it contains the existing lhs and thus the index order
# if lhs is a tuple, the result is a temporary object and the tuple (lind, rind) gives a suggestion for the preferred index order
function _extract_contraction_pairs(rhs, lhs, pre, temporaries)
    if TO.isscalarexpr(rhs)
        return rhs
    elseif TO.isgeneraltensor(rhs)
        if TO.hastraceindices(rhs) && lhs isa Tuple
            s = gensym()
            newlhs = Expr(:typed_vcat, s, Expr(:tuple, lhs[1]...), Expr(:tuple, lhs[2]...))
            push!(temporaries, s)
            push!(pre, Expr(:(:=), newlhs, rhs))
            return newlhs
        else
            return rhs
        end
    elseif rhs.head == :call && rhs.args[1] == :*
        @assert length(rhs.args) == 3

        if lhs isa Expr
            _, leftind, rightind = TO.decomposetensor(lhs)
        else
            leftind, rightind = lhs
        end
        lhs_ind = vcat(leftind, reverse(rightind))

        # find possible planar order
        rhs_inds = Any[]
        for ind1 in get_possible_planar_indices(rhs.args[2])
            for ind2 in get_possible_planar_indices(rhs.args[3])
                for (oind1, oind2, cind1, cind2) in possible_planar_complements(ind1, ind2)
                    if iscyclicpermutation(vcat(oind1, oind2), lhs_ind)
                        push!(rhs_inds, (ind1, ind2, oind1, oind2, cind1, cind2))
                    end
                    isempty(rhs_inds) || break
                end
                isempty(rhs_inds) || break
            end
            isempty(rhs_inds) || break
        end
        ind1, ind2, oind1, oind2, cind1, cind2 = only(rhs_inds) # inds_rhs should hold exactly one match
        if all(in(leftind), oind2) && all(in(rightind), oind1) # reverse order
            a1 = _extract_contraction_pairs(rhs.args[3], (oind2, reverse(cind2)), pre, temporaries)
            a2 = _extract_contraction_pairs(rhs.args[2], (cind1, reverse(oind1)), pre, temporaries)
            oind1, oind2 = oind2, oind1
            cind1, cind2 = cind2, cind1
        else
            a1 = _extract_contraction_pairs(rhs.args[2], (oind1, reverse(cind1)), pre, temporaries)
            a2 = _extract_contraction_pairs(rhs.args[3], (cind2, reverse(oind2)), pre, temporaries)
        end

        if TO.isscalarexpr(a1) || TO.isscalarexpr(a2)
            rhs = Expr(:call, :*, a1, a2)
            s = gensym()
            newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...),
                                        Expr(:tuple, reverse(oind2)...))
            push!(temporaries, s)
            push!(pre, Expr(:(:=), newlhs, rhs))
            return newlhs
        end

        # note that index order in _extract... is only a suggestion, now we have actual index order
        _, l1, r1, = TO.decomposegeneraltensor(a1)
        _, l2, r2, = TO.decomposegeneraltensor(a2)
        if all(in(r1), oind1) && all(in(l2), oind2) # reverse order
            a1, a2 = a2, a1
            ind1, ind2 = ind2, ind1
            oind1, oind2 = oind2, oind1
        end
        if lhs isa Tuple
            rhs = Expr(:call, :*, a1, a2)
            s = gensym()
            newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...),
                                        Expr(:tuple, reverse(oind2)...))
            push!(temporaries, s)
            push!(pre, Expr(:(:=), newlhs, rhs))
            return newlhs
        else
            if leftind == oind1 && rightind == reverse(oind2)
                rhs = Expr(:call, :*, a1, a2)
                return rhs
            elseif leftind == oind2 && rightind == reverse(oind1) # probably this can not happen anymore
                rhs = Expr(:call, :*, a2, a1)
                return rhs
            else
                rhs = Expr(:call, :*, a1, a2)
                s = gensym()
                newlhs = Expr(:typed_vcat, s, Expr(:tuple, oind1...),
                                            Expr(:tuple, reverse(oind2)...))
                push!(temporaries, s)
                push!(pre, Expr(:(:=), newlhs, rhs))
                return newlhs
            end
        end
    elseif rhs.head == :call && rhs.args[1] ∈ (:+, :-)
        args = [_extract_contraction_pairs(a, lhs, pre, temporaries) for
                    a in rhs.args[2:end]]
        return Expr(rhs.head, rhs.args[1], args...)
    else
        throw(ArgumentError("unknown tensor expression"))
    end
end

function _findindex(i, list) # finds an index i in a list of tensor expressions
    for t in list
        obj, l, r = TO.decomposetensor(t)
        pos = findfirst(==(i), l)
        isnothing(pos) || return (obj, pos)
        pos = findfirst(==(i), r)
        isnothing(pos) || return (obj, pos + length(l))
    end
    return nothing
end


text:
# Build the real docs first.
include(joinpath(dirname(@__FILE__), "..", "docs", "make.jl"))

# docs module
# ===========

module Mod

"""
    func(x)

[`T`](@ref)
"""
func(x) = x

"""
    T

[`func(x)`](@ref)
"""
type T end

end

# autodocs module
# ===============

"`AutoDocs` module."
module AutoDocs

module Pages

include("pages/a.jl")
include("pages/b.jl")
include("pages/c.jl")
include("pages/d.jl")
include("pages/e.jl")

end

"Function `f`."
f(x) = x

"Constant `K`."
const K = 1

"Type `T`."
type T end

"Macro `@m`."
macro m() end

"Module `A`."
module A

"Function `A.f`."
f(x) = x

"Constant `A.K`."
const K = 1

"Type `B.T`."
type T end

"Macro `B.@m`."
macro m() end

end

"Module `B`."
module B

"Function `B.f`."
f(x) = x

"Constant `B.K`."
const K = 1

"Type `B.T`."
type T end

"Macro `B.@m`."
macro m() end

end

end

# tests module
# ============

module Tests

using Documenter
using Base.Test
using Compat

# Unit tests for module internals.

module UnitTests

module SubModule end

# Does `submodules` collect *all* the submodules?
module A
module B
module C
module D
end
end
end
end

type T end

"Documenter unit tests."
Base.length(::T) = 1

end

let doc = @doc(length)
    a = Documenter.Utilities.filterdocs(doc, Set{Module}())
    b = Documenter.Utilities.filterdocs(doc, Set{Module}([UnitTests]))
    c = Documenter.Utilities.filterdocs(doc, Set{Module}([Base]))
    d = Documenter.Utilities.filterdocs(doc, Set{Module}([Tests]))

    @test !isnull(a)
    @test get(a) === doc
    @test !isnull(b)
    @test contains(stringmime("text/plain", get(b)), "Documenter unit tests.")
    @test !isnull(c)
    @test !contains(stringmime("text/plain", get(c)), "Documenter unit tests.")
    @test isnull(d)
end

# Documenter.Utilities.issubmodule
@test Documenter.Utilities.issubmodule(Main, Main) === true
@test Documenter.Utilities.issubmodule(UnitTests, UnitTests) === true
@test Documenter.Utilities.issubmodule(UnitTests.SubModule, Main) === true
@test Documenter.Utilities.issubmodule(UnitTests.SubModule, UnitTests) === true
@test Documenter.Utilities.issubmodule(UnitTests.SubModule, Base) === false
@test Documenter.Utilities.issubmodule(UnitTests, UnitTests.SubModule) === false

@test UnitTests.A in Documenter.Utilities.submodules(UnitTests.A)
@test UnitTests.A.B in Documenter.Utilities.submodules(UnitTests.A)
@test UnitTests.A.B.C in Documenter.Utilities.submodules(UnitTests.A)
@test UnitTests.A.B.C.D in Documenter.Utilities.submodules(UnitTests.A)

## NavNode tests

include("navnode.jl")

# DocSystem unit tests.

import Documenter: DocSystem

## Bindings.

@test_throws ArgumentError DocSystem.binding(9000)
let b = Docs.Binding(current_module(), :DocSystem)
    @test DocSystem.binding(b) == b
end
let b = DocSystem.binding(Documenter.Documents.Document)
    @test b.mod === Documenter.Documents
    @test b.var === :Document
end
let b = DocSystem.binding(Documenter)
    @test b.mod === Main
    @test b.var === :Documenter
end
let b = DocSystem.binding(:Main)
    # @test b.mod === Main
    @test b.var === :Main
end
let b = DocSystem.binding(DocSystem.binding)
    @test b.mod === DocSystem
    @test b.var === :binding
end
let b = DocSystem.binding(getfield(Core.Intrinsics, :ccall))
    @test b.mod === Core.Intrinsics
    @test b.var === :ccall
end
let b = DocSystem.binding(Documenter, :Documenter)
    @test b.mod === Main
    @test b.var === :Documenter
end

## `MultiDoc` object.

@test isdefined(DocSystem, :MultiDoc)
@test fieldnames(DocSystem.MultiDoc) == [:order, :docs]

## `DocStr` object.

@test isdefined(DocSystem, :DocStr)
@test fieldnames(DocSystem.DocStr) == [:text, :object, :data]

## `getdocs`.

const alias_of_getdocs = DocSystem.getdocs

let b   = DocSystem.binding(DocSystem, :getdocs),
    d_0 = DocSystem.getdocs(b, Tuple{}),
    d_1 = DocSystem.getdocs(b),
    d_2 = DocSystem.getdocs(b, Union{Tuple{ANY}, Tuple{ANY, Type}}; compare = (==)),
    d_3 = DocSystem.getdocs(b; modules = Module[Main]),
    d_4 = DocSystem.getdocs(DocSystem.binding(current_module(), :alias_of_getdocs)),
    d_5 = DocSystem.getdocs(DocSystem.binding(current_module(), :alias_of_getdocs); aliases = false)

    @test length(d_0) == 0
    @test length(d_1) == 2
    @test length(d_2) == 1
    @test length(d_3) == 0
    @test length(d_4) == 2
    @test length(d_5) == 0

    @test d_1[1].data[:binding] == b
    @test d_1[2].data[:binding] == b
    @test d_1[1].data[:typesig] == Union{Tuple{Docs.Binding}, Tuple{Docs.Binding, Type}}
    @test d_1[2].data[:typesig] == Union{Tuple{ANY}, Tuple{ANY, Type}}
    @test d_1[1].data[:module]  == DocSystem
    @test d_1[2].data[:module]  == DocSystem

    @test d_2[1].data[:binding] == b
    @test d_2[1].data[:typesig] == Union{Tuple{ANY}, Tuple{ANY, Type}}
    @test d_2[1].data[:module]  == DocSystem

    @test d_1 == d_4
    @test d_1 != d_5
end

## DOM Tests.

module DOMTests

using Base.Test
import Documenter.Utilities.DOM: DOM, @tags

@tags div ul li p

for tag in (:div, :ul, :li, :p)
    TAG = getfield(current_module(), tag)
    @test isdefined(tag)
    @test isa(TAG, DOM.Tag)
    @test TAG.name === tag
end

@test div().name === :div
@test div().text == ""
@test isempty(div().nodes)
@test isempty(div().attributes)

@test div("...").name === :div
@test div("...").text == ""
@test length(div("...").nodes) === 1
@test div("...").nodes[1].text == "..."
@test div("...").nodes[1].name === Symbol("")
@test isempty(div("...").attributes)

@test div[".class"]("...").name === :div
@test div[".class"]("...").text == ""
@test length(div[".class"]("...").nodes) === 1
@test div[".class"]("...").nodes[1].text == "..."
@test div[".class"]("...").nodes[1].name === Symbol("")
@test length(div[".class"]("...").attributes) === 1
@test div[".class"]("...").attributes[1] == (:class => "class")
@test div[:attribute].attributes[1] == (:attribute => "")
@test div[:attribute => "value"].attributes[1] == (:attribute => "value")

let d = div(ul(map(li, [string(n) for n = 1:10])))
    @test d.name === :div
    @test d.text == ""
    @test isempty(d.attributes)
    @test length(d.nodes) === 1
    let u = d.nodes[1]
        @test u.name === :ul
        @test u.text == ""
        @test isempty(u.attributes)
        @test length(u.nodes) === 10
        for n = 1:10
            let v = u.nodes[n]
                @test v.name === :li
                @test v.text == ""
                @test isempty(v.attributes)
                @test length(v.nodes) === 1
                @test v.nodes[1].name === Symbol("")
                @test v.nodes[1].text == string(n)
                @test !isdefined(v.nodes[1], :attributes)
                @test !isdefined(v.nodes[1], :nodes)
            end
        end
    end
end

@tags script style img

@test string(div(p("one"), p("two"))) == "

one

two

"
@test string(div[:key => "value"])    == "

"
@test string(p(" < > & ' \" "))       == "
 < > & ' "
"
@test string(img[:src => "source"])   == ""
@test string(img[:none])              == ""
@test string(script(" < > & ' \" "))  == ""
@test string(style(" < > & ' \" "))   == ""
@test string(script)                  == "

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Basics (Modia3D/Basics/_module.jl)
#

"""
    path = PTP_path(names;
                    positions = [zeros(size(names,1))'; ones(size(names,1))'],
                    startTime=0.0,
                    v_max = ones(size(names,1)),
                    a_max = ones(size(names,2))))

Generate a new path object.
"""
mutable struct PTP_path
   names::Vector{AbstractString}
   startTime::Float64
   v_max::Vector{Float64}
   a_max::Vector{Float64}
   positions::Matrix{Float64}

   delta::Matrix{Float64}
   hasPath::Vector{Bool}
   sd_max::Vector{Float64}
   sdd_max::Vector{Float64}
   Ta1::Vector{Float64}
   Ta2::Vector{Float64}
   noWphase::Vector{Bool}
   Tv::Vector{Float64}
   Te::Vector{Float64}
   Ta1s::Vector{Float64}
   Ta2s::Vector{Float64}
   Tvs::Vector{Float64}
   Tes::Vector{Float64}
   sd_max2::Vector{Float64}
   s1::Vector{Float64}
   s2::Vector{Float64}
   s3::Vector{Float64}
   Tend::Float64
   posTemp::Vector{Float64}   # Temporary storage that can be used by the functions operation on PTP_path

   function PTP_path(names::AbstractVector;
                     positions::Matrix{Float64} = [zeros(size(names,1))'; ones(size(names,1))'],
                     startTime=0.0,
                     v_max = ones(size(names,1)),
                     a_max = ones(size(names,2)))

        #@assert(size(positions,1) > 1)
        #@assert(size(positions,2) == Base.length(names))
        #@assert(Base.length(v_max) == Base.length(names))
        #@assert(Base.length(a_max) == Base.length(names))
        np    = Base.length(names)
        npath = size(positions,1) - 1   # number of path points

        for i in eachindex(v_max)
            @assert(v_max[i] > 0.0)
            @assert(a_max[i] > 0.0)
        end

        delta    = zeros(npath, np)
        hasPath  = fill(false,npath)
        sd_max   = zeros(npath)
        sdd_max  = zeros(npath)
        Ta1      = zeros(npath)
        Ta2      = zeros(npath)
        noWphase = fill(false,npath)
        Tv       = zeros(npath)
        Te       = zeros(npath)
        Ta1s     = zeros(npath)
        Ta2s     = zeros(npath)
        Tvs      = zeros(npath)
        Tes      = zeros(npath)
        sd_max2  = zeros(npath)
        s1       = zeros(npath)
        s2       = zeros(npath)
        s3       = zeros(npath)
        aux1     = zeros(np)
        aux2     = zeros(np)

        for i in 1:npath
            delta[i,:] = positions[i+1,:] - positions[i,:]

            for j in 1:np
                aux1[j] = delta[i,j]/v_max[j]
                aux2[j] = delta[i,j]/a_max[j]
            end

            sd_max_inv  = maximum(abs.(aux1))
            sdd_max_inv = maximum(abs.(aux2))

            small = 1000*eps()
            if sd_max_inv <= small || sdd_max_inv <= small
                hasPath[i] = false
            else
                hasPath[i] = true
                sd_max[i]  = 1/sd_max_inv
                sdd_max[i] = 1/sdd_max_inv

                Ta1[i] = sqrt(1/sdd_max[i])
                Ta2[i] = sd_max[i]/sdd_max[i]
                noWphase[i] = Ta2[i] >= Ta1[i]

                Tv[i]   = noWphase[i] ? Ta1[i] : 1/sd_max[i]
                Te[i]   = noWphase[i] ? Ta1[i] + Ta1[i] : Tv[i] + Ta2[i]

                Tbegin  = i==1 ? startTime : Tes[i-1]
                Ta1s[i] = Ta1[i] + Tbegin
                Ta2s[i] = Ta2[i] + Tbegin
                Tvs[i]  = Tv[i]  + Tbegin
                Tes[i]  = Te[i]  + Tbegin
                sd_max2[i] = sdd_max[i]*Ta1[i]
                s1[i] = sdd_max[i]*(noWphase[i] ? Ta1[i]*Ta1[i] : Ta2[i]*Ta2[i])/2
                s2[i] = s1[i] + (noWphase[i] ? sd_max2[i]*(Te[i] - Ta1[i]) - (sdd_max[i]/2)*(Te[i] - Ta1[i])^2 : sd_max[i]*(Tv[i] - Ta2[i]))
                s3[i] = s2[i] + sd_max[i]*(Te[i] - Tv[i]) - (sdd_max[i]/2)*(Te[i] - Tv[i])^2
            end
        end

        Tend = Tes[end]

        new(names, startTime, v_max, a_max, positions, delta,
            hasPath, sd_max, sdd_max, Ta1, Ta2, noWphase, Tv, Te, Ta1s, Ta2s, Tvs, Tes, sd_max2, s1, s2, s3, Tend, zeros(np))
    end
end

"""
    getPosition!(path, time, position)

Given a `path::PTP_path` and a time instant `time`, return the actual
position at time `time` in vector `position`.
"""
function getPosition!(path::PTP_path, time::Number, position::Vector{Float64})
    time = ustrip(time)
    @assert(length(position) == size(path.positions,2))
    npath = length(path.hasPath)
    np    = length(position)

    # Search correct time interval
    i = 0
    if time <= path.startTime
        i = 1
        s = 0

    else
        while i < npath
            i = i+1
            if time <= path.Tes[i]
               break
            end
        end
        if time >= path.Tes[end]
            i = npath
            s = path.noWphase[i] ? path.s2[end] : path.s3[end]
            #println("... time=$time i=$i s=$s qbegin=", path.positions[i,1], ", qdelta = ", path.delta[i,1])
        else
            Tbegin = i==1 ? path.startTime : path.Tes[i-1]
            if path.noWphase[i]
                if time < path.Ta1s[i]
                    s = (path.sdd_max[i]/2)*(time - Tbegin)^2
                elseif time < path.Tes[i]
                    s = path.s1[i] + path.sd_max2[i]*(time - path.Ta1s[i]) - (path.sdd_max[i]/2)*(time - path.Ta1s[i])^2
                else
                    s = path.s2[i]
                end
            elseif time < path.Ta2s[i]
                s = (path.sdd_max[i]/2)*(time - Tbegin)^2
            elseif time < path.Tvs[i]
                s = path.s1[i] + path.sd_max[i]*(time - path.Ta2s[i])
            elseif time < path.Tes[i]
                s = path.s2[i] + path.sd_max[i]*(time - path.Tvs[i]) - (path.sdd_max[i]/2)*(time - path.Tvs[i])*(time - path.Tvs[i])
            else
                s = path.s3[i]
            end
        end
    end

    for j in 1:np
        # println("... i = $i, j= $j, np = $np, s=$s, positions[i,1]=",path.positions[i,1])
        position[j]= path.positions[i,j] + path.delta[i,j]*s
    end
end

"""
    pos = getPosition(path, index, time)

Given a `path::PTP_path`, the `index` of a signal, and a time instant `time`, return the actual
position at time `time`.
"""
function getPosition(path::PTP_path, index, time::Number)
    getPosition!(path, time, path.posTemp)
    return path.posTemp[index]
end

"""
    index = getIndex(path, name)

Return the index of `name` in `path` or trigger an error, if not present.
"""
function getIndex(path, name)
    index = findfirst(x -> x==name, path.names)
    if typeof(index) == Nothing
        error("getIndex(path,name): \"", name, "\" not in path")
    end
    return index
end

"""
    plotPath(path;
             names=path.names, heading="PTP plots",
             tend=1.1*path.Tend, figure=1, ntime=101)

Given a `path::PTP_path`, plot the path over `time` up to `tend` for all points
identified by the vector or tuple `names` to figure `figure`
using `ntime` time points.
"""
function plotPath(path::PTP_path; names=path.names, heading="PTP plots", figure=1, ntime=101, tend = 1.1*path.Tend)
    time = range(0u"s",(tend)u"s",length=ntime)
    indices = indexin(names, path.names)
    names2  = deepcopy(names)
    for i in eachindex(indices)
        if typeof(i) == Nothing
            @warn "plotPath(path, ...): \""*names[i]*"\" is ignored, because not in path"
            deleteat!(indices,i)
            deleteat!(names2,i)
        end
    end

    np   = length(indices)
    q    = zeros(length(time), np)
    qt   = zeros(length(path.names))
    for i in eachindex(time)
        getPosition!(path, time[i], qt)
        q[i,:] = qt[indices]
    end

    series = Dict{AbstractString,Any}()
    series["time"] = time
    for i in eachindex(names2)
        series[names2[i]] = q[:,i]
    end

    ModiaMath.plot(series, Tuple(names2), heading=heading, figure=figure)
end


text:
# targets for phase in / phase out of policies

pinft = 0.375.*pinf_o # infection prob at meeting

socialmaxyyt = 2 # max no. people met outside firm young_young
socialmaxoyt = 1 # max no. people met outside firm old_young
socialmaxoot = 0.5 # max no. people met outside firm old_old
socialmaxyot = 0.5 # max no. people met outside firm old_old

phomeofficet = 1

pshopt = [[1,1,1],[1,1,1]] - alpha_lockdown*[[0.1,0.33,0],[0.1,0.33,0]]
pshop = [[1,1,1],[1,1,1]] - alpha_lockdown*[[0.1,0.33,0],[0.1,0.33,0]]

workmeett = [4,5,8,2] # max number of people met at work: man / ser / food /public
shopmeett = [5,20,10] # max number of people met at each shopping instant  man / ser / food

# direct change in variable through policy
tauadj = 0

genhomeoffice = true

ptravel = 0

bailoutprogram = true # initially no bailouts

shorttimeprogram = true # initially no short term program

# If this is the first (non-adaptive) lockdown, set policy indicator to 'bad' and counter to 1
if polswitchcount == 0
    currentadaptivepolicy = "BAD"
    polswitchvount = 1
end


text:
n = 5 # rand(1:5)
A = matrixdepot("wathen", n)

@test typeof(A) <: SparseMatrixCSC
@test issymmetric(Matrix(A))
@test isposdef(Matrix(A))

println("'wathen' passed test...")


text:
"""
素数判定
"""
function is_prime_simple(n::Int64)
    if n < 2
        return false
    elseif n == 2
        return true
    else
        for i in 2:n-1
            #println(i, "...")
            if n % i == 0
                return false
            end
            #println("わりきれない")
        end
    end
    return true
end

function do_2(n₀::Int64, n₁::Int64)
    ns = n₀:n₁
    prime_list = is_prime_simple.(ns)
    num = prime_list' * prime_list
    println(n₀, " から ", n₁ ," に素数は", num, " 個ある")
end

#is_prime_simple(999961)
#@time is_prime_simple(999961)  # 100万以下で最大の素数

#do_2(1, 10)
@time do_2(1, 100000)


text:
Base.@kwdef struct SmoothAnalysis
    npoints::Int=200
    span::Float64=0.75
    degree::Int=2
end

function (l::SmoothAnalysis)(input::ProcessedLayer)
    output = map(input) do p, _
        x, y = p
        model = Loess.loess(x, y; l.span, l.degree)
        x̂ = range(extrema(x)..., length=l.npoints)
        ŷ = Loess.predict(model, x̂)
        return (x̂, ŷ), (;)
    end
    plottype = Makie.plottype(output.plottype, Lines)
    return ProcessedLayer(output; plottype)
end

"""
    smooth(; span=0.75, degree=2, npoints=200)

Fit a loess model. `span` is the degree of smoothing, typically in `[0,1]`.
Smaller values result in smaller local context in fitting.
`degree` is the polynomial degree used in the loess model.
`npoints` is the number of points used by Makie to draw the line
"""
smooth(; options...) = transformation(SmoothAnalysis(; options...))


text:
function apodizer(data::Array{Float64, 2})
    (Nx, Ny) = size(data)
    Amat = wind_2d(Nx)
    datad_w = fweights(Amat)
    meanVal = mean(data, datad_w) #<AF>
    temp2d_a = (data.-meanVal).*Amat.+meanVal #A(F-μ) + μ
    return temp2d_a
end

function wind_2d(nx)
    dx   = nx/2-1
    filter = zeros(Float64, nx, nx)
    A = DSP.tukey(nx, 0.3)
    itp = extrapolate(interpolate(A,BSpline(Linear())),0)
    @inbounds for x = 1:nx
        sx = x-dx-1    # define sx,sy so that no fftshift() needed
        for y = 1:nx
            sy = y-dx-1
            r  = sqrt.((sx).^2 + (sy).^2) + nx/2
            filter[x,y] = itp(r)
        end
    end
    return filter
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
#===============================================================================
# Practice with the julia function syntax and construction.
#
# Author: Carl Cortright
# Date: 12/26/2016
===============================================================================#

# Function that completes one operation, adding x and y
function f(x, y)
  x+y
end

println(f(1,2))

# In assignment form
f(x, y) = x + y

# f can be passed around like a value, it refers to the function object
g = f
println(g(200, 21))

# Oh, and this is awesome
∑(x, y, z) = x + y + z
println(∑(1,2,3))

# Example of where return is useful
function hypotCustom(x, y)
  a = abs(x)
  b = abs(y)
  if x > y
    r = y/x
    return x * sqrt(1 + r^2)
  elseif y === 0
    return 0
  else
    r = x/y
    return y * sqrt(1 + r^2)
  end
end

if hypotCustom(3,2) === sqrt(13)
  println("got it")
else
  println(hypotCustom(3,2))
end

# Using the + operator as a function
println(+(1,2,3))
# Oooh, and you can pass operators
h = +
println(h(1,2,3))

# Playing with multiple return values
function foo(x, y)
  return x+y, x*y
end

# Creating functions with a variable number of arguments
function variableArgs(a,b,x...)
  return (a,b,x)
end

# Optional arguments
function optional(debug=true)
  if debug === true
    println("Debug is on")
  end
end
optional()
optional(false)

#
# In this case, c is optional and can be specified using a keyword
#
function keywordSum(a, b; c=0)
  return a+b+c
end
println(keywordSum(1,2; c=8))

max_issues_repo_name: carlcortright/Playground


text:
using Base.Test
using LinearResponseVariationalBayes.ExponentialFamilies
import Distributions

println("Testing ExponentialFamilies.")

function rel_diff(x, y)
	(x - y) ./ y
end

function test_wishart()
	srand(42)
	p = 2
	p_ud = round(Int, p * (p + 1) / 2)
	mu_sigma = eye(p) + repmat([0.2], p, p)
	lambda = inv(mu_sigma)
	n_draws = 500_000

	wn = 100.
	v0 = inv(lambda)

	w = Distributions.Wishart(wn, v0)
	w_draws = [ rand(w) for i=1:n_draws ];
	w_logdet_draws = Float64[ logdet(w_draw) for w_draw in w_draws ];

	# Test the means.
	@test_approx_eq_eps(rel_diff(mean(w_logdet_draws), wishart_e_log_det(wn, inv(v0))),
	                    0.0, 1e-3)

	@test_approx_eq_eps(rel_diff(mean(w_draws), (v0 * wn)), zeros(p, p), 1e-2)

	@test_approx_eq(ExponentialFamilies.wishart_entropy(wn, lambda, p),
	                Distributions.entropy(w))

	# Test the covariances.
	ud_ind = make_ud_index_matrix(p)
	w_lin = Array(Float64, n_draws, p_ud)
	for row = 1:n_draws
		w_lin[row, :] = linearize_matrix(w_draws[row], ud_ind)
	end

	wishart_cov =
		ExponentialFamilies.get_wishart_variational_covariance(v0, wn, ud_ind)
	@test_approx_eq_eps(rel_diff(cov(w_lin), wishart_cov), zeros(p_ud, p_ud), 5e-2)

	wishart_suff_stats = hcat(w_lin, w_logdet_draws);
	wishart_suff_stats_sample_cov = cov(wishart_suff_stats)
	wishart_cross_sample_cov = wishart_suff_stats_sample_cov[p_ud + 1,:][1:p_ud]
	wishart_cross_cov = ExponentialFamilies.get_wishart_cross_variance(v0, ud_ind)
	@test_approx_eq_eps(rel_diff(wishart_cross_sample_cov, wishart_cross_cov),
	                    zeros(p_ud), 1e-1)

	wishart_log_det_var = ExponentialFamilies.get_wishart_log_det_variance(wn, p)
	@test_approx_eq_eps(rel_diff(var(w_logdet_draws), wishart_log_det_var),
											0.0, 5e-2)

	wishart_suff_stats_cov = full(sparse_mat_from_tuples(
		get_wishart_sufficient_stats_variational_covariance(v0, wn, collect(1:p_ud),
			p_ud + 1, ud_ind)))

	@test_approx_eq_eps(rel_diff(wishart_suff_stats_sample_cov, wishart_suff_stats_cov),
	                    zeros(p_ud + 1, p_ud + 1), 1e-1)
end

function test_normal()
	srand(42)
	p = 2
	mu_mean = convert(Array{Float64}, collect(1:p))
	mu_sigma = 0.5 * (eye(p) + repmat([0.2], p, p))

	mu_dist = Distributions.MvNormal(mu_mean, mu_sigma)
	n_draws = 1_000_000
	mu_draws = rand(mu_dist, n_draws)'
	second_order_sample_cov =
		Float64[ cov(mu_draws[:, k1] .* mu_draws[:, k2],
		             mu_draws[:, k3] .* mu_draws[:, k4])
		         for k1=1:p, k2=1:p, k3=1:p, k4=1:p ]
	second_order_cov =
 		Float64[ ExponentialFamilies.get_mvn_fourth_order_cov(mu_mean, mu_sigma,
		                                                      k1, k2, k3, k4)
						 for k1=1:p, k2=1:p, k3=1:p, k4=1:p ]
	@test_approx_eq_eps(rel_diff(second_order_sample_cov, second_order_cov),
	                    zeros(p, p, p, p), 5e-2)

	suff_stat_draws = hcat(mu_draws,
	                       mu_draws[:, 1] .* mu_draws[:, 1],
												 mu_draws[:, 1] .* mu_draws[:, 2],
												 mu_draws[:, 2] .* mu_draws[:, 2])

  suff_stat_cov = full(sparse_mat_from_tuples(
		ExponentialFamilies.get_mvn_variational_covariance(
			mu_mean, mu_sigma, collect(1:2), Int64[ 3 4; 4 5])))
	suff_stat_sample_cov = cov(suff_stat_draws)
	@test_approx_eq_eps(rel_diff(suff_stat_sample_cov, suff_stat_cov),
											zeros(5, 5), 5e-2)

	suff_stat_draws_1d = hcat(mu_draws[:, 1], mu_draws[:, 1] .^ 2)
	suff_stat_1d_cov = full(sparse_mat_from_tuples(
		ExponentialFamilies.get_normal_variational_covariance(
			mu_mean[1], mu_mean[1] ^ 2 + mu_sigma[1, 1], 1, 2)))
	suff_stat_1d_sample_cov = cov(suff_stat_draws_1d)
	@test_approx_eq_eps(rel_diff(suff_stat_1d_sample_cov, suff_stat_1d_cov),
											zeros(2, 2), 5e-2)

end

function test_gamma()
	alpha = 5.0
	beta = 10.0
	gamma_dist = Distributions.Gamma(alpha, 1 / beta)

	n_draws = 100_000
	gamma_draws = rand(gamma_dist, n_draws)
	gamma_suff_draws = hcat(gamma_draws, log(gamma_draws))

	gamma_suff_cov = full(sparse_mat_from_tuples(
		ExponentialFamilies.get_gamma_variational_covariance(alpha, beta, 1, 2)))
	gamma_suff_sample_cov = cov(gamma_suff_draws)

	@test_approx_eq_eps(rel_diff(gamma_suff_sample_cov, gamma_suff_cov), zeros(2, 2), 1e-2)
	@test_approx_eq(ExponentialFamilies.gamma_entropy(alpha, beta),
	                Distributions.entropy(gamma_dist))
end

test_wishart()
test_normal()
test_gamma()


text:
using Test, Unitful, Plots
using Unitful: m, s, cm, DimensionError
using UnitfulRecipes

# Some helper functions to access the subplot labels and the series inside each test plot
xguide(plt, idx=length(plt.subplots)) = plt.subplots[idx].attr[:xaxis].plotattributes[:guide]
yguide(plt, idx=length(plt.subplots)) = plt.subplots[idx].attr[:yaxis].plotattributes[:guide]
zguide(plt, idx=length(plt.subplots)) = plt.subplots[idx].attr[:zaxis].plotattributes[:guide]
xseries(plt, idx=length(plt.series_list)) = plt.series_list[idx].plotattributes[:x]
yseries(plt, idx=length(plt.series_list)) = plt.series_list[idx].plotattributes[:y]
zseries(plt, idx=length(plt.series_list)) = plt.series_list[idx].plotattributes[:z]

macro isplot(ex) # @isplot macro to streamline tests
    :(@test $(esc(ex)) isa Plots.Plot)
end

@testset "heatmap" begin
   x = (1:3)m
   @isplot heatmap(x*x', clims=(1, 7)) # unitless
   @isplot heatmap(x*x', clims=(2m^2, 8m^2)) # units
   @isplot heatmap(x*x', clims=(2e6u"mm^2", 7e-6u"km^2")) # conversion
   @isplot heatmap(1:3, (1:3)m, x*x', clims=(1m^2, 7e-6u"km^2")) # mixed
end

@testset "plot(y)" begin
    y = rand(3)m

    @testset "no keyword argument" begin
        @test yguide(plot(y)) == "m"
        @test yseries(plot(y)) ≈ ustrip.(y)
    end

    @testset "ylabel" begin
        @test yguide(plot(y, ylabel="hello")) == "hello (m)"
        @test yguide(plot(y, ylabel=P"hello")) == "hello"
        @test yguide(plot(y, ylabel="")) == ""
        pl = plot(y; ylabel="hello")
        plot!(pl, -y)
        @test yguide(pl) == "hello (m)"
    end

    @testset "yunit" begin
        @test yguide(plot(y, yunit=cm)) == "cm"
        @test yseries(plot(y, yunit=cm)) ≈ ustrip.(cm, y)
    end

    @testset "ylims" begin # Using all(lims .≈ lims) because of uncontrolled type conversions?
        @test all(ylims(plot(y, ylims=(-1,3))) .≈ (-1,3))
        @test all(ylims(plot(y, ylims=(-1m,3m))) .≈ (-1,3))
        @test all(ylims(plot(y, ylims=(-100cm,300cm))) .≈ (-1,3))
        @test all(ylims(plot(y, ylims=(-100cm,3m))) .≈ (-1,3))
    end

    @testset "keyword combinations" begin
        @test yguide(plot(y, yunit=cm, ylabel="hello")) == "hello (cm)"
        @test yseries(plot(y, yunit=cm, ylabel="hello")) ≈ ustrip.(cm, y)
        @test all(ylims(plot(y, yunit=cm, ylims=(-1,3))) .≈ (-1,3))
        @test all(ylims(plot(y, yunit=cm, ylims=(-1,3))) .≈ (-1,3))
        @test all(ylims(plot(y, yunit=cm, ylims=(-100cm,300cm))) .≈ (-100,300))
        @test all(ylims(plot(y, yunit=cm, ylims=(-100cm,3m))) .≈ (-100,300))
    end
end

@testset "plot(x,y)" begin
    x, y = randn(3)m, randn(3)s

    @testset "no keyword argument" begin
        @test xguide(plot(x,y)) == "m"
        @test xseries(plot(x,y)) ≈ ustrip.(x)
        @test yguide(plot(x,y)) == "s"
        @test yseries(plot(x,y)) ≈ ustrip.(y)
    end

    @testset "labels" begin
        @test xguide(plot(x, y, xlabel= "hello")) == "hello (m)"
        @test xguide(plot(x, y, xlabel=P"hello")) == "hello"
        @test yguide(plot(x, y, ylabel= "hello")) == "hello (s)"
        @test yguide(plot(x, y, ylabel=P"hello")) == "hello"
        @test xguide(plot(x, y, xlabel= "hello", ylabel= "hello")) == "hello (m)"
        @test xguide(plot(x, y, xlabel=P"hello", ylabel=P"hello")) == "hello"
        @test yguide(plot(x, y, xlabel= "hello", ylabel= "hello")) == "hello (s)"
        @test yguide(plot(x, y, xlabel=P"hello", ylabel=P"hello")) == "hello"
    end

    @testset "unitformat" begin
        args = (x, y)
        kwargs = (:xlabel=>"hello", :ylabel=>"hello")
        @test yguide(plot(args...; kwargs..., unitformat=nothing)) == "hello s"
        @test yguide(plot(args...; kwargs..., unitformat=(l, u) -> string(u, " is the unit of ", l))) == "s is the unit of hello"
        @test yguide(plot(args...; kwargs..., unitformat=", dear ")) == "hello, dear s"
        @test yguide(plot(args...; kwargs..., unitformat=(", dear ", " esq."))) == "hello, dear s esq."
        @test yguide(plot(args...; kwargs..., unitformat=("well ", ", dear ", " esq."))) == "well hello, dear s esq."
        @test yguide(plot(args...; kwargs..., unitformat='?')) == "hello ? s"
        @test yguide(plot(args...; kwargs..., unitformat=('<', '>'))) == "hello <s>"
        @test yguide(plot(args...; kwargs..., unitformat=('A', 'B', 'C'))) == "Ahello BsC"
        @test yguide(plot(args...; kwargs..., unitformat=false)) == "hello s"
        @test yguide(plot(args...; kwargs..., unitformat=true)) == "hello (s)"
        @test yguide(plot(args...; kwargs..., unitformat=:round)) == "hello (s)"
        @test yguide(plot(args...; kwargs..., unitformat=:square)) == "hello [s]"
        @test yguide(plot(args...; kwargs..., unitformat=:curly)) == "hello {s}"
        @test yguide(plot(args...; kwargs..., unitformat=:angle)) == "hello <s>"
        @test yguide(plot(args...; kwargs..., unitformat=:slash)) == "hello / s"
        @test yguide(plot(args...; kwargs..., unitformat=:slashround)) == "hello / (s)"
        @test yguide(plot(args...; kwargs..., unitformat=:slashsquare)) == "hello / [s]"
        @test yguide(plot(args...; kwargs..., unitformat=:slashcurly)) == "hello / {s}"
        @test yguide(plot(args...; kwargs..., unitformat=:slashangle)) == "hello / <s>"
        @test yguide(plot(args...; kwargs..., unitformat=:verbose)) == "hello in units of s"
    end
end

@testset "With functions" begin
    x, y = randn(3), randn(3)
    @testset "plot(f, x) / plot(x, f)" begin
        f(x) = x^2
        @test plot(  f, x*m) isa Plots.Plot
        @test plot(x*m,   f) isa Plots.Plot
        g(x) = x*m # If the unit comes from the function only then it throws
        @test_throws DimensionError plot(x, g) isa Plots.Plot
        @test_throws DimensionError plot(g, x) isa Plots.Plot
    end
    @testset "plot(x, y, f)" begin
        f(x,y) = x*y
        @test plot(x*m, y*s, f) isa Plots.Plot
        @test plot(x*m,   y, f) isa Plots.Plot
        @test plot(  x, y*s, f) isa Plots.Plot
        g(x,y) = x*y*m # If the unit comes from the function only then it throws
        @test_throws DimensionError plot(x, y, g) isa Plots.Plot
    end
    @testset "plot(f, u)" begin
        f(x) = x^2
        pl = plot(x*m, f.(x*m))
        @test plot!(pl, f, m) isa Plots.Plot
        @test_throws DimensionError plot!(pl, f, s) isa Plots.Plot
        pl = plot(f, m)
        @test xguide(pl) == string(m)
        @test yguide(pl) == string(m^2)
        g(x) = exp(x/(3m))
        @test plot(g, u"m") isa Plots.Plot
    end
end

@testset "More plots" begin
    @testset "data as $dtype" for dtype in [:Vectors, :Matrices, Symbol("Vectors of vectors")]
        if dtype == :Vectors
            x, y, z = randn(10), randn(10), randn(10)
        elseif dtype == :Matrices
            x, y, z = randn(10,2), randn(10,2), randn(10,2)
        else
            x, y, z = [rand(10), rand(20)], [rand(10), rand(20)], [rand(10), rand(20)]
        end

        @testset "One array" begin
            @test plot(x*m)                    isa Plots.Plot
            @test plot(x*m, ylabel="x")        isa Plots.Plot
            @test plot(x*m, ylims=(-1,1))      isa Plots.Plot
            @test plot(x*m, ylims=(-1,1) .* m) isa Plots.Plot
            @test plot(x*m, yunit=u"km")       isa Plots.Plot
            @test plot(x*m, xticks=(1:3)*m)    isa Plots.Plot
        end

        @testset "Two arrays" begin
            @test plot(x*m, y*s)                    isa Plots.Plot
            @test plot(x*m, y*s, xlabel="x")        isa Plots.Plot
            @test plot(x*m, y*s, xlims=(-1,1))      isa Plots.Plot
            @test plot(x*m, y*s, xlims=(-1,1) .* m) isa Plots.Plot
            @test plot(x*m, y*s, xunit=u"km")       isa Plots.Plot
            @test plot(x*m, y*s, ylabel="y")        isa Plots.Plot
            @test plot(x*m, y*s, ylims=(-1,1))      isa Plots.Plot
            @test plot(x*m, y*s, ylims=(-1,1) .* s) isa Plots.Plot
            @test plot(x*m, y*s, yunit=u"ks")       isa Plots.Plot
            @test plot(x*m, y*s, yticks=(1:3)*s)    isa Plots.Plot
            @test scatter(x*m, y*s)                 isa Plots.Plot
            if dtype ≠ Symbol("Vectors of vectors")
                @test scatter(x*m, y*s, zcolor=z*(m/s)) isa Plots.Plot
            end
        end

        @testset "Three arrays" begin
            @test plot(x*m, y*s, z*(m/s))                        isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), xlabel="x")            isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), xlims=(-1,1))          isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), xlims=(-1,1) .* m)     isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), xunit=u"km")           isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), ylabel="y")            isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), ylims=(-1,1))          isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), ylims=(-1,1) .* s)     isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), yunit=u"ks")           isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), zlabel="z")            isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), zlims=(-1,1))          isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), zlims=(-1,1) .* (m/s)) isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), zunit=u"km/hr")        isa Plots.Plot
            @test plot(x*m, y*s, z*(m/s), zticks=(1:2)*m/s)      isa Plots.Plot
            @test scatter(x*m, y*s, z*(m/s))                     isa Plots.Plot
        end

        @testset "Unitful/unitless combinations" begin
            mystr(x::Array{<:Quantity}) = "Q"
            mystr(x::Array) = "A"
            @testset "plot($(mystr(xs)), $(mystr(ys)))" for xs in [x, x*m], ys in [y, y*s]
                @test plot(xs, ys) isa Plots.Plot
            end
            @testset "plot($(mystr(xs)), $(mystr(ys)), $(mystr(zs)))" for xs in [x, x*m], ys in [y, y*s], zs in [z, z*(m/s)]
                @test plot(xs, ys, zs) isa Plots.Plot
            end
        end
    end

    @testset "scatter(x::$(us[1]), y::$(us[2]))" for us in collect(Iterators.product(fill([1, u"m", u"s"], 2)...))
        x, y = rand(10)*us[1], rand(10)*us[2]
        @test scatter(x,y)  isa Plots.Plot
        @test scatter(x,y, markersize=x)  isa Plots.Plot
        @test scatter(x,y, line_z=x)  isa Plots.Plot
    end

    @testset "contour(x::$(us[1]), y::$(us[2]))" for us in collect(Iterators.product(fill([1, u"m", u"s"], 2)...))
        x, y = (1:0.01:2)*us[1], (1:0.02:2)*us[2]
        z = x' ./ y
        @test contour(x,y,z)  isa Plots.Plot
        @test contourf(x,y,z) isa Plots.Plot
    end

    @testset "ProtectedString" begin
        y = rand(10)*u"m"
        @test plot(y, label=P"meters") isa Plots.Plot
    end
end

@testset "Comparing apples and oranges" begin
    x1 = rand(10) * u"m"
    x2 = rand(10) * u"cm"
    x3 = rand(10) * u"s"
    plt = plot(x1)
    plt = plot!(plt, x2)
    @test yguide(plt) == "m"
    @test yseries(plt) ≈ ustrip.(x2) / 100
    @test_throws DimensionError plot!(plt, x3) # can't place seconds on top of meters!
end

@testset "Bare units" begin
    plt = plot(u"m", u"s")
    @test xguide(plt) == "m"
    @test yguide(plt) == "s"
    @test iszero(length(plt.series_list[1].plotattributes[:y]))
    hline!(plt, [1u"hr"])
    @test yguide(plt) == "s"
end

@testset "Inset subplots" begin
    x1 = rand(10) * u"m"
    x2 = rand(10) * u"s"
    plt = plot(x1)
    plt = plot!(x2, inset=bbox(0.5, 0.5, 0.3, 0.3), subplot=2)
    @test yguide(plt,1) == "m"
    @test yguide(plt,2) == "s"
end

@testset "Missing values" begin
    x = 1:5
    y = [1.0*u"s", 2.0*u"s", missing, missing, missing]
    @show typeof(y)
    plt = plot(x,y)
    @test yguide(plt,1) == "s"
end

@testset "Errors" begin
    x = rand(10) * u"mm"
    ex = rand(10) * u"μm"
    y = rand(10) * u"s"
    ey = rand(10) * u"ms"
    plt =  plot(x, y, xerr=ex, yerr=ey)
    @test plt isa Plots.Plot
    @test xguide(plt) == "mm"
    @test yguide(plt) == "s"
end

@testset "Ribbon" begin
    x = rand(10) * u"mm"
    y = rand(10) * u"s"
    ribbon = rand(10) * u"ms"
    plt =  plot(x, y, ribbon=ribbon)
    @test plt isa Plots.Plot
    @test xguide(plt) == "mm"
    @test yguide(plt) == "s"
end

@testset "Fillrange" begin
    x = rand(10) * u"mm"
    y = rand(10) * u"s"
    fillrange = rand(10) * u"ms"
    plt =  plot(x, y, fillrange=fillrange)
    @test plt isa Plots.Plot
    @test xguide(plt) == "mm"
    @test yguide(plt) == "s"
end

# https://github.com/jw3126/UnitfulRecipes.jl/issues/60
@testset "Start with empty plot" begin
    plt = plot()
    plot!(plt, (1:3)m)
    @test yguide(plt) == "m"
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
a = [1.0, 2.0]
b = [1.0, 2.0]
c = zeros(2)

f!(a, b, c) = @. c = b - a/b^2
@time f!(a, b, c)
@time f!(a, b, c)
@time f!(a, b, c)

# %%
a = [1.0, 2.0]
b = [1.0, 2.0]
c = zeros(2)

f2!(a, b, c) = @. c = -a/b^2 + b
@time f2!(a, b, c)
@time f2!(a, b, c)
@time f2!(a, b, c)

# %%
using BenchmarkTools
@btime f!($a, $b, $c)
@btime f2!($a, $b, $c)

# %%
g!(a, b, c) = c .= .- a ./ (b .^ 2) .+ b
@btime g!($a, $b, $c)

# %%

max_stars_repo_path: 0010/weired behavior of broadcasting.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0010/weired behavior of broadcasting.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0010/weired behavior of broadcasting.jl

max_forks_repo_name: genkuroki/public


text:
using Distributions
using Klara

ρ = Hyperparameter(:ρ, 1)

p1 = BasicContUnvParameter(
  :p1,
  setpdf=(p::Float64, v::Vector) -> Normal(v[1]*v[3], sqrt(1-abs2(v[1]))),
  nkeys=3
)

p2 = BasicContUnvParameter(
  :p2,
  setpdf=(p::Float64, v::Vector) -> Normal(v[1]*v[2], sqrt(1-abs2(v[1]))),
  nkeys=3
)

model = GenericModel([ρ, p1, p2], isindexed=false)

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:ρ=>0.8, :p1=>5.1, :p2=>2.3)

job = BasicGibbsJob(model, Dict(), mcrange, v0)

@time run(job)

output(job)

chains = Dict(job)

[mean(chains[k]) for k in keys(chains)]

cor(chains[:p1].value, chains[:p2].value)


text:
using Flux
using Flux: logitcrossentropy, normalise, onecold, onehotbatch
using Statistics: mean
using MLFlowLogger
using Logging
using Parameters: @with_kw

@with_kw mutable struct Args
    lr::Float64 = 0.5
    repeat::Int = 110
end

function get_processed_data(args)
    labels = Flux.Data.Iris.labels()
    features = Flux.Data.Iris.features()

    # Subract mean, divide by std dev for normed mean of 0 and std dev of 1.
    normed_features = normalise(features, dims=2)

    klasses = sort(unique(labels))
    onehot_labels = onehotbatch(labels, klasses)

    # Split into training and test sets, 2/3 for training, 1/3 for test.
    train_indices = [1:3:150 ; 2:3:150]

    X_train = normed_features[:, train_indices]
    y_train = onehot_labels[:, train_indices]

    X_test = normed_features[:, 3:3:150]
    y_test = onehot_labels[:, 3:3:150]

    #repeat the data `args.repeat` times
    train_data = Iterators.repeated((X_train, y_train), args.repeat)
    test_data = (X_test,y_test)

    return train_data, test_data
end

# Accuracy Function
accuracy(x, y, model) = mean(onecold(model(x)) .== onecold(y))

# Function to build confusion matrix
function confusion_matrix(X, y, model)
    ŷ = onehotbatch(onecold(model(X)), 1:3)
    y * transpose(ŷ)
end

function train(; kws...)
    # Initialize hyperparameter arguments
    args = Args(; kws...)

    # Setting up the MLFlowLogger
    logger = MLFLogger()

    #Loading processed data
    train_data, test_data = get_processed_data(args)

    # Declare model taking 4 features as inputs and outputting 3 probabiltiies,
    # one for each species of iris.
    model = Chain(Dense(4, 3))

    # Defining loss function to be used in training
    # For numerical stability, we use here logitcrossentropy
    loss(x, y) = logitcrossentropy(model(x), y)

    # Training
    # Gradient descent optimiser with learning rate `args.lr`
    optimiser = Descent(args.lr)

    function MLFCallback()
        with_logger(logger) do
          @info "train" loss=loss(first(iterate(train_data))...) acc=accuracy(first(iterate(train_data))..., model) log_step_increment=0
          @info "test" loss=loss(test_data...) acc=accuracy(test_data..., model)
        end
    end

    println("Starting training.")
    Flux.@epochs 15 Flux.train!(loss, params(model), train_data, optimiser, cb = Flux.throttle(MLFCallback, 5))

    return model, test_data
end

function test(model, test)
    # Testing model performance on test data
    X_test, y_test = test
    accuracy_score = accuracy(X_test, y_test, model)

    println("\nAccuracy: $accuracy_score")

    # Sanity check.
    @assert accuracy_score > 0.8

    # To avoid confusion, here is the definition of a Confusion Matrix: https://en.wikipedia.org/wiki/Confusion_matrix
    println("\nConfusion Matrix:\n")
    display(confusion_matrix(X_test, y_test, model))
end

cd(@__DIR__)
model, test_data = train()
test(model, test_data)

max_stars_repo_name: denniseik/azureml-examples

max_issues_repo_name: denniseik/azureml-examples

max_forks_repo_name: denniseik/azureml-examples


text:
import Base: Pair, Tuple, show, ==

struct Edge{T<:Integer}
    src::T
    dst::T
end

Edge(t::Tuple) = Edge(t[1], t[2])
Edge(p::Pair) = Edge(p.first, p.second)

# Accessors
src(e::Edge) = e.src
dst(e::Edge) = e.dst

# I/O
show(io::IO, e::Edge) = print(io, "Edge: $(e.src) => $(e.dst)")

==(e1::Edge, e2::Edge) = (src(e1) == src(e2) && dst(e1) == dst(e2))


text:
"""
    Chain(layers...)

Chain multiple layers / functions together, so that they are called in sequence
on a given input.

```julia
m = Chain(x -> x^2, x -> x+1)
m(5) == 26

m = Chain(Dense(10, 5), Dense(5, 2))
x = rand(10)
m(x) == m[2](m[1](x))
```

`Chain` also supports indexing and slicing, e.g. `m[2]` or `m[1:end-1]`.
`m[1:3](x)` will calculate the output of the first three layers.
"""
type Chain
  layers::Vector{Any}
  Chain(xs...) = new([xs...])
end

@forward Chain.layers Base.getindex, Base.first, Base.last, Base.endof, Base.push!
@forward Chain.layers Base.start, Base.next, Base.done

children(c::Chain) = c.layers
mapchildren(f, c::Chain) = Chain(f.(c.layers)...)

(c::Chain)(x) = foldl((x, m) -> m(x), x, c.layers)

Base.getindex(c::Chain, i::AbstractArray) = Chain(c.layers[i]...)

function Base.show(io::IO, c::Chain)
  print(io, "Chain(")
  join(io, c.layers, ", ")
  print(io, ")")
end

"""
    Dense(in::Integer, out::Integer, σ = identity)

Creates a traditional `Dense` layer with parameters `W` and `b`.

    y = σ.(W * x .+ b)

The input `x` must be a vector of length `in`, or a batch of vectors represented
as an `in × N` matrix. The out `y` will be a vector or batch of length `out`.

```julia
julia> d = Dense(5, 2)
Dense(5, 2)

julia> d(rand(5))
Tracked 2-element Array{Float64,1}:
  0.00257447
  -0.00449443
```
"""
struct Dense{F,S,T}
  σ::F
  W::S
  b::T
end

Dense(in::Integer, out::Integer, σ = identity; init = initn) =
  Dense(σ, param(init(out, in)), param(init(out)))

treelike(Dense)

function (a::Dense)(x)
  W, b, σ = a.W, a.b, a.σ
  σ.(W*x .+ b)
end

function Base.show(io::IO, l::Dense)
  print(io, "Dense(", size(l.W, 2), ", ", size(l.W, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end


text:
isdefined(:LazyAlgebra) || include("../src/LazyAlgebra.jl")

module LazyAlgebraBenchmarks

using BenchmarkTools
using LazyAlgebra
import LazyAlgebra: vdot

vdot(::Type{Val{:BLAS}}, x, y) =
    LazyAlgebra.blas_vdot(x,y)

vdot(::Type{Val{:Julia}}, x, y) =
    dot(reshape(x, length(x)), reshape(y, length(y)))

function vdot(::Type{Val{:basic}},
              x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tx, Ty, N}
    @assert size(x) == size(y)
    T = promote_type(Tx, Ty)
    s = zero(T)
    for i in eachindex(x, y)
        s += convert(T, x[i])*convert(T, y[i])
    end
    return s
end

function vdot(::Type{Val{:inbounds}},
              x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tx, Ty, N}
    @assert size(x) == size(y)
    T = promote_type(Tx, Ty)
    s = zero(T)
    @inbounds for i in eachindex(x, y)
        s += convert(T, x[i])*convert(T, y[i])
    end
    return s
end

function vdot(::Type{Val{:simd}},
              x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tx, Ty, N}
    @assert size(x) == size(y)
    T = promote_type(Tx, Ty)
    s = zero(T)
    @inbounds @simd for i in eachindex(x, y)
        s += convert(T, x[i])*convert(T, y[i])
    end
    return s
end

function vdot(::Type{Val{:simdlinear}},
              x::DenseArray{Tx,N},
              y::DenseArray{Ty,N}) where {Tx, Ty, N}
    @assert size(x) == size(y)
    T = promote_type(Tx, Ty)
    s = zero(T)
    @inbounds @simd for i in 1:length(x)
        s += x[i]*y[i]
    end
    return s
end

function vdot(::Type{Val{:basic}},
              w::AbstractArray{Tw,N},
              x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tw, Tx, Ty, N}
    @assert size(x) == size(y)
    T = promote_type(Tw, Tx, Ty)
    s = zero(T)
    for i in eachindex(x, x, y)
        s += convert(T, w[i])*convert(T, x[i])*convert(T, y[i])
    end
    return s
end

function vdot(::Type{Val{:inbounds}},
              w::AbstractArray{Tw,N},
              x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tw, Tx, Ty, N}
    @assert size(w) == size(x) == size(y)
    T = promote_type(Tw, Tx, Ty)
    s = zero(T)
    @inbounds for i in eachindex(w, x, y)
        s += convert(T, w[i])*convert(T, x[i])*convert(T, y[i])
    end
    return s
end

function vdot(::Type{Val{:simd}},
              w::AbstractArray{Tw,N},
              x::AbstractArray{Tx,N},
              y::AbstractArray{Ty,N}) where {Tw, Tx, Ty, N}
    @assert size(w) == size(x) == size(y)
    T = promote_type(Tw, Tx, Ty)
    s = zero(T)
    @inbounds @simd for i in eachindex(w, x, y)
        s += convert(T, w[i])*convert(T, x[i])*convert(T, y[i])
    end
    return s
end

function vdot(::Type{Val{:simdlinear}},
              w::DenseArray{Tw,N},
              x::DenseArray{Tx,N},
              y::DenseArray{Ty,N}) where {Tw, Tx, Ty, N}
    @assert size(w) == size(x) == size(y)
    T = promote_type(Tw, Tx, Ty)
    s = zero(T)
    @inbounds @simd for i in 1:length(x)
        s += w[i]*x[i]*y[i]
    end
    return s
end

testdot() = testdot(33, 33)

function testdot(_dims::Integer...)
    global dims, w, x, y, z

    dims = Int.(_dims)

    #show(stdout, MIME"text/plain"(), @benchmark $p(dat, a, img, b, rois))

    println("\n\nDot products of $dims elements")

    println("\\begin{tabular}{lrr}")
    println("\\hline")
    println(" & \\multicolumn{2}{c}{Median time (ns) for ",
            dims, " elements} \\\\")
    println("Operation & \\texttt{Float32} & \\texttt{Float64}\\\\")
    println("\\hline")
    println("\\hline")
    for p in (:Julia, :BLAS, :basic, :inbounds, :simd, :simdlinear)
        @printf(" dot %-12s", string("(",p,")"))
        for T in (Float32, Float64)
            w = randn(T, dims)
            x = randn(T, dims)
            y = randn(T, dims)
            z = randn(T, dims)
            s = vdot(Val{p}, x, y)
            t = @benchmark vdot($(Val{p}), x, y)
            @printf(" & %4.0f", median(t.times))
        end
        println(" \\\\")
    end
    println("\\hline")
    for p in (:basic, :inbounds, :simd, :simdlinear)
        @printf("wdot %-12s", string("(",p,")"))
        for T in (Float32, Float64)
            w = randn(T, dims)
            x = randn(T, dims)
            y = randn(T, dims)
            z = randn(T, dims)
            s = vdot(Val{p}, x, y)
            t = @benchmark vdot($(Val{p}), w, x, y)
            @printf(" & %4.0f", median(t.times))
        end
        println(" \\\\")
    end
    println("\\hline")
    println("\\end{tabular}")
end

norminf(x::AbstractArray) = ((xmn, xmx) = extrema(x); return max(-xmn, xmx))
norm1(x::AbstractArray) = sum(abs.(x))
norm2(x::AbstractArray) = sum(abs2.(x))

end


text:
for func in (:daily, :weekly, :monthly)
    x = "digital_currency_$(func)"
    fname = Symbol(x)
    @eval begin
        function ($fname)(symbol::String, market::String="USD"; client = GLOBAL[], outputsize::String="compact", datatype::Union{String, Nothing}=nothing, parser = "default")
            @argcheck in(outputsize, ["compact", "full"])
            @argcheck in(datatype, ["json", "csv", nothing])
            params = Dict(
                "function"=>uppercase($x),
                "symbol"=>symbol,
                "market"=>market,
                "outputsize"=>outputsize,
                "datatype"=>isnothing(datatype) ? "csv" : datatype,
                "apikey"=>key(client)
            )
            uri = _build_uri(client.scheme, client.host, "query", params)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end
        export $fname
    end
end

function digital_currency_intraday(symbol::String, market::String="USD", interval::String="1min"; client = GLOBAL[], outputsize::String="compact", datatype::Union{String, Nothing}=nothing, parser = "default")
    @argcheck in(outputsize, ["compact", "full"])
    @argcheck in(datatype, ["json", "csv", nothing])
    @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min"])
    params = Dict(
        "function"=>"CRYPTO_INTRADAY",
        "symbol"=>symbol,
        "market"=>market,
        "interval"=>interval,
        "outputsize"=>outputsize,
        "datatype"=>isnothing(datatype) ? "csv" : datatype,
        "apikey"=>key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

function crypto_rating(symbol::String; client = GLOBAL[], parser = "default")
    params = Dict(
        "function"=>"CRYPTO_RATING",
        "symbol"=>symbol,
        "apikey"=>key(client)
    )
    uri = _build_uri(client.scheme, client.host, "query", params)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, "json")
    return p(data)
end


text:
"""
`module  JAC.DoubleAuger`  ... a submodel of JAC that contains all methods for computing double Auger and autoionization amplitudes
                               and rates; it is using JAC, JAC.ManyElectron, JAC.ImpactExcitation, JAC.AutoIonization.
"""
module DoubleAuger

    using JAC, JAC.ManyElectron, JAC.ImpactExcitation, JAC.AutoIonization

    """
    `struct  DoubleAuger.Settings`  ... defines a type for the settings in estimating double-Auger and autoionization rates.

        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
        + minAugerEnergy          ::Float64                      ... Minimum energy of free (Auger) electrons to be included.
        + maxAugerEnergy          ::Float64                      ... Maximum energy of free (Auger) electrons to be included.
        + maxKappa                ::Int6464                      ... Maximum kappa value of partial waves to be included.
    """
    struct Settings
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
        minAugerEnergy            ::Float64
        maxAugerEnergy            ::Float64
        maxKappa                  ::Int64
    end

    """
    `JAC.DoubleAuger.Settings()`  ... constructor for the default values of DoubleAuger line computations
    """
    function Settings()
        Settings(false, false, Array{Tuple{Int64,Int64},1}[], 0., 10e5, 100)
    end

    # `Base.show(io::IO, settings::DoubleAuger.Settings)`  ... prepares a proper printout of the variable settings::DoubleAuger.Settings.
    function Base.show(io::IO, settings::DoubleAuger.Settings)
        println(io, "printBeforeComputation:       $(settings.printBeforeComputation)  ")
        println(io, "selectLines:                  $(settings.selectLines)  ")
        println(io, "selectedLines:                $(settings.selectedLines)  ")
        println(io, "minAugerEnergy:               $(settings.minAugerEnergy)  ")
        println(io, "maxAugerEnergy:               $(settings.maxAugerEnergy)  ")
        println(io, "maxKappa:                     $(settings.maxKappa)  ")
    end

    """
    `struct  Channel`  ... defines a type for a DoubleAuger channel to help characterize a scattering (continuum) state of many
                           electron-states with two free electrons.

        + kappa1         ::Int64                ... partial-wave of the free electron_1
        + kappa2         ::Int64                ... partial-wave of the free electron_2
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase1         ::Float64              ... phase of the partial wave_1
        + phase2         ::Float64              ... phase of the partial wave_2
        + amplitude      ::Complex{Float64}     ... Auger amplitude associated with the given channel.
    """
    struct  Channel
        kappa1           ::Int64
        kappa2           ::Int64
        symmetry         ::LevelSymmetry
        phase1           ::Float64
        phase2           ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  Line`  ... defines a type for a Auger line that may include the definition of sublines and their corresponding amplitudes.

        + initialLevel   ::Level         ... initial-(state) level
        + finalLevel     ::Level         ... final-(state) level
        + electronEnergy ::Float64       ... Energy of the (incoming free) electron.
        + totalRate      ::Float64       ... Total rate of this line.
        + hasChannels    ::Bool          ... Determines whether the individual scattering (sub-) channels are defined in terms of their free-
                                             electron energies, kappas and the total angular momentum/parity as well as the amplitude, or not.
        + channels       ::Array{DoubleAuger.Channel,1}  ... List of DoubleAuger channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        electronEnergy   ::Float64
        totalRate        ::Float64
        hasChannels      ::Bool
        channels         ::Array{DoubleAuger.Channel,1}
    end

    # `Base.show(io::IO, line::DoubleAuger.Line)`  ... prepares a proper printout of the variable line::DoubleAuger.Line.
    function Base.show(io::IO, line::DoubleAuger.Line)
        println(io, "initialLevel:           $(line.initialLevel)  ")
        println(io, "finalLevel:             $(line.finalLevel)  ")
        println(io, "electronEnergy:         $(line.electronEnergy)  ")
        println(io, "totalRate:              $(line.totalRate)  ")
        println(io, "hasChannels:            $(line.hasChannels)  ")
        println(io, "channels:               $(line.channels)  ")
    end

end # module


text:
using FourVectors
using Test

const mπ0 = 0.135
#
pπ0 = Particle(1,1,30; msq = mπ0^2)
#
pγ1 = mπ0/2 * Particle(E=1,p=[sin(0.3),0,cos(0.3)])
Bz!(pγ1,boostfactor(pπ0))
cosθ,ϕ = sphericalangles(pπ0)
Ry!(pγ1,acos(cosθ))
Rz!(pγ1,ϕ)
#
pγ2 = pπ0 - pγ1
#
invmasssq(pγ1)
@test invmasssq(pγ1) < 1e-10
@test invmasssq(pγ2) < 1e-10
@test invmasssq(pγ1+pγ2) == invmasssq(pπ0)
@test pγ1+pγ2 == pπ0
#
p = Particle(160,3,1; msq = 1.1);
α = 1.4; γ = 10;
@test sum(Rx(Rx(p,α),-α) .≈ p) == 4
@test sum(Rz(Rz(p,α),-α) .≈ p) == 4
@test sum(Ry(Ry(p,α),-α) .≈ p) == 4

@test sum(Bz(Bz(p,-γ),γ) .≈ p) == 4


text:
inverse(side::HorizonSide) = HorizonSide(mod(Int(side) + 2, 4))

left(side::HorizonSide) = HorizonSide(mod(Int(side) + 1, 4))

right(side::HorizonSide) = HorizonSide(mod(Int(side) - 1, 4))

function Back!(r::Robot, side::HorizonSide, num_steps::Int)
    while num_steps > 0
        move_if_possible!(r, side)
        num_steps -= 1
    end
end

function movements!(r::Robot, side::HorizonSide)
    while !isborder(r, side)
        move!(r, side)
    end
end

function putmarkers!(r::Robot, side::HorizonSide)
    while !isborder(r, side)
        move!(r, side)
        putmarker!(r)
    end
end

function movements!(r::Robot, side::HorizonSide, num_steps::Int)
    for i in 1:num_steps
        move!(r, side)
    end
end

# on particular digits of steps move to the.....
function digit_movements!(r::Robot, side::HorizonSide)
    num_steps = 0
    while !isborder(r, side)
        move!(r, side)
        num_steps += 1
    end
    return num_steps
end

function digit_movements!(r::Robot, side::HorizonSide,num_steps::Int)
    num_steps = 0
    while !isborder(r, side)
        move!(r, side)
        num_steps += 1
    end
    return num_steps
end

function putmarkers!(r::Robot, side::HorizonSide,num::Int)
    for _ in 1:num
        putmarker!(r)
        move!(r,side)
    end
    putmarker!(r)
end

# Перемещает робота в заданном направлении, если это возможно, и возвращает true,
# если перемещение состоялось; в противном случае - false.

function move_if_possible!(r::Robot, now_side::HorizonSide)::Bool
    # в направлении now_side не встретилось перегородки
    if !isborder(r, now_side)
        move!(r, now_side)
        return true
    end

    # в направлении now_side есть перегородка => начинается ее обход
    left_side = left(now_side)
    right_side = inverse(left_side)
    num_steps = 0

    while isborder(r, now_side)
        #print(num_steps)
        if !isborder(r, left_side)
            move!(r, left_side)
            num_steps += 1
        else
            movements!(r, right_side, num_steps)
            return false # робот пришел в угол
        end
    end

    move!(r, now_side)

    while isborder(r, right_side)
        if !isborder(r, now_side)
            move!(r, now_side)
        else
            movements!(r, left_side, num_steps)
            return false # робот пришел в угол
        end
    end

    for i in 1:num_steps
        move!(r, right_side)
    end

    return true
end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
module QuantileSummaries

export QuantileStat, QuantileSummary, DataItem,
    fit!,
    value,
    nobs,
    summarize,
    percentile,
    qindex,
    qvalue,
    ndata

import StatsBase: nobs, fit!, merge!, value
import OnlineStatsBase: value, OnlineStat, _fit!, _merge!

include("QuantileStat.jl")

end # module


text:
import MLJBase
import MLJModelInterface
import Soss
import Statistics
import NamedTupleTools

function MLJModelInterface.fit(sm::SossMLJModel, verbosity::Integer, X, y, w = nothing)
    # construct the model
    args = merge(sm.transform(X), sm.hyperparams)

    jd = sm.model(args)

    y_namedtuple = NamedTupleTools.namedtuple(sm.response)(tuple(y))

    post = sm.infer(jd, y_namedtuple) # this is the slow part: doing inference with e.g. MCMC

    # TODO: Allow w to be included

    cache = nothing
    report = NamedTuple{}()

    newargs = setdiff(Soss.sampled(jd.model),(:y,))
    pred = Soss.predictive(jd.model, newargs...)
    fitresult = (model=sm, post=post)
    return (fitresult, cache, report)
end

function MLJModelInterface.clean!(smm::SossMLJModel)
    warning = ""
    return warning
end

function MLJModelInterface.predict(sm::SossMLJModel{<:SossMLJPredictor}, fitresult, Xnew)
    m = sm.model
    post = fitresult.post
    pred = Soss.predictive(m, keys(post[1])...)

    return map(Tables.rowtable(Xnew)) do xrow
        args = merge(sm.transform([xrow]), sm.hyperparams)
        return SossMLJPredictor(sm, post, pred, args)
    end
end

# This method MUST return a `UnivariateFiniteVector`
# Currently, this is a hacky and incorrect implementation.
# TODO: Implement this method correctly.
function MLJModelInterface.predict(sm::SossMLJModel{<:MLJBase.UnivariateFinite}, fitresult, Xnew; response = sm.response)
    return _predict_marginals_incorrectly(sm, fitresult, Xnew; response = sm.response)
end

# TODO: remove this hacky and incorrect method.
# Here's how it works:
# 1. Construct the joint distribution
# 2. Draw five thousand samples from the joint distribution
# 3. Use these five thousand samples to empirically estimate the probabilities of each class
# 4. Return a `UnivariateFiniteVector` using these probabilities
function _predict_marginals_incorrectly(sm::SossMLJModel{<:MLJBase.UnivariateFinite}, fitresult, Xnew; response = sm.response)
    predictor_joint = MLJModelInterface.predict_joint(sm, fitresult, Xnew)
    num_samples = 5_000
    samples = hcat([rand(predictor_joint; response = response) for sample = 1:num_samples]...)
    pool = sm.hyperparams.pool
    levels = pool.levels
    num_levels = length(levels)
    num_rows = size(Xnew, 1)
    probs = Matrix{Float64}(undef, num_rows, num_levels)
    for i = 1:num_rows
        for j = 1:num_levels
            level = levels[j]
            level_prob = sum(samples[i, :] .== level)/num_samples
            probs[i, j] = level_prob
        end
    end
    support = levels
    marginal_distributions = MLJBase.UnivariateFinite(support, probs; pool = pool)
    return marginal_distributions
end

function MLJModelInterface.predict_joint(sm::SossMLJModel, fitresult, Xnew)
    m = sm.model
    post = fitresult.post
    pred = Soss.predictive(m, keys(post[1])...)
    args = merge(sm.transform(Xnew), sm.hyperparams)
    return SossMLJPredictor(sm, post, pred, args)
end

function MLJModelInterface.predict_mean(sm::SossMLJModel{<:SossMLJPredictor}, fitresult, Xnew; response = sm.response)
    predictor_joint = MLJModelInterface.predict_joint(sm, fitresult, Xnew)
    return Statistics.mean(getproperty(predict_particles(predictor_joint, Xnew), response))
end

function MLJModelInterface.predict_mode(sm::SossMLJModel{<:MLJBase.UnivariateFinite}, fitresult, Xnew; response = sm.response)
    marginal_distributions = MLJModelInterface.predict(sm, fitresult, Xnew; response = sm.response) # `marginal_distributions` is of type `UnivariateFiniteVector`
    modes = [MLJBase.mode(marginal_dist) for marginal_dist in marginal_distributions]
    return modes
end


text:
module HOI4CombatSimulator

export get_unit_list, Unit, DivisionTemplate, Division, fire, fire!

include("defs.jl")
include("get_data.jl")
include("data_model.jl")
include("resolve.jl")

include("inference/Inference.jl")

include("quick.jl")

end # module


text:
export create,
    retrieve,
    update,
    delete

"""
    create(c::Client, ::Type{T}, args...; kwargs...) -> Future{Response}

Create, add, send, etc.

## Examples
Sending a [`Message`](@ref):
```julia
create(c, Message, channel; content="foo")
```
Creating a new [`DiscordChannel`](@ref):
```julia
create(c, DiscordChannel, guild; name="bar")
```
Banning a [`Member`](@ref):
```julia
create(c, Ban, guild, member; reason="baz")
```
"""
function create end

"""
    retrieve(c::Client, ::Type{T}, args...; kwargs...) -> Future{Response}

Retrieve, get, list, etc.

## Examples
Getting the [`Client`](@ref)'s [`User`](@ref):
```julia
retrieve(c, User)
```
Getting a [`Guild`](@ref)'s [`DiscordChannel`](@ref)s:
```julia
retrieve(c, DiscordChannel, guild)
```
Getting an [`Invite`](@ref) to a [`Guild`](@ref) by code:
```julia
retrieve(c, Invite, "abcdef")
```
"""
function retrieve end

"""
    update(c::Client, x::T, args...; kwargs...) -> Future{Response}

Update, edit, modify, etc.

## Examples
Editing a [`Message`](@ref):
```julia
update(c, message; content="foo2")
```
Modifying a [`Webhook`](@ref):
```julia
update(c, webhook; name="bar2")
```
Updating a [`Role`](@ref):
```julia
update(c, role, guild; permissions=8)
```
"""
function update end

"""
    delete(c::Client, x::T, args...) -> Future{Response}

Delete, remove, discard, etc.

## Examples
Kicking a [`Member`](@ref):
```julia
delete(c, member)
```
Unbanning a [`Member`](@ref):
```julia
delete(c, ban, guild)
```
Deleting all [`Reaction`](@ref)s from a [`Message`](@ref) (note: this is the only
update/delete method which takes a type parameter):
```julia
delete(c, Reaction, message)
```
"""
function delete end

include("audit_log.jl")
include("ban.jl")
include("channel.jl")
include("emoji.jl")
include("guild_embed.jl")
include("guild.jl")
include("integration.jl")
include("invite.jl")
include("member.jl")
include("message.jl")
include("overwrite.jl")
include("reaction.jl")
include("role.jl")
include("user.jl")
include("voice_region.jl")
include("webhook.jl")

#=
Functions not covered here:
- trigger_typing_indicator
- group_dm_add_recipient
- group_dm_remove_recipient
- modify_guild_channel_positions
- modify_current_user_nick
- modify_guild_role_positions
- get_guild_prune_count
- begin_guild_prune
- sync_guild_integration (combine with update?)
- get_guild_vanity_url
- get_guild_widget_image
- leave_guild (combine with delete?)
- execute_webhook
- execute_slack_compatible_webhook
- execute_github_compatible_webhook
=#


text:
#these are all the control data, they hold the flags that guide the calculation
function pdos(calculation::Calculation, args...)
    @error "pdos reading not implemented for package $(eltype(calculation))."
end

function Emin_from_projwfc(calculation::Calculation, args...)
    @error "Emin_from_projwfc is not implemented for package $(eltype(calculation))."
end

function readoutput(c::Calculation; kwargs...)
    @error "Output parsing for package $(eltype(c)) not implemented."
end

"""
    outputdata(calculation::Calculation, file; extra_parse_funcs=[], print=true, overwrite=true)

If `file` exists this will parse it and return a `Dict` with the parsed data.
If `overwrite=false` and `calculation.outputdata` is not empty, this will be returned instead of reparsing the
output file.

`extra_parse_funcs` should be `Vector{Pair{String,Function}}`, where the string will be used to match `occursin(str, line)`
for each line of the output file. If this returns `true`, the function will be called as `func(results_dict, line, file)`.
The purpose is to allow for additional parsing that is not implemented, or for temporary non-standard values that are printed
while working on the DFT code, e.g. for debugging.

!!! note
    This only works for files where not all information is pulled out for yet,
    e.g. projwfc.x outputs are fully parsed already.

Example (from src/qe/fileio.jl):
```julia

function qe_parse_nat(results, line, f)
    results[:nat] = parse(Int, split(line)[end])
end

outputdata(job["scf"],joinpath(job, job["scf"].outfile), extra_parse_funcs = ["number of atoms/cell" => qe_parse_nat])
```
"""
function outputdata(calculation::Calculation, file;
                    extra_parse_funcs::Vector{<:Pair{String}} = Pair{String}[],
                    print = true, overwrite = true)
    if ispath(file)
        t = FileIO.readoutput(calculation, file; parse_funcs = extra_parse_funcs)
        return t === nothing ?
                              FileIO.parse_file(file,
                                                extra_parse_funcs) : t
    end
    print &&
        (@warn "File $file does not exist.")
    return Dict{Symbol,Any}()
end

verify_exec(args...) = Calculations.verify_exec(args...)

known_execs(args...) = Calculations.known_execs(args...)
load_exec(args...) = Calculations.load_exec(args...)

function save(e::Exec)
    @assert Calculations.isrunnable(e) "Exec is not runnable..."
    Calculations.save(e)
    return e
end


text:
function is_deadend(cell::Cell)
    return length(cell.connections) == 1
end

function get_deadends(grid::Grid)
    des = []
    for row in 1:(grid.nrows)
        for col in 1:(grid.ncols)
            cell = grid[row, col]
            if is_deadend(cell)
                push!(des, cell)
            end
        end
    end
    return des
end

function remove_deadends!(grid::Grid, p::Float64)
    @assert (p <= 1.0) & (p >= 0.0)
    deadends = shuffle(get_deadends(grid))
    for cell in deadends
        if is_deadend(cell) & (rand() < p)
            nbs = filter(x -> !is_connected(cell, x), neighbors(grid, cell))
            if isempty(nbs)
                break
            end
            best = filter(x -> is_deadend(x), nbs)
            if isempty(best)
                best = nbs
            end
            neighbor = rand(best)
            connect!(cell, neighbor)
        end
    end
end


text:
module GLOM_RV_Example

const n_kep_parms = 6
include("general_functions.jl")
include("glo_rv_functions.jl")
export GLO_RV
include("RV_functions.jl")
include("fit_GLOM_functions.jl")
include("keplerian_derivatives.jl")
include("prior_functions.jl")
include("diagnostic_functions.jl")

end # module


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, DiffEqCallbacks,
      Test

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

prob = prob_ode_linear
ts = Vector{Float64}()
cb = FunctionCallingCallback((u,t,integrator)->push!(ts,t))
sol = solve(prob, Tsit5(), callback=cb)
@test !isempty(ts)

ts = Vector{Float64}()
cb = FunctionCallingCallback((u,t,integrator)->push!(ts,t),funcat= 0.0:0.25:1.0)
sol = solve(prob, Tsit5(), callback=cb)
@test collect(0.0:0.25:1.0) == ts

ts = Vector{Float64}()
cb = FunctionCallingCallback((u,t,integrator)->push!(ts,t),funcat= 0.0:0.25:1.0,func_everystep=true)
sol = solve(prob, Tsit5(), callback=cb)
@test intersect(collect(0.0:0.25:1.0),ts) == collect(0.0:0.25:1.0)


text:
function accept_reject_entry!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    i::Int,
    mcmc::T,
    model::MultigraphSNF,
    x_mode::Vector{Int} # Must also pass vectorised mode
    ) where {T<:Union{SnfMcmcSysGibbs,SnfMcmcRandGibbs}}

    # Proposal generation
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:mcmc.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    # Eval acceptance probability
    γ, d = (model.γ, model.d)
    log_α = γ * (d(x_curr, x_mode) - d(x_prop, x_mode))

    if log(rand()) < log_α
        @inbounds x_curr[i] = w_prop
        return true
    else
        @inbounds x_prop[i] = w_curr
        return false
    end
end

# Optimised version for the hamming distances between mutlisets
function accept_reject_entry_hamming!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    i::Int,
    mcmc::T,
    model::MultigraphSNF,
    x_mode::Vector{Int} # Must also pass vectorised mode
    ) where {T<:Union{SnfMcmcRandGibbs,SnfMcmcSysGibbs}}

    # Proposal generation
    @inbounds w_curr = x_curr[i]
    d = 1-2rand(0:1)
    w_tmp = w_curr + d * rand(1:mcmc.ν)
    w_prop = w_tmp ≥ 0 ? w_tmp : -w_tmp
    @inbounds x_prop[i] = w_prop

    # Eval acceptance probability
    γ, d = (model.γ, model.d)
    @inbounds w_mode = x_mode[i]
    log_α = γ * (abs(w_curr - w_mode) - abs(w_prop - w_mode))

    if log(rand()) < log_α
        @inbounds x_curr[i] = w_prop
        return true
    else
        @inbounds x_prop[i] = w_curr
        return false
    end
end

# Systematic Gibbs scan optimised for Hamming distance
function draw_sample_hamming!(
    out::Vector{Vector{Int}},
    mcmc::SnfMcmcSysGibbs,
    model::MultigraphSNF;
    burn_in::Int=mcmc.burn_in, lag::Int=mcmc.lag,
    init::Vector{Int}=adj_mat_to_vec(model.mode, directed=model.directed, self_loops=model.self_loops)
    )

    x_curr = copy(init)
    x_prop = copy(x_curr)
    x_mode = adj_mat_to_vec(
        model.mode,
        directed=model.directed,
        self_loops=model.self_loops
    )
    M = length(x_mode) # Number of possible unqiue edges
    sample_count = 1
    accept_count = 0
    iter_count = 0
    while sample_count ≤ length(out)
        # Store value if beyond burn-in and coherent with lags
        if (iter_count > burn_in) & (((iter_count-1) % lag)==0)
            @inbounds copy!(out[sample_count], x_curr)
            sample_count += 1
        end
        # Gibbs scan
        for i in 1:M
            was_accepted = accept_reject_entry_hamming!(
                x_curr, x_prop,
                i,
                mcmc, model,
                x_mode
            )
            accept_count += was_accepted
            iter_count += 1
        end

    end
    return accept_count / iter_count
end

function draw_sample!(
    out::Vector{Vector{Int}},
    mcmc::SnfMcmcSysGibbs,
    model::MultigraphSNF;
    burn_in::Int=mcmc.burn_in, lag::Int=mcmc.lag,
    init::Vector{Int}=adj_mat_to_vec(model.mode, directed=model.directed, self_loops=model.self_loops)
    )
    # For Cityblock (what I see as hamming for multisets) we can optimise
    x_curr = copy(init)
    x_prop = copy(x_curr)
    x_mode = adj_mat_to_vec(
        model.mode,
        directed=model.directed,
        self_loops=model.self_loops
    )
    M = length(x_mode) # Number of possible unqiue edges
    sample_count = 1
    accept_count = 0
    iter_count = 0
    while sample_count ≤ length(out)
        # Store value if beyond burn-in and coherent with lags
        if (iter_count > burn_in) & (((iter_count-1) % lag)==0)
            @inbounds copy!(out[sample_count], x_curr)
            sample_count += 1
        end
        # Gibbs scan
        for i in 1:M
            was_accepted = accept_reject_entry!(
                x_curr, x_prop,
                i,
                mcmc, model,
                x_mode
            )
            accept_count += was_accepted
            iter_count += 1
        end

    end
    return accept_count / iter_count
end

# Random scann Gibbs
# -------------------

function draw_sample_hamming!(
    out::Vector{Vector{Int}},
    mcmc::SnfMcmcRandGibbs,
    model::MultigraphSNF;
    burn_in::Int=mcmc.burn_in, lag::Int=mcmc.lag,
    init::Vector{Int}=adj_mat_to_vec(model.mode, directed=model.directed, self_loops=model.self_loops)
    )

    x_curr = copy(init)
    x_prop = copy(x_curr)
    x_mode = adj_mat_to_vec(
        model.mode,
        directed=model.directed,
        self_loops=model.self_loops
    )
    M = length(x_mode) # Number of possible unqiue edges
    sample_count = 1
    accept_count = 0
    iter_count = 0

    while sample_count ≤ length(out)
        # Store value if beyond burn-in and coherent with lags
        if (iter_count > burn_in) & (((iter_count-1) % lag)==0)
            @inbounds copy!(out[sample_count], x_curr)
            sample_count += 1
        end
        # Gibbs updates random component
        i = rand(1:M)
        was_accepted = accept_reject_entry_hamming!(
            x_curr, x_prop,
            i,
            mcmc, model,
            x_mode
        )
        accept_count += was_accepted
        iter_count += 1

    end
    return accept_count / iter_count
end

function draw_sample!(
    out::Vector{Vector{Int}}, # Vectorised implementation
    mcmc::SnfMcmcRandGibbs,
    model::MultigraphSNF;
    burn_in::Int=mcmc.burn_in, lag::Int=mcmc.lag,
    init::Vector{Int}=adj_mat_to_vec(model.mode, directed=model.directed, self_loops=model.self_loops)
    )

    x_curr = copy(init)
    x_prop = copy(x_curr)
    x_mode = adj_mat_to_vec(
        model.mode,
        directed=model.directed,
        self_loops=model.self_loops
    )
    M = length(x_mode) # Number of possible unqiue edges
    sample_count = 1
    accept_count = 0
    iter_count = 0

    while sample_count ≤ length(out)
        # Store value if beyond burn-in and coherent with lags
        if (iter_count > burn_in) & (((iter_count-1) % lag)==0)
            @inbounds copy!(out[sample_count], x_curr)
            sample_count += 1
        end
        # Gibbs updates random component
        i = rand(1:M)
        was_accepted = accept_reject_entry!(
            x_curr, x_prop,
            i,
            mcmc, model,
            x_mode
        )
        accept_count += was_accepted
        iter_count += 1

    end
    return accept_count / iter_count
end


text:
macro FZ_RR()#(N, F_yf, F_yr, V, Ax, r)
	code = quote
		# define the lateral tire forces
		F_yf = @F_YF()
		F_yr = @F_YR()
		FZ_rr = zeros(length(Ax),1)
		for ii in eachindex(Ax)
			FZ_rr[ii] = 0.5*(FzR0 + KZX*(Ax[ii] - V[ii]*R[ii])) + KZYR*((F_yf[ii] + F_yr[ii])/m)
		end
		FZ_rr
	end
	return esc(code)
end


text:
function colorbar!(fg::FigureGrid; position=:right,
                   vertical=default_isvertical(position), kwargs...)

    guide_pos = guides_position(fg.figure, position)
    return colorbar!(guide_pos, fg; vertical, kwargs...)
end

"""
    colorbar!(figpos, grid; kwargs...)

Compute colorbar for `grid` (which should be the output of [`draw!`](@ref)) and draw it in
position `figpos`. Attributes allowed in `kwargs` are the same as `MakieLayout.Colorbar`.
"""
function colorbar!(figpos, grid; kwargs...)
    colorbar = compute_colorbar(grid)
    return isnothing(colorbar) ? nothing : Colorbar(figpos; colorbar..., kwargs...)
end

compute_colorbar(fg::FigureGrid) = compute_colorbar(fg.grid)
function compute_colorbar(grid::Matrix{AxisEntries})
    colorscales = filter(!isnothing, [get(ae.continuousscales, :color, nothing) for ae in grid])
    isempty(colorscales) && return
    colorscale = reduce(mergescales, colorscales)
    label = getlabel(colorscale)
    limits = colorscale.extrema
    colormap = current_default_theme().colormap[]
    # FIXME: handle separate colorbars
    for entry in entries(grid)
        colormap = to_value(get(entry, :colormap, colormap))
    end
    return (; label, limits, colormap)
end