text:
#==============================================================================#
# AWSIAM.jl
#
# IAM API. See http://aws.amazon.com/documentation/iam/
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

__precompile__()

module AWSIAM

export iam,
       iam_whoami,
       iam_create_user, iam_create_access_key,
       user_arn, role_arn

using AWSCore
using Retry
using SymDict

role_arn(aws::AWSConfig, role_name) = arn(aws, "iam", "role/$role_name")
user_arn(aws::AWSConfig, user_name) = arn(aws, "iam", "user/$user_name")

function iam(aws::AWSConfig; Action="", args...)

    @repeat 4 try

        AWSCore.Services.iam(aws, Action, stringdict(args))

    catch e
        @retry if ecode(e) == "NoSuchEntity" end
    end
end

function sts(aws::AWSConfig; Action="", args...)

    AWSCore.Services.sts(aws, Action, stringdict(args))
end

function iam_whoami(aws::AWSConfig)

    sts(aws, Action = "GetCallerIdentity")["Arn"]
end

function iam_create_user(aws::AWSConfig, user_name)

    iam(aws, Action = "CreateUser", UserName = user_name)

    iam_create_access_key(aws, user_name)
end

#=

function iam_delete_access_key(aws::AWSConfig, user_name)

    r = iam(aws, "ListAccessKeys", Dict("UserName" => user_name))
    # ListAccessKeysResult AccessKeyMetadata] {
 #       set key [get $key AccessKeyId]
    iam(aws, "DeleteAccessKey", Dict("UserName" => user_name, "AccessKeyId" => key))
end

=#

function iam_create_access_key(aws::AWSConfig, user_name)

    r = iam(aws, Action = "CreateAccessKey", UserName = user_name)
    r = r["AccessKey"]
    AWSCredentials(r["AccessKeyId"], r["SecretAccessKey"])
end

#=
function iam_delete_user(aws, user_name)

    r = iam(aws, "ListUserPolicies", Dict("UserName" => user_name))
#    for {- policy} in
#                    ListUserPoliciesResult PolicyNames] {
        iam(aws, "DeleteUserPolicy", Dict("UserName" => user_name,
                                          "PolicyName" => policy))
#    }

    iam_delete_access_key(aws::AWSConfig, user_name)

    r = iam(aws, "ListMFADevices", Dict("UserName" => user_name))
#    for {- key} in
#                    ListMFADevicesResult MFADevices] {
        sn = get(r, "SerialNumber")
        iam(aws, "DeactivateMFADevice", Dict("UserName" => user_name,
                                             "SerialNumber" => sn))
        iam(aws, "DeleteVirtualMFADevice", Dict("SerialNumber" => sn))
#    end

    iam(aws, "DeleteUser", Dict("UserName" => user_name))
end

function iam_put_user_policy(aws::AWSConfig, user_name, policy_name, policy)

    iam(aws, "PutUserPolicy",
             Dict("UserName" => user_name,
                  "PolicyName" => policy_name,
                  "PolicyDocument" => iam_format_policy(policy)))
end

function iam_format_policy(policy_statement)

#    json [dict create Version 2012-10-17 Statement $policy_statement]
end

=#

function iam_create_role(aws::AWSConfig, name; path="/")

    policy = """{
      "Version": "2012-10-17",
      "Statement": [
        {
          "Effect": "Allow",
          "Principal": {
            "Service": "ec2.amazonaws.com"
          },
          "Action": "sts:AssumeRole"
        }
      ]
    }"""

    r = iam(aws, Action = "CreateRole",
                 AssumeRolePolicyDocument = policy,
                 Path = path,
                 RoleName = name)
end

#=

proc assume_aws_sts_role {aws duration_s name role {policy {}} {mfa {}}} {

    dset args Version 2011-06-15
    dset args DurationSeconds $duration_s
    dset args RoleArn [aws_iam_role_arn $aws $role]
    dset args RoleSessionName $name

    if {$policy != {}} {
        dset args Policy [aws_iam_policy_format $policy]
    }
    if {$mfa != {}} {
        dset args SerialNumber [lindex $mfa 0]
        dset args TokenCode [lindex $mfa 1]
    }

    set response [aws_request $aws sts Action AssumeRole {*}$args]

    set TokenArn [get $response AssumeRoleResult AssumedRoleUser Arn]
    set creds [get $response AssumeRoleResult Credentials]
    dict with creds {}
    subst {
        AWSAccessKeyId $AccessKeyId
        AWSSecretKey   $SecretAccessKey
        AWSToken       $SessionToken
        AWSUserArn     $TokenArn
        Expiration     $Expiration
    }
}

proc create_aws_iam_instance_profile {aws name {path /}} {
    # Create an Instance Profile for use with and EC2 instance.

    try {

        aws_iam $aws DeleteInstanceProfile InstanceProfileName $name

    } trap NoSuchEntity {} {}

    set response [aws_iam $aws CreateInstanceProfile \
                               InstanceProfileName $name \
                               Path $path]

    get $response CreateInstanceProfileResult InstanceProfile Arn
}

proc add_role_to_aws_iam_instance_profile {aws ip_name role_name} {
    # Add "role_name to "ip_name".

    set response [aws_iam $aws AddRoleToInstanceProfile \
                               InstanceProfileName  $ip_name \
                               RoleName $role_name]
}

proc create_aws_iam_role {aws name {path /} {options {}}} {
    # Create a Role.

    puts "Creating Role \"$name\"..."

    # Allow EC2 to assume this role...
    # Allow this account number to assume this role...
    if {"-require-mfa" in $options} {
        set assume_role_policy [aws_iam_policy_format [tcl_subst {
            Effect Allow
            Action sts:AssumeRole
            Principal {AWS "arn:aws:iam::[aws_account_number $aws]:root"}
            Condition {Null {JSONDict: aws:MultiFactorAuthAge false}}
        }]]
    } else {
        set assume_role_policy [aws_iam_policy_format [tcl_subst {
            Effect Allow
            Action sts:AssumeRole
            Principal {
                Service "ec2.amazonaws.com"
                AWS "arn:aws:iam::[aws_account_number $aws]:root"
            }
        }]]
    }

    # Clean up old role policies...
    try {

        set response [aws_iam $aws ListRolePolicies RoleName $name]
        set policy_names [get $response ListRolePoliciesResult PolicyNames]
        foreach {member policy_name} $policy_names {
            aws_iam $aws DeleteRolePolicy RoleName $name PolicyName $policy_name
        }
    } trap NoSuchEntity {} {}

    # Remove role from instance profiles...
    try {

        set response [aws_iam $aws ListInstanceProfilesForRole RoleName $name]
        set ip_names [get $response ListInstanceProfilesForRoleResult \
                                         InstanceProfiles]
        foreach {member ip} $ip_names {
            aws_iam $aws RemoveRoleFromInstanceProfile \
                         InstanceProfileName [get $ip InstanceProfileName]\
                         RoleName $name
        }
    } trap NoSuchEntity {} {}

    # Delete role...
    try {

        aws_iam $aws DeleteRole RoleName $name

    } trap NoSuchEntity {} {}

    set response [aws_iam $aws CreateRole \
                               AssumeRolePolicyDocument $assume_role_policy \
                               Path $path \
                               RoleName $name]

    get $response CreateRoleResult Role Arn
}

proc put_aws_iam_role_policy {aws role_name policy_name policy} {
    # Attach "policy" as "policy_name" to "role_name".

    set response [aws_iam $aws PutRolePolicy \
                               PolicyDocument [aws_iam_policy_format $policy] \
                               PolicyName $policy_name \
                               RoleName $role_name]
}

proc create_aws_iam_mfa {aws name {path /}} {

    aws_iam $aws DeleteVirtualMFADevice \
                 SerialNumber [aws_arn $aws iam mfa$path$name]

    : [aws_iam $aws CreateVirtualMFADevice \
                    VirtualMFADeviceName $name \
                    Path $path] \
    | get CreateVirtualMFADeviceResult VirtualMFADevice
}

proc enable_aws_iam_mfa {aws mfa_name user_name code1 code2} {

    aws_iam $aws EnableMFADevice \
                 UserName $user_name \
                 SerialNumber [aws_arn $aws iam mfa/$mfa_name] \
                 AuthenticationCode1 $code1 \
                 AuthenticationCode2 $code2
}

=#

end # module AWSIAM

#==============================================================================#
# End of file.
#==============================================================================#


text:
module VSCodeServer

using REPL, Sockets, Base64, Pkg, UUIDs
import Base: display, redisplay
import Dates

function __init__()
    atreplinit() do repl
        @async try
            hook_repl(repl)
        catch err
            Base.display_error(err, catch_backtrace())
        end
    end

    push!(Base.package_callbacks, pkgload)
end

include("../../JSON/src/JSON.jl")

module JSONRPC
    import ..JSON
    import ..UUIDs

    include("../../JSONRPC/src/core.jl")
end

include("misc.jl")
include("trees.jl")
include("repl.jl")
include("../../../debugger/debugger.jl")
include("gridviewer.jl")

const INLINE_RESULT_LENGTH = 100
const MAX_RESULT_LENGTH = 10_000

function get_variables()
    M = Main
    variables = []
    clear_lazy()

    for n in names(M, all=true, imported=true)
        !isdefined(M, n) && continue
        Base.isdeprecated(M, n) && continue

        x = getfield(M, n)
        x === vscodedisplay && continue
        x === VSCodeServer && continue
        x === Main && continue

        n_as_string = string(n)
        startswith(n_as_string, "#") && continue
        t = typeof(x)

        rendered = treerender(x)

        push!(variables, Dict(
            "type" => string(t),
            "value" => get(rendered, :head, "???"),
            "name" => n_as_string,
            "id" => get(rendered, :id, get(get(rendered, :child, Dict()), :id, false)),
            "haschildren" => get(rendered, :haschildren, false),
            "lazy" => get(rendered, :lazy, false),
            "icon" => get(rendered, :icon, ""),
            "canshow" => can_display(x)
        ))
    end

    return variables
end

struct InlineDisplay <: AbstractDisplay end

function ends_with_semicolon(x)
    return REPL.ends_with_semicolon(split(x,'\n',keepempty = false)[end])
end

function sendDisplayMsg(kind, data)
    JSONRPC.send_notification(conn_endpoint[], "display", Dict{String,String}("kind"=>kind, "data"=>data))
end

"""
    render(x)

Produce a representation of `x` that can be displayed by a UI. Must return a dictionary with
the following fields:
- `inline`: Short one-line plain text representation of `x`. Typically limited to `INLINE_RESULT_LENGTH` characters.
- `all`: Plain text string (that may contain linebreaks and other signficant whitespace) to further describe `x`.
- `iserr`: Boolean. The frontend may style the UI differently depending on this value.
"""
function render(x)
    str = sprintlimited(MIME"text/plain"(), x, limit = MAX_RESULT_LENGTH)

    return Dict(
        "inline" => strlimit(first(split(str, "\n")), limit = INLINE_RESULT_LENGTH),
        "all" => str,
        "iserr" => false
    )
end

function render(::Nothing)
    return Dict(
        "inline" => "✓",
        "all" => "nothing",
        "iserr" => false
    )
end

struct EvalError
    err
    bt
end

function render(err::EvalError)
    str = sprintlimited(err.err, err.bt, func = Base.display_error, limit = MAX_RESULT_LENGTH)

    return Dict(
        "inline" => strlimit(first(split(str, "\n")), limit = INLINE_RESULT_LENGTH),
        "all" => str,
        "iserr" => true
    )
end
"""
    safe_render(x)

Calls `render`, but catches errors in the display system.
"""
function safe_render(x)
    try
        render(x)
    catch err
        out = render(EvalError(err, catch_backtrace()))
        out["inline"] = string("Display Error: ", out["inline"])
        out["all"] = string("Display Error: ", out["all"])
    end
end

function module_from_string(mod)
    ms = split(mod, '.')

    out = Main

    loaded_module = findfirst(==(first(ms)), string.(Base.loaded_modules_array()))

    if loaded_module !== nothing
        out = Base.loaded_modules_array()[loaded_module]
        popfirst!(ms)
    end

    for m in Symbol.(ms)
        if isdefined(out, m)
            resolved = getfield(out, m)

            if resolved isa Module
                out = resolved
            else
                return out
            end
        end
    end

    return out
end

is_module_loaded(mod) = mod == "Main" || module_from_string(mod) !== Main

function get_modules(toplevel = nothing, mods = Set(Module[]))
    top_mods = toplevel === nothing ? Base.loaded_modules_array() : [toplevel]

    for mod in top_mods
        push!(mods, mod)

        for name in names(mod, all=true)
            if !Base.isdeprecated(mod, name) && isdefined(mod, name)
                thismod = getfield(mod, name)
                if thismod isa Module && thismod !== mod && !(thismod in mods)
                    push!(mods, thismod)
                    get_modules(thismod, mods)
                end
            end
        end
    end
    mods
end

const conn_endpoint = Ref{Union{Nothing,JSONRPC.JSONRPCEndpoint}}(nothing)

isactive() = conn_endpoint[] !== nothing

function serve(args...; is_dev = false, crashreporting_pipename::Union{AbstractString,Nothing}=nothing)
    conn = connect(args...)
    conn_endpoint[] = JSONRPC.JSONRPCEndpoint(conn, conn)
    run(conn_endpoint[])

    if is_dev
        @async while true
            try
                Base.invokelatest(handle_message)
            catch err
                Base.display_error(err, catch_backtrace())
            end
        end
    else
        @async try
            while true
                Base.invokelatest(handle_message; crashreporting_pipename = crashreporting_pipename)
            end
        catch err
            Base.display_error(err, catch_backtrace())
        end
    end
end

function handle_message(; crashreporting_pipename::Union{AbstractString,Nothing}=nothing)
    msg = JSONRPC.get_next_message(conn_endpoint[])
    if msg["method"] == "repl/runcode"
        params = msg["params"]

        source_filename = params["filename"]
        code_line = params["line"]
        code_column = params["column"]
        source_code = params["code"]
        mod = params["module"]

        resolved_mod = try
            module_from_string(mod)
        catch err
            # maybe trigger error reporting here
            Main
        end

        show_code = params["showCodeInREPL"]
        show_result = params["showResultInREPL"]

        JSONRPC.send_notification(conn_endpoint[], "repl/starteval", nothing)
        hideprompt() do
            if isdefined(Main, :Revise) && isdefined(Main.Revise, :revise) && Main.Revise.revise isa Function
                let mode = get(ENV, "JULIA_REVISE", "auto")
                    mode == "auto" && Main.Revise.revise()
                end
            end
            if show_code
                for (i,line) in enumerate(eachline(IOBuffer(source_code)))
                    if i==1
                        printstyled("julia> ", color=:green)
                        print(' '^code_column)
                    else
                        # Indent by 7 so that it aligns with the julia> prompt
                        print(' '^7)
                    end

                    println(line)
                end
            end

            withpath(source_filename) do
                res = try
                    ans = Base.invokelatest(include_string, resolved_mod, '\n'^code_line * ' '^code_column *  source_code, source_filename)
                    @eval Main ans = $(QuoteNode(ans))
                catch err
                    EvalError(err, catch_backtrace())
                end

                if show_result
                    if res isa EvalError
                        Base.display_error(stderr, res.err, res.bt)
                    elseif res !== nothing && !ends_with_semicolon(source_code)
                        Base.invokelatest(display, res)
                    end
                else
                    try
                        Base.invokelatest(display, InlineDisplay(), res)
                    catch err
                        if !(err isa MethodError)
                            printstyled(stderr, "Display Error: ", color = Base.error_color(), bold = true)
                            Base.display_error(stderr, err, catch_backtrace())
                        end
                    end
                end

                JSONRPC.send_success_response(conn_endpoint[], msg, safe_render(res))
            end
        end
        JSONRPC.send_notification(conn_endpoint[], "repl/finisheval", nothing)
    elseif msg["method"] == "repl/getvariables"
        vars = get_variables()

        JSONRPC.send_success_response(conn_endpoint[], msg, vars)
    elseif msg["method"] == "repl/getlazy"
        JSONRPC.send_success_response(conn_endpoint[], msg, get_lazy(msg["params"]))
    elseif msg["method"] == "repl/showingrid"
        try
            var = Core.eval(Main, Meta.parse(msg["params"]))

            Base.invokelatest(internal_vscodedisplay, var)
        catch err
            Base.display_error(err, catch_backtrace())
        end
    elseif msg["method"] == "repl/loadedModules"
        JSONRPC.send_success_response(conn_endpoint[], msg, string.(collect(get_modules())))
    elseif msg["method"] == "repl/isModuleLoaded"
        mod = msg["params"]

        is_loaded = is_module_loaded(mod)

        JSONRPC.send_success_response(conn_endpoint[], msg, is_loaded)
    elseif msg["method"] == "repl/startdebugger"
        hideprompt() do
            debug_pipename = msg["params"]
            try
                VSCodeDebugger.startdebug(debug_pipename)
            catch err
                VSCodeDebugger.global_err_handler(err, catch_backtrace(), crashreporting_pipename, "Debugger")
            end
        end
    end
end

function display(d::InlineDisplay, ::MIME{Symbol("image/png")}, x)
    payload = stringmime(MIME("image/png"), x)
    sendDisplayMsg("image/png", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/png")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)
    sendDisplayMsg("image/svg+xml", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("text/html")}, x)
    payload = stringmime(MIME("text/html"), x)
    sendDisplayMsg("text/html", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("text/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}, x)
    payload = stringmime(MIME("juliavscode/html"), x)
    sendDisplayMsg("juliavscode/html", payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("juliavscode/html")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v2+json"), x)
    sendDisplayMsg("application/vnd.vegalite.v2+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v3+json"), x)
    sendDisplayMsg("application/vnd.vegalite.v3+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v3+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v4+json"), x)
    sendDisplayMsg("application/vnd.vegalite.v4+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v3+json"), x)
    sendDisplayMsg("application/vnd.vega.v3+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v3+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v4+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v4+json"), x)
    sendDisplayMsg("application/vnd.vega.v4+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v4+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v5+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v5+json"), x)
    sendDisplayMsg("application/vnd.vega.v5+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v5+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.plotly.v1+json")}, x)
    payload = stringmime(MIME("application/vnd.plotly.v1+json"), x)
    sendDisplayMsg("application/vnd.plotly.v1+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.dataresource+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.dataresource+json")}, x)
    payload = stringmime(MIME("application/vnd.dataresource+json"), x)
    sendDisplayMsg("application/vnd.dataresource+json", payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("application/vnd.dataresource+json")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.plotly.v1+json")}) = true

const DISPLAYABLE_MIMES = [
    "application/vnd.vegalite.v4+json",
    "application/vnd.vegalite.v3+json",
    "application/vnd.vegalite.v2+json",
    "application/vnd.vega.v5+json",
    "application/vnd.vega.v4+json",
    "application/vnd.vega.v3+json",
    "application/vnd.plotly.v1+json",
    "juliavscode/html",
    # "text/html",
    "image/svg+xml",
    "image/png"
]

function can_display(x)
    for mime in DISPLAYABLE_MIMES
        if showable(mime, x)
            return true
        end
    end

    if showable("application/vnd.dataresource+json", x)
        return true
    end

    istable = Base.invokelatest(_isiterabletable, x)

    if istable === missing || istable === true || x isa AbstractVector || x isa AbstractMatrix
        return true
    end

    return false
end

function Base.display(d::InlineDisplay, x)
    for mime in DISPLAYABLE_MIMES
        if showable(mime, x)
            return display(d, mime, x)
        end
    end

    throw(MethodError(display,(d,x)))
end

function _display(d::InlineDisplay, x)
    if showable("application/vnd.dataresource+json", x)
        display(d, "application/vnd.dataresource+json", x)
    else
        try
            display(d, x)
        catch err
            if err isa MethodError
                @warn "Cannot display values of type $(typeof(x)) in VS Code."
            else
                rethrow(err)
            end
        end
    end
end

const tabletraits_uuid = UUIDs.UUID("3783bdb8-4a98-5b6b-af9a-565f29a5fe9c")
const datavalues_uuid = UUIDs.UUID("e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5")

global _isiterabletable = i -> false
global _getiterator = i -> i

function pkgload(pkg)
    if pkg.uuid==tabletraits_uuid
        x = Base.require(pkg)

        global _isiterabletable = x.isiterabletable
        global _getiterator = x.getiterator
    elseif pkg.uuid==datavalues_uuid
        x = Base.require(pkg)

        eval(
            quote
                function JSON_print_escaped(io, val::$(x.DataValue))
                    $(x.isna)(val) ? print(io, "null") : JSON_print_escaped(io, val[])
                end

                julia_type_to_schema_type(::Type{T}) where {S, T<:$(x.DataValue){S}} = julia_type_to_schema_type(S)
            end
        )
    end
end

function remove_lln!(ex::Expr)
    for i in length(ex.args):-1:1
        if ex.args[i] isa LineNumberNode
            deleteat!(ex.args, i)
        elseif ex.args[i] isa Expr
            remove_lln!(ex.args[i])
        end
    end
end

function internal_vscodedisplay(x)
    if showable("application/vnd.dataresource+json", x)
        _display(InlineDisplay(), x)
    elseif _isiterabletable(x)===true
        buffer = IOBuffer()
        io = IOContext(buffer, :compact=>true)
        printdataresource(io, _getiterator(x))
        buffer_asstring = CachedDataResourceString(String(take!(buffer)))
        _display(InlineDisplay(), buffer_asstring)
    elseif _isiterabletable(x)===missing
        try
            buffer = IOBuffer()
            io = IOContext(buffer, :compact=>true)
            printdataresource(io, _getiterator(x))
            buffer_asstring = CachedDataResourceString(String(take!(buffer)))
            _display(InlineDisplay(), buffer_asstring)
        catch err
            _display(InlineDisplay(), x)
        end
    elseif x isa AbstractVector || x isa AbstractMatrix
        buffer = IOBuffer()
        io = IOContext(buffer, :compact=>true)
        print_array_as_dataresource(io, _getiterator(x))
        buffer_asstring = CachedDataResourceString(String(take!(buffer)))
        _display(InlineDisplay(), buffer_asstring)
    else
        _display(InlineDisplay(), x)
    end
end

vscodedisplay(x) = internal_vscodedisplay(x)
vscodedisplay() = i -> vscodedisplay(i)

macro enter(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/enter", $(string(command))))
end

macro run(command)
    remove_lln!(command)
    :(JSONRPC.send_notification(conn_endpoint[], "debugger/run", $(string(command))))
end

export vscodedisplay, @enter, @run

end  # module

ext: jl


text:
module Comm

using Funs, MPI

## import Base.LPROC, Base.PGRP

export ProcID, procid
# export myproc, procs, nprocs
export run_main
export rexec   # @rexec
export rexec_everywhere   # @rexec_everywhere

typealias ProcID Int

const OPTIMIZE_SELF_COMMUNICATION = false   #TODO true
const USE_MPI_MANAGER = false

type CommInfo
    was_initialized::Bool
    comm::MPI.Comm
    rank::Int
    size::Int
    manager::MPIManager
    CommInfo() = new()
end
const comminfo = CommInfo()

function init()
    comminfo.was_initialized = MPI.Initialized()
    if !comminfo.was_initialized
        MPI.Init()
    end
    # comminfo.comm = MPI.Comm_dup(MPI.COMM_WORLD)
    comminfo.comm = MPI.COMM_WORLD
    comminfo.rank = MPI.Comm_rank(comminfo.comm)
    comminfo.size = MPI.Comm_size(comminfo.comm)
    if USE_MPI_MANAGER
        # This call returns only on the root process
        @assert comminfo.comm == MPI.COMM_WORLD
        comminfo.manager = MPI.start(MPI_TRANSPORT_ALL)
    else
        ## # Override some Base settings to describe our multi-process setup
        ## LPROC.id = myproc()
        ## PGRP.workers = collect(procs())
    end
end

function finalize()
    if USE_MPI_MANAGER
        if !comminfo.was_initialized && !MPI.Finalized()
            @mpi_do comminfo.manager MPI.Finalize()
        end
    else
        ## # Undo override above to prevent errors during shutdown
        ## LPROC.id = 1
        ## PGRP.workers = []
        if !comminfo.was_initialized && !MPI.Finalized()
            MPI.Finalize()
        end
    end
end

myproc() = comminfo.rank+1
nprocs() = comminfo.size
procs() = 1:nprocs()

type CommState
    use_recv_loop::Bool
    stop_sending::Bool
    stop_receiving::Bool
    CommState() = new()
end
const commstate = CommState()

const TAG = 0
const META_TAG = 1

# Prevent unnecessary specialization
immutable Item
    item::Any
end
(item::Item)() = item.item()

# TODO: Use TestSome
function send_item(p::Int, t::Int, item::Item)
    @assert commstate.use_recv_loop
    req = MPI.isend(item, p-1, t, comminfo.comm)
    while t==META_TAG || !commstate.stop_sending
        flag, status = MPI.Test!(req)
        flag && return
        yield()
    end
    # MPI.cancel(req)
    nothing
end

function recv_item(p::Int, t::Int)
    @assert commstate.use_recv_loop
    while t==META_TAG || !commstate.stop_receiving
        flag, item, status =
            MPI.irecv(p==0 ? MPI.ANY_SOURCE : p-1, t, comminfo.comm)
        flag && return item::Item
        yield()
    end
    Item(nothing)
end

function terminate()
    # TODO: Use MPI.Ibarrier
    @assert commstate.use_recv_loop
    # Determine parents and children for a binary tree
    pmax = div(myproc(), 2)
    pmin = max(pmax, 1)
    cmin = 2*myproc()
    cmax = min(2*myproc()+1, nprocs())
    # Stage 1: Stop sending
    # Wait for termination message from parent
    @sync for p in pmin:pmax
        @async recv_item(p, META_TAG)
    end
    # Stop sending
    commstate.stop_sending = true
    # Send termination message to children
    @sync for c in cmin:cmax
        @async send_item(c, META_TAG, Item(nothing))
    end
    # Wait for termination confirmation from children
    @sync for c in cmin:cmax
        @async recv_item(c, META_TAG)
    end
    # Send termination confirmation to parent
    @sync for p in pmin:pmax
        @async send_item(p, META_TAG, Item(nothing))
    end
    # Stage 2: Stop receiving
    # Wait for second termination message from parent
    @sync for p in pmin:pmax
        @async recv_item(p, META_TAG)
    end
    commstate.stop_receiving = true
    # Send second termination message to children
    @sync for c in cmin:cmax
        @async send_item(c, META_TAG, Item(nothing))
    end
end

const main_running = Ref{Bool}(false)
function run_main(main; run_main_everywhere::Bool=false)
    @assert !main_running[]
    main_running[] = true
    init()
    @assert !USE_MPI_MANAGER
    commstate.use_recv_loop = !(OPTIMIZE_SELF_COMMUNICATION && nprocs()==1)
    r = nothing
    @sync begin
        if commstate.use_recv_loop
            commstate.stop_sending = false
            commstate.stop_receiving = false
            @async recv_loop()
        end
        if run_main_everywhere || myproc()==1
            r = main()
        end
        if commstate.use_recv_loop
            terminate()
        end
    end
    finalize()
    @assert main_running[]
    main_running[] = false
    r
end

function recv_loop()
    @assert commstate.use_recv_loop
    try
        while !commstate.stop_receiving
            run_task(recv_item(0, TAG))
        end
    catch ex
        info(ex)
    end
end

function run_task(item::Item)
    if !commstate.stop_sending
        @schedule item()
    end
    nothing
end

rexec(f, p::Int) = rexec(Item(f), p)
function rexec(item::Item, p::Int)
    if OPTIMIZE_SELF_COMMUNICATION && p == myproc()
        run_task(item)
    else
        @schedule send_item(p, TAG, item)
    end
    nothing
end

function _rexec_tree(item::Item)
    pmin = 2*myproc()
    pmax = min(2*myproc()+1, nprocs())
    for p in pmin:pmax
        rexec(p) do
            _rexec_tree(item)
        end
    end
    item()
    nothing
end
rexec_everywhere(f) = rexec_everywhere(Item(f))
function rexec_everywhere(item::Item)
    rexec(1) do
        _rexec_tree(item)
    end
    nothing
end

# macro rexec(p, expr)
#     :(rexec(()->$(esc(expr)), $(esc(p))))
# end
#
# macro rexec_everywhere(expr)
#     :(rexec_everywhere(()->$(esc(expr))))
# end

end


text:
abstract type AbstractObjective end

# Given callables to calculate objectives and partial first derivatives
# create a function that calculates both.
function make_fdf(x, F, f!, j!)
    function fj!(fx, jx, x)
        j!(jx, x)
        return f!(fx, x)
    end
end
function make_fdf(x, F::Number, f, g!)
    function fg!(gx, x)
        g!(gx, x)
        return f(x)
    end
end

# Initialize an n-by-n Jacobian
alloc_DF(x, F) = eltype(x)(NaN) .* vec(F) .* vec(x)'

# Initialize a gradient shaped like x
alloc_DF(x, F::T) where T<:Number = x_of_nans(x, promote_type(eltype(x), T))
# Initialize an n-by-n Hessian
function alloc_H(x, F::T) where T<:Number
  eltype(x)(NaN).*x*x'
end


text:
# Some built-in functions
println("println is a function")
println("parse(Int64, \"123\"): $(parse(Int64, "123"))")
println("parse(Float64, \"123.456\"): $(parse(Float64, "123.456"))")
println("trunc(Int64, 3.67123): $(trunc(Int64, 3.67123))")
println("trunc(Int64, -3.67123): $(trunc(Int64, -3.67123))")
println("Other examples: ", float(32), "\t", string(32), "\t", string(123.456))
println("Natural log and base, log(exp(1)): $(log(exp(1)))")
println("Other log bases, log2(10), log10(10): $(log2(10)), $(log10(10))")
println("Trigonometry, sin(pi), cos(pi), tan(pi): $(sin(pi)), $(cos(pi)), $(tan(pi))")
println("sqrt(2) = $(sqrt(2))")

# Defning functions
function printlines()
    println("This is the first line.")
    println("This is the second line.")
end

# printlines()

function printmorelines()
    printlines()
    printlines()
end

# printmorelines()

function printargument(arg)
    println(arg)
end

printargument(123)
printargument("123456")
printargument(cos(pi/2))

# Locality of variables and parameters
function sumthree(a, b, c)
    a = a + b
    total = a
    total += c
    return total
end
d, e, f = 1, 2, 3
println("sumthree(d, e, f): $(sumthree(d, e, f))")
println("d, e, f: $d, $e, $f")
# These lines will cause an error:
# println(a, b, c)
# println(total)

# Example of void function that returns nothing:
println(printlines())

# Example of a docstring
"""
polyline(t, n, len, angle)

Draws n line segments with the given length and
angle (in degrees) between them.  t is a turtle.
"""
function polyline(t, n, len, angle)
end

# Example of including errors
function fact(n)
    if !(n isa Int64)
        error("Factorial is only defined for integers.")
    elseif n < 0
        error("Factorial is not defined for negative integers.")
    elseif n == 0
        return 1
    else
        return n * fact(n-1)
    end
end

# println(fact(1.5))
# println(fact(-1))

# Shorthand notation for function definition
function func_ex(a=1, b=2)
    return a + 2b
end

# The above function definition is equivalent to the following:
func(a, b) = a + 2b
func(a) = func(a, 2)
func() = func(1, 2)
println("func(7, 5): $(func(7, 5))")
println("func(2): $(func(2))")
println("func(): $(func())")


text:
## ------------------------------------------------------------------
# USER CONFIG
_rbook_toml_name() = "ResearchBook.toml"
config_file(bookdir) = joinpath(bookdir, _rbook_toml_name())

function read_configfile(bookdir)
    cfile = config_file(bookdir)
    return TOML.parsefile(cfile)
end

## ------------------------------------------------------------------
function _find_bookdir(dir0::String)

    # search up
    root = homedir()
    bookdir = ""
    walkup(dir0) do path
        if isdir(path)
            (path == root) && return true
            bookfile = config_file(path)
            if isfile(bookfile)
                bookdir = path
                return true
            end
        end
    end
    return bookdir
end


text:
module TestFoldl

using Transducers: Map, ProductRF, TeeRF, __foldl__

using ..Utils: check_rrule

function test_simple()
    check_rrule(__foldl__, +, 0.0, 1:10)
    check_rrule(__foldl__, *, 1.0, 1:5)
    check_rrule(__foldl__, max, 0, [2, 6, 4, 9, 2, 7, 9, 5, 4, 9])
end

function test_tee()
    check_rrule(__foldl__, TeeRF(+, *), (0.0, 1.0), 1:5)
    check_rrule(__foldl__, TeeRF(+, *, max), (0.0, 1.0, 0.0), [2, 9, 6, 4])
    # TODO: Make sure that it can have non-differentiable parts in the accumulator?
    #=
    check_rrule(
        __foldl__,
        TeeRF(+, *, Map(string)'(*), max),
        (0.0, 1.0, "", 0),
        [2, 6, 4, 9, 2, 7, 9, 5, 4, 9],
    )
    =#
end

dup3(x) = (x, x, x)

function test_product()
    check_rrule(__foldl__, Map(dup3)'(ProductRF(+, *, *)), (0.0, 1.0, 0.0), [2, 9, 6, 4])
    # TODO: support zip-of-arrays
    #=
    check_rrule(__foldl__, ProductRF(+, *, *), (0.0, 1.0, 0.0), zip(1:5, 1:5, 'a':'e'))
    =#
end

end  # module


text:
using StatsBase

##########################################
# Uncertain theoretical distributions
##########################################
uncertain_theoretical_distributions = [
    UncertainValue(Uniform, -1, 1),
    UncertainValue(Normal, 0, 1),
    UncertainValue(Gamma, 1, 2),
    UncertainValue(Beta, 1, 2),
    UncertainValue(BetaPrime, 4, 5),
    UncertainValue(Frechet, 1, 1),
    UncertainValue(Binomial, 10, 0.3),
    UncertainValue(BetaBinomial, 100, 2, 3)
]

##########################################
# Kernel density estimates
##########################################
n = 10

uncertain_kde_estimates = [
    UncertainValue(rand(100))
]

##########################################
# Fitted theoretical distributions
##########################################
n = 10

uncertain_fitted_distributions = [
    UncertainValue(Uniform, rand(Uniform(-2, 2), n)),
    UncertainValue(Normal, rand(Normal(0, 1), n))
]

########################
# Uncertain populations
########################
pop1 = UncertainValue(
        [3.0, UncertainValue(Normal, 0, 1),
        UncertainValue(Gamma, 2, 3),
        UncertainValue(Uniform, rand(1000))],
        [0.5, 0.5, 0.5, 0.5]
    )

pop2 = UncertainValue([1, 2, 3], rand(3))
pop3 = UncertainValue([1.0, 2.0, 3.0], Weights(rand(3)))

# Uncertain population consisting of uncertain populations and other stuff
pop4 = UncertainValue([pop1, pop2], [0.1, 0.5])
pop5 = UncertainValue([pop1, pop2, 2, UncertainValue(Normal, -2, 3)], Weights(rand(4)));

uncertain_scalar_populations = [pop1, pop2, pop3, pop4, pop5]

##########################################
# Gather all examples
##########################################
uvals = [
    uncertain_scalar_populations;
    uncertain_theoretical_distributions;
    uncertain_kde_estimates;
    uncertain_fitted_distributions
];

##########################################
# Apply functions to datasets `n` times
##########################################
@test x1 = resample(median, uvals, 1000) isa Vector{T} where T <: Real


text:
using Distributed

## Define a new univariate Distribution type for MCPhylo.
## The definition must be placed within an unevaluated quote block.
@everywhere extensions = quote

  ## Load needed packages and import methods to be extended
  using Distributions
  import Distributions: minimum, maximum, logpdf

  ## Type declaration
  mutable struct NewUnivarDist <: ContinuousUnivariateDistribution
    mu::Float64
    sigma::Float64
  end

  ## The following method functions must be implemented

  ## Minimum and maximum support values
  minimum(d::NewUnivarDist) = -Inf
  maximum(d::NewUnivarDist) = Inf

  ## Normalized or unnormalized log-density value
  function logpdf(d::NewUnivarDist, x::Real)
    -log(d.sigma) - 0.5 * ((x - d.mu) / d.sigma)^2
  end

end

## Test the extensions in a temporary module (optional)
module Testing end
Core.eval(Testing, extensions)
d = Testing.NewUnivarDist(0.0, 1.0)
Testing.minimum(d)
Testing.maximum(d)
Testing.insupport(d, 2.0)
Testing.logpdf(d, 2.0)

## Add the extensions
@everywhere using MCPhylo
@everywhere eval(extensions)

## Implement a MCPhylo model using the new distribution
model = Model(

  y = Stochastic(1,
    (mu, s2) ->
      begin
        sigma = sqrt(s2)
        UnivariateDistribution[
          NewUnivarDist(mu[i], sigma) for i in 1:length(mu)
        ]
      end,
    false
  ),

  mu = Logical(1,
    (xmat, beta) ->
      xmat * beta,
    false
  ),

  beta = Stochastic(1,
    () -> MvNormal(2, sqrt(1000))
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Sampling Scheme
scheme = [NUTS(:beta),
          Slice(:s2, 3.0)]

## Sampling Scheme Assignment
setsamplers!(model, scheme)

## Data
line = Dict{Symbol, Any}(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]

## Initial Values
inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta => rand(Normal(0, 1), 2),
    :s2 => rand(Gamma(1, 1))
  )
  for i in 1:3
]

## MCMC Simulation
sim = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
describe(sim)

@test true


text:
"""
    IPR = Inverse Participation Ratio
"""
function get_IPRs(problem)
    n_modes = get_number_modes(problem)
    IPRs = zeros(n_modes)
    Threads.@threads for n = 1:n_modes
        Ψ² = get_ψ²(problem, n)
        Ψ⁴ = Ψ² .^ 2
        IPRs[n] = sum(Ψ⁴) / sum(Ψ²)
    end
    return IPRs
end
function get_PRs(problem)
    IPRs = get_IPRs(problem)
    return 1.0 ./ IPRs
end

function get_number_modes(problem::LinearOptics{Scalar})
    return problem.atoms.N
end
"""
    get_spectrum(problem; forceComputation=false)
return ωₙ, Γₙ
"""
function get_spectrum(problem; forceComputation = false)
    @debug "start: get spectrum"
    ωₙ, Γₙ = zeros(problem.atoms.N), zeros(problem.atoms.N)

    if is_spectrum_NOT_available(problem) || forceComputation
        H = get_interaction_matrix(problem)
        spectrum = eigen(H)

        problem.spectrum[:λ] = spectrum.values
        problem.spectrum[:ψ] = spectrum.vectors
        make_spectrum_available(problem)

        ωₙ, Γₙ = imag.(spectrum.values), -real.(spectrum.values)
        if any(Γₙ .< 0)
            @warn "some Γₙ were negatives and ignored without further investigation"
            Γₙ = abs.(Γₙ)
        end
    else
        ωₙ, Γₙ = imag.(problem.spectrum[:λ]), -real.(problem.spectrum[:λ])
    end

    @debug "end  : get spectrum"
    return ωₙ, Γₙ
end

function is_spectrum_NOT_available(problem)
    if problem.spectrum[:isSpectrumAvailable]
        return false
    else
        return true
    end
end
function make_spectrum_available(problem)
    problem.spectrum[:isSpectrumAvailable] = true
end

function get_ψ(problem::LinearOptics, n::Integer)
    ψ = view(problem.spectrum[:ψ], :, n)
    return ψ
end
function get_ψ²(problem::LinearOptics, n::Integer)
    ψ = problem.spectrum[:ψ][:, n]
    return abs2.(ψ)
end

"""
    get_localization_length(LinearOptics; forceComputation=false)

Main function to obtain the `Localization Length`, ξ, and its `Coefficient of Determination`, R1 (R¹).
"""
function get_localization_length(problem::LinearOptics; forceComputation = false)
    if is_localization_NOT_available(problem) && forceComputation == false
        return problem.data[:ξ], problem.data[:R1]
    end
    N = problem.atoms.N
    ξₙ = zeros(N)
    R¹ₙ = zeros(N)

    # create spectrum if neeeded
    get_spectrum(problem)

    for n = 1:N # faster without Threads.@threads
        DCM, ψ² = get_spatial_profile_single_mode(problem, n)
        ξₙ[n], R¹ₙ[n] = get_single_ξ_and_R1(DCM, ψ²)
    end
    problem.data[:ξ], problem.data[:R1] = ξₙ, R¹ₙ
    return ξₙ, R¹ₙ
end

function is_localization_NOT_available(problem)
    if haskey(problem.data, :ξ) && haskey(problem.data, :R1)
        return true
    else
        return false
    end

end

function get_spatial_profile_single_mode(problem, mode_index::Integer)
    r = problem.atoms.r
    try
        ψ²ₙ = get_ψ²(problem, mode_index)
        r_cm = get_coordinates_of_center_of_mass(r, ψ²ₙ)

        # Specific for the problem (ones needs to define this function)
        DCM = get_Distances_from_r_to_CM(r, r_cm)

        sort_spatial_profile!(DCM, ψ²ₙ)
        return DCM, ψ²ₙ
    catch
        @error(
            "Probably you tried to access a spectrum not created (run `get_spectrum(problem)`). Or `mode_index` is not valid"
        )
    end
end
"""
    get_coordinates_of_center_of_mass(r, Ψ²_mode)
"""
function get_coordinates_of_center_of_mass(r, Ψ²_mode)
    N = length(Ψ²_mode)
    dimensions = size(r, 1)
    r_cm = zeros(dimensions)
    j = 1
    for c in eachcol(r)
        r_cm .+= c .* Ψ²_mode[j]
        j += 1
    end
    return r_cm ./ sum(Ψ²_mode)
end

"""
    get_Distances_from_r_to_CM(r, r_CM)

YOU need to specify for how to get the distance from points.
By default, I compute the Canonical Euclidian Distance
"""
get_Distances_from_r_to_CM(r, r_CM) = get_Distance_A_to_b(r, r_CM)

function sort_spatial_profile!(DCM, ψ²ₙ)
    idx = sortperm(DCM)
    DCM[:] = DCM[idx]
    ψ²ₙ[:] = ψ²ₙ[idx]
    return nothing
end

### --------------- Classification r---------------
"""
    classify_modes(problem)
returns a tuple `(loc, sub, super)` with indices
"""
function classify_modes(problem; forceComputation = false)
    ωₙ, Γₙ = get_spectrum(problem; forceComputation = forceComputation)
    ξₙ, R¹ₙ = get_localization_length(problem; forceComputation = forceComputation)

    localized_modes = findall((Γₙ .< Γ) .* (R¹ₙ .≥ R1_threshold))
    sub_radiant_modes = findall((Γₙ .< Γ) .* (R¹ₙ .< R1_threshold))
    super_radiant_modes = findall(Γₙ .> Γ)
    return (loc = localized_modes, sub = sub_radiant_modes, super = super_radiant_modes)
end


text:
## functions for writing the midi codes as byte sequences

brev(v::UInt8) = v
brev(v::UInt16) = UInt8[(v>>8)%UInt8,(v)%UInt8]
brev(v::UInt32) = UInt8[(v>>24)%UInt8,(v>>16)%UInt8,(v>>8)%UInt8,(v)%UInt8]
brev(v::UInt64) = UInt8[(v>>56)%UInt8,(v>>48)%UInt8,(v>>40)%UInt8,(v>>32)%UInt8,(v>>24)%UInt8,(v>>16)%UInt8,(v>>8)%UInt8,(v)%UInt8]

msb14(v::UInt16) = (v>>7)%UInt8&0x7f
lsb14(v::UInt16) = (v)%UInt8&0x7f

function write14(out::IOStream, v::Integer)
	@assert 0 <= v < 2^14
	write(out,msb14(v%UInt16)) ;  write(out,lsb14(v%UInt16))
end

function varb(v::Integer)
	@assert 0 <= v < 2^28
	if v<2^7
		return UInt8[(v)%UInt8]
	elseif v<2^14
		return UInt8[(v>>7)%UInt8|0x80,(v)%UInt8&0x7f]
	elseif v<2^21
		return UInt8[(v>>14)%UInt8|0x80,(v>>7)%UInt8|0x80,(v)%UInt8&0x7f]
	else
		return UInt8[(v>>21)%UInt8|0x80,(v>>14)%UInt8|0x80,(v>>7)%UInt8|0x80,(v)%UInt8&0x7f]
	end
end

function writeFileHeaderT0(out::IOStream, tpq::Int)
	@assert 0 < tpq < 2^15
	write(out,ascii("MThd"))
 	write(out,UInt8[0x00,0x00,0x00,0x06,0x00,0x00,0x00,0x01])
	write(out,brev(tpq%UInt16))
end

function writeTrackHeader(out::IOStream, sz::Int)
	@assert 0 <= sz < 2^31
	write(out,ascii("MTrk"))
	write(out,brev(sz%UInt32))
end

function bytesTimeSig(num::Int, den::Int, clocks::Int, sub::Int)
	@assert 0 < num <= 32 ;  @assert 0 <= den <= 5
	@assert 0 < clocks < 128 ;  @assert 0 < sub <= 128
	UInt8[0xff,0x58,0x04,num%UInt8,den%UInt8,clocks%UInt8,sub%UInt8]
end

bytesDt(v::Int)  = varb(v)

function bytesTempo(bpm::Int)
	@assert 0 < bpm
	micspq = (div(10^6*60,bpm))%UInt32
	@assert 0 < micspq < 2^24
	vcat(UInt8[0xff,0x51,0x03],brev(micspq%UInt32)[2:4])
end

function bytesProgSelect(chan::Int, prog::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= prog < 128
	UInt8[0xC0|chan%UInt8,prog%UInt8]
end

function bytesNoteOn(chan::Int, pch::Int, vel::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= pch < 128 ;  @assert 0 <= vel < 128
	UInt8[0x90|chan%UInt8,pch%UInt8,vel%UInt8]
end

function bytesNoteOff(chan::Int, pch::Int, vel::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= pch < 128 ;  @assert 0 <= vel < 128
	UInt8[0x80|chan%UInt8,pch%UInt8,vel%UInt8]
end

function bytesPitchWheel(chan::Int, val::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= val < 2^14
	UInt8[0xe0|chan%UInt8,lsb14(val%UInt16),msb14(val%UInt16)]
end

function bytesChanAfter(chan::Int, val::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= val < 128
	UInt8[0xd0|chan%UInt8,val%UInt8]
end

function bytesPolyAfter(chan::Int, pitch::Int, val::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= pitch < 128 ;  @assert 0 <= val < 128
	UInt8[0xa0|chan%UInt8,pitch%UInt8,val%UInt8]
end

function bytesControl7(chan::Int, num::Int, val::Int)
	@assert 0 <= chan < 16 ;  @assert 64 <= num < 102 || 120 <= num < 128 ;  @assert 0 <= val < 128
	UInt8[0xb0|chan%UInt8,num%UInt8,val%UInt8]
end

function bytesControlMsb14(chan::Int, num::Int, val::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= num < 32 ;  @assert 0 <= val < 2^14
	UInt8[0xb0|chan%UInt8,num%UInt8,msb14(val%UInt16)]
end

function bytesControlLsb14(chan::Int, num::Int, val::Int)
	@assert 0 <= chan < 16 ;  @assert 0 <= num < 32 ;  @assert 0 <= val < 2^14
	UInt8[0xb0|chan%UInt8,(num+32)%UInt8,lsb14(val%UInt16)]
end

function bytesTrackEnd()
	UInt8[0xff,0x2f,0x00]
end

function bytesTrackMarker(mes::ASCIIString)
	len = length(mes) ;  @assert 0<len<128
	bs = UInt8[0xff,0x06,len%UInt8]
	for x in mes
		push!(bs,x%UInt8)
	end
	bs
end


text:
using TupleVectors
using Sobol
using Test
using TupleVectors: chainvec, unwrap
using MeasureTheory
using LinearAlgebra

@testset "TupleVectors.jl" begin
    ## Testing constructors
    w = randn(10)
    y = randn(10)
    nt = (w=w, y=y)
    A = TupleVector(nt)
    B = TupleVector(; w=w, y=y)
    @test A.w == B.w

    @test chainvec(3,5)[1] == 3

    x = randn(2,2)
    A = chainvec((x = x,),100);

    @test A.x[1] == x
    @test A[1].x == A.x[1]

    for i in 2:100
        A[i] = (x=randn(2,2),)
    end

    B = @with A begin
        (x = x' + x,)
    end

    C = @with B begin
        (eigs = eigvals!(x),)
    end

    using Sobol

    ω = SobolHypercube(2)

    πapprox = @with (;ω) 10000 begin
        x = rand(ω)
        y = rand(ω)
        val = x^2 + y^2 < 1 ? 4 : 0
        (;val)
    end

    @test mean(πapprox.val) ≈ π atol=0.1

end


text:
"""
TestSolution

"""
type TestSolution{hasinterp} <: DESolution
  t
  u
  interp
  dense
end
(T::TestSolution)(t) = T.interp(t)
TestSolution(t,u) = TestSolution{false}(t,u,nothing,false)
TestSolution(t,u,interp) = TestSolution{true}(t,u,interp,true)
TestSolution(interp::DESolution) = TestSolution{true}(nothing,nothing,interp,true)

"""
`appxtrue(sol::AbstractODESolution,sol2::TestSolution)`

Uses the interpolant from the higher order solution sol2 to approximate
errors for sol. If sol2 has no interpolant, only the final error is
calculated.
"""
function appxtrue{hasinterp}(sol::AbstractODESolution,sol2::TestSolution{hasinterp})
  if sol2.u == nothing && hasinterp
    sol2.u = sol2(sol.t)
    sol2.t = sol.t
  end
  errors = Dict(:final=>mean(abs.(sol[end]-sol2[end])))
  if sol2.dense
    timeseries_analytic = sol2(sol.t)
    errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol[:]-timeseries_analytic))
    errors[:l2] = sqrt(mean(vecvecapply((x)->float(x).^2,sol[:]-timeseries_analytic)))
    if !(typeof(sol) <: SDESolution) && sol.dense
      densetimes = collect(linspace(sol.t[1],sol.t[end],100))
      interp_u = sol(densetimes)
      interp_analytic = sol2(densetimes)
      interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),
                           :L2=>sqrt(mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
      errors = merge(errors,interp_errors)
    end
  end
  build_solution(sol,sol2.u,errors)
end

"""
`appxtrue!(sol::FEMSolution,sol2::FEMSolution)`

Adds the solution from `sol2` to the `FEMSolution` object `sol`.
Useful to add a quasi-true solution when none is known by
computing once at a very small time/space step and taking
that solution as the "true" solution
"""
function appxtrue!(sol::AbstractFEMSolution,sol2::AbstractFEMSolution)
  sol.u_analytic = sol2.u
  sol.errors = Dict(:l∞=>maximum(abs.(sol.u-sol.u_analytic)),:l2=>norm(sol.u-sol.u_analytic,2))
  nothing
end

"""
`appxtrue(sol::AbstractODESolution,sol2::AbstractODESolution)`

Uses the interpolant from the higher order solution sol2 to approximate
errors for sol. If sol2 has no interpolant, only the final error is
calculated.
"""
function appxtrue(sol::AbstractODESolution,sol2::AbstractODESolution)
  errors = Dict(:final=>mean(abs.(sol[end]-sol2[end])))
  if !(typeof(sol2) <: SDESolution) && sol2.dense
    timeseries_analytic = sol2(sol.t)
    errors = Dict(:final=>mean(abs.(sol[end]-sol2[end])),:l∞=>maximum(vecvecapply((x)->abs.(x),sol[:]-timeseries_analytic)),:l2=>sqrt(mean(vecvecapply((x)->float(x).^2,sol[:]-timeseries_analytic))))
    if !(typeof(sol) <: SDESolution) && sol.dense
      densetimes = collect(linspace(sol.t[1],sol.t[end],100))
      interp_u = sol(densetimes)
      interp_analytic = sol2(densetimes)
      interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),:L2=>sqrt(mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
      errors = merge(errors,interp_errors)
    end
  end
  build_solution(sol,sol2.u,errors)
end

max_stars_repo_name: davidmuaddib/DiffEqDevTools

max_issues_repo_name: davidmuaddib/DiffEqDevTools

max_forks_repo_name: davidmuaddib/DiffEqDevTools


text:
using Pkg

Pkg.build(verbose=true)

Pkg.add("GaussianProcess")
Pkg.add("Distributions")
Pkg.add("Arpack")
Pkg.add("CausalInference")
Pkg.add("Match")
Pkg.add("RecursiveArrayTools")

# for LightGraphs
Pkg.add("FileIO")
Pkg.add("LightGraphs")

Pkg.add("GraphPlot")
# a compatibility bug with Cairo
# https://github.com/GiovineItalia/Compose.jl/pull/360
# Pkg.add(PackageSpec(name="Compose", rev="master"))
Pkg.add("Compose")
Pkg.add("Cairo")
# Pkg.add("Fontconfig")
Pkg.add("MetaGraphs")

Pkg.add("Optim")
Pkg.add("NLopt")
Pkg.add("LineSearches")
Pkg.add("CSV")
Pkg.add("BSON")
Pkg.add("ForwardDiff")
Pkg.add("Tracker")

Pkg.add("Images")
Pkg.add("ProgressMeter")
Pkg.add("MLDatasets")
Pkg.add("ImageMagick")

Pkg.add("TensorBoardLogger")
Pkg.add("Plots")
Pkg.add("PyPlot")

Pkg.add("Interpolations")
Pkg.add("GaussianMixtures")
Pkg.add("StatsFuns")
Pkg.add("TensorOperations")
Pkg.add(Pkg.PackageSpec(name="TensorOperations", rev="master"))
Pkg.add(Pkg.PackageSpec(url="https://github.com/mcabbott/TensorGrad.jl", rev="master"))
Pkg.free("TensorOperations")
Pkg.status()
Pkg.update()
Pkg.build("CUDA")
Pkg.rm("Zygote")
Pkg.add("Zygote")
Pkg.build()
import GPUArrays
Pkg.add("FillArrays")
Pkg.add("DataFrames")

# FIXME I have to manually specify the version
# Pkg.update does not update the package
Pkg.add(Pkg.PackageSpec(name="Flux", version="0.11"))

#!!!
Pkg.add("Distributions")

##############################
## debugging (not needed)
Pkg.add("Colors")
Pkg.add("BenchmarkTools")
Pkg.add("GPUArrays")

Pkg.develop("Zygote")
Pkg.update()
Pkg.status()
Pkg.build(verbose=true)

##############################
## Using Flux#master and Zygote

Pkg.add("Flux")
Pkg.rm("Flux")
Pkg.rm("Tracker")
Pkg.add(PackageSpec(name="Flux", rev="master"))
Pkg.free("Flux")
Pkg.add(Pkg.PackageSpec(name="Zygote", rev="master"))

##############################
## Or using Flux 0.9 and Tracker
Pkg.rm("Flux")
Pkg.rm("GPUArrays")
Pkg.rm("Distributions")
Pkg.add(PackageSpec(name="Flux", version="0.9"))
Pkg.pin(PackageSpec(name="Flux", version="0.9"))
Pkg.add("Tracker")
Pkg.add("Distributions")

Pkg.add("HDF5")
Pkg.build("HDF5")

Pkg.add("PackageCompiler")

Pkg.Registry.update()

Pkg.instantiate()

Pkg.status()
Pkg.update()

# Pkg.build("Arpack")

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

"""
    share_rpc!(d::Dict, set::Dict)
"""
function share_rpc!(d::Dict, set::Dict)
    println("  Calculating regional purchase coefficient")
    cols = [:r,:g,:value]
    !(:ng   in collect(keys(d))) && _share_ng!(d, set)
    !(:d0   in collect(keys(d))) && _share_d0!(d, set)
    !(:mrt0 in collect(keys(d))) && _share_mrt0!(d)
    !(:xn0  in collect(keys(d))) && _share_xn0!(d, set)
    !(:mng  in collect(keys(d))) && _share_mn0!(d, set)

    df = dropnan(d[:d0] / (d[:d0] + d[:mn0]))

    df = edit_with(df, Drop(:g,"uti","=="))
    df_uti = fill_with((r = set[:r], g = "uti"), 0.9)

    d[:rpc] = sort([dropnan(df); df_uti])
end

"`ng`: Sectors not included in the CFS."
function _share_ng!(d::Dict, set::Dict)
    df = combine_over(copy(d[:cfs]), [:orig_state, :dest_state, :sg])
    df = fill_zero((g = set[:g],), df)
    df[!,:value] .= df[:,:value] .== 0.0
    d[:ng] = df
    set[:ng] = df[df[:,:value],:g]
end

"`d0`: Local supply-demand. Trade that remains within the same region."
function _share_d0!(d::Dict, set::Dict)
    cols = [:r,:g,:value]
    df = copy(d[:cfs])
    df = df[df[:,:orig_state] .== df[:,:dest_state],:]
    df = sort(edit_with(df, Rename(:orig_state,:r))[:,cols])
    d[:d0_temp] = copy(df)
    d[:d0] = _avg_ng(df, d, set)[:,cols]
end

"`mrt0`: Interstate trade (CFS)"
function _share_mrt0!(d::Dict)
    cols = [:r,:g,:value]
    df = copy(d[:cfs])
    df = df[df[:,:orig_state] .!= df[:,:dest_state],:]
    d[:mrt0] = df
end

"`mn0`: National demand (CFS)"
function _share_mn0!(d::Dict, set::Dict)
    :mrt0 in collect(keys(d)) && _share_mrt0!(d)

    cols = [:r,:g,:value]
    df = copy(d[:mrt0])
    df = edit_with(combine_over(df, :orig_state), Rename(:dest_state,:r))
    d[:mn0_temp] = copy(df)

    # Average and stuff.
    d[:mn0] = _avg_ng(df, d, set)[:,cols]
end

"`xn0`: National exports (CFS)"
function _share_xn0!(d::Dict, set::Dict)
    !(:mrt0 in collect(keys(d))) && _share_mrt0!(d)

    cols = [:r,:g,:value]
    df = copy(d[:mrt0])
    df = edit_with(combine_over(df, :dest_state), Rename(:orig_state,:r))
    d[:xn0_temp] = copy(df)

    d[:xn0] = _avg_ng(df, d, set)[:,cols]
end

" !!!! "
function _avg_ng(df::DataFrame, d::Dict, set::Dict)
    !(:ng in collect(keys(d))) && _share_ng!(d, set)

    df = copy(df)
    not_ng = edit_with(d[:ng], Replace(:value,Bool,Not))
    df = fill_zero((g = set[:g], r = set[:r]), df)

    df_ng = transform_over(df, :g) / transform_over(not_ng, :g)

    df = SLiDE._join_to_operate(df, df_ng, copy(d[:ng]))
    df = edit_with(df, Rename.([:x1,:x2,:x3], [:value, :value_ng, :ng]))
    df[!,:ng] .= convert_type.(Bool, df[:,:ng])

    df[!,:value_calc] .= df[:,:value]
    df[df[:,:ng], :value] .= df[df[:,:ng], :value_ng]
    return df
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
import Base
using Printf

@enum TokenType ILLEGAL = 1 EOF IDENT INT ASSIGN PLUS MINUS BANG ASTERISK SLASH LT GT COMMA SEMICOLON LPAREN RPAREN LBRACE RBRACE FUNCTION LET TRUE FALSE RETURN IF ELSE EQUAL NOTEQUAL LTE GTE

struct Token
  type::TokenType
  literal::String
  line::UInt32
  pos::UInt32
end

Token(type::TokenType, chr::Char, line::UInt32, pos::UInt32) =
  Token(type, String([chr]), line, pos)

function Base.show(io::IO, token::Token)
  str = @sprintf(
    """Token %s("%s") at line %d pos %d""",
    token.type,
    token.literal,
    token.line,
    token.pos
  )
  print(io, str)
end

max_stars_repo_name: viert/monkeylang-julia

max_issues_repo_name: viert/monkeylang-julia

max_forks_repo_name: viert/monkeylang-julia


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/bound.jl
# Computes interval bounds of various functions.
#############################################################################

@inline function _get_lo(m::GlobalOptimizer{N,T,S}, t::SAT) where {N,T<:AbstractFloat,S}
    c = t.coefficient
    i = t.variable_index.value
    xL = _lower_bound(FullVar, m, i)
    xU = _upper_bound(FullVar, m, i)
    return (c > 0.0) ? c*xL : c*xU
end

@inline function _get_bnds(m::GlobalOptimizer{N,T,S}, t::SAT) where {N,T<:AbstractFloat,S}
    c = t.coefficient
    i = t.variable_index.value
    xL = _lower_bound(FullVar, m, i)
    xU = _upper_bound(FullVar, m, i)
    l = (c > 0.0) ? c*xL : c*xU
    u = (c > 0.0) ? c*xU : c*xL
    return l, u
end

@inline function _get_bnds(m::GlobalOptimizer{N,T,S}, t::SQT) where {N,T<:AbstractFloat,S}
    c = t.coefficient
    v1 = t.variable_index_1.value
    v2 = t.variable_index_2.value
    xL1 = _lower_bound(FullVar, m, v1)
    xU1 = _upper_bound(FullVar, m, v1)

    if v1 == v2
        if c > zero(T)
            bnds = (zero(T) < xL1) ? 0.5*c*xL1*xL1 : ((xU1 <= zero(T)) ? 0.5*c*xU1*xU1 : zero(T))
        else
            bnds = (xL1 < xU1) ? 0.5*c*xU1*xU1 : 0.5*c*xL1*xL1
        end
    else
        xL2 = _lower_bound(FullVar, m, v2)
        xU2 = _upper_bound(FullVar, m, v2)
        bnds = c*Interval{Float64}(xL, xU)*Interval{Float64}(xL2, xU2)
    end
    return bnds
end

function lower_interval_bound(m::GlobalOptimizer{N,T,S}, t::Vector{SAT}, c) where {N,T<:AbstractFloat,S}
    mapreduce(i -> _get_lo(m, i), +, t; init = c)
end
function interval_bound(m::GlobalOptimizer{N,T,S}, t::Vector{SQT}, c) where {N,T<:AbstractFloat,S}
    mapreduce(i -> _get_bnds(m, i), +, t; init = c)
end

function interval_bound(m::GlobalOptimizer{N,T,S}, f::AffineFunction{T}) where {N,T<:Real,S}
    mapreduce(i -> _get_bnds(m, i), +, f.func.terms; init = (f.func.constant, f.func.constant))
end

function lower_interval_bound(m::GlobalOptimizer{N,T,S}, f::AffineFunction{T}) where {N,T<:Real,S}
    lower_interval_bound(m, f.func.terms, f.func.constant)
end

function interval_bound(m::GlobalOptimizer{N,T,S}, f::BufferedQuadraticEq) where {N,T<:AbstractFloat,S}
    bnds = lower_interval_bound(m, f.func.affine_terms, Interval{Float64}(f.func.constant))
    interval_bound(m, f.func.quadratic_terms, lower)
end
function lower_interval_bound(m::GlobalOptimizer{N,T,S}, f::BufferedQuadraticIneq) where {N,T<:AbstractFloat,S}
    bnds = interval_bound(m, f.func.quadratic_terms, Interval{Float64}(f.func.constant))
    return bnds.lo
end

###
### SECOND-ORDER CONE
###
function lower_interval_bound(m::GlobalOptimizer{N,T,S}, d::BufferedSOC) where {N,T<:AbstractFloat,S}
    vec_of_vi = d.variables.variables
    norm_bound = zero(Interval{T})
    for i = 2:length(vec_of_vi)
        vi = vec_of_vi[i].value
        x = Interval{T}(_lower_bound(FullVar, m, vi), _upper_bound(FullVar, m, vi))
        norm_bound += pow(x, 2)
    end
    norm_bound = sqrt(norm_bound)
    return norm_bound.lo - _upper_bound(FullVar, m, vec_of_vi[1].value)
end

###
### NONLINEAR FUNCTIONS
###
function lower_interval_bound(m::GlobalOptimizer{N,T,S}, d::BufferedNonlinearFunction{V}) where {N,T<:AbstractFloat,S,V}
    !d.has_value && forward_pass!(m._working_problem._relaxed_evaluator, d)
    ex = d.expr
    v = ex.isnumber[1] ? ex.numberstorage[1] : ex.setstorage[1].Intv.lo
    return v
end
function interval_bound(m::GlobalOptimizer{N,T,S}, d::BufferedNonlinearFunction{V}) where {N,T<:AbstractFloat,S,V}
    !d.has_value && forward_pass!(d.evaluator, d)
    ex = d.expr
    v = ex.isnumber[1] ? Interval(ex.numberstorage[1]) : ex.setstorage[1].Intv
    return v.lo, v.hi
end

function _is_feas(::Type{Val{:leq}}, m::GlobalOptimizer{N,T,S}, x::AffineFunction{T}) where {N,T<:Real,S}
    lower_interval_bound(m, x) <= zero(T)
end
function _is_feas(::Type{Val{:eq}}, m::GlobalOptimizer{N,T,S}, x::AffineFunction{T}) where {N,T<:Real,S}
    lower_value, upper_value = interval_bound(m, x)
    return lower_value <= zero(T) <= upper_value
end

_is_feas(m::GlobalOptimizer{N,T,S}, x::BufferedQuadraticIneq) where {N,T<:AbstractFloat,S} = lower_interval_bound(m, x) <= zero(T)
function _is_feas(m::GlobalOptimizer{N,T,S}, x::BufferedQuadraticEq) where {N,T<:AbstractFloat,S}
    lower_value, upper_value = interval_bound(m, x)
    return lower_value <= 0.0 <= upper_value
end
function _is_feas(m::GlobalOptimizer{N,T,S}, x::NonlinearExpression) where {N,T<:AbstractFloat,S}
    lower_value, upper_value = interval_bound(m, x)
    feasible_flag &= upper_value < x.lower_bound
    feasible_flag &= lower_value > x.upper_bound
    return !feasible_flag
end


text:
# For more information about how the error handling in the Julia wrapper is implemented,
# see the Internals section of the docs.

"""
    @check_warning

Check if the GRAMPC solver has triggered a warning and propagate it to the user as a
warning log message with the location information as the current line of code.
"""
macro check_warning()
    quote
        _check_warning( $(String(__source__.file)), $(__source__.line) )
    end
end

function _check_warning( file, line )
    # This buffer size matches the buffer in the code
    warnCString = Vector{UInt8}(undef, 1024)

    if( ccall( (:grampcjl_get_last_warning, libgrampc_julia), Cint, (Ptr{UInt8},), warnCString ) == 1 )
        # Convert the string
        warnString = unsafe_string( pointer( warnCString ) )

        @warn warnString _line=line _file=file
    end
end

"""
    _clear_warning()

Clears the current warning stored in the buffer inside the library.
"""
function _clear_warning()
    ccall( (:grampcjl_clear_last_warning, libgrampc_julia), Cvoid, () )
end


text:
import AIBECS: units, @units
import AIBECS: prior, @prior
import AIBECS: initial_value, @initial_value
import AIBECS: flattenable, @flattenable
import AIBECS: description, @description
@initial_value @units @description @prior @flattenable struct TestParameters{T} <: AbstractParameters{T}
    xgeo::T | 2.17 | u"mmol/m^3" | "Geological mean P concentration"             | LogNormal(0,1) | true
    τgeo::T |  1.0 | u"Myr"      | "Geological restoring timescale"              | LogNormal(0,1) | false
    k::T    | 10.0 | u"μmol/m^3" | "Half-saturation constant (Michaelis-Menten)" | LogNormal(0,1) | true
    z₀::T   | 80.0 | u"m"        | "Depth of the euphotic layer base"            | LogNormal(0,1) | false
    w₀::T   |  1.0 | u"m/d"      | "Sinking velocity at surface"                 | LogNormal(0,1) | true
    w′::T   | 0.22 | u"m/d/m"    | "Vertical gradient of sinking velocity"       |    Normal(0,1) | true
    τDOP::T | 0.25 | u"yr"       | "Remineralization timescale (DOP to DIP)"     | LogNormal(0,1) | true
    τPOP::T | 5.25 | u"d"        | "Dissolution timescale (POP to DOP)"          | LogNormal(0,1) | true
    σ::T    |  1/3 | unit(1)     | "Fraction of quick local uptake recycling"    |   Uniform(0,1) | true
    τDIP::T | 30.0 | u"d"        | "Uptake maximum timescale (DIP to POP)"       | LogNormal(0,1) | true
end
Distributions.gradlogpdf(::Uniform, x::T) where {T<:Real} = zero(T) # Required for these tests because undefined in Distributions...

p = TestParameters()

@testset "Parameters" begin
    println(p)
    @test p isa AIBECS.AbstractParameters{Float64}
    @testset "Adding a dual-valued vector" begin
        using DualNumbers
        dualp = p + ε * ones(sum(flattenable(p)))
        @test dualp isa AIBECS.AbstractParameters{Dual{Float64}}
        println(dualp)
    end
    @testset "Non-exported functions" begin
        @test AIBECS.symbols(p) isa NTuple{length(fieldnames(typeof(p))), Symbol}
        @test AIBECS.flattenable_symbols(p) isa NTuple{length(p), Symbol}
        @test length(values(p)) == length(fieldnames(typeof(p)))
        @test AIBECS.flattenable_values(p) isa Vector{Float64}
        @test length(AIBECS.flattenable_values(p)) == length(p)
    end
    @testset "Exported functions" begin
        @test values(p) isa Vector{Float64}
        @test vec(p) isa Vector{Float64}
        @test size(p) == size(vec(p))
        @test length(p) == length(vec(p))
        println(latex(p))
    end
    @testset "Unpacking" begin
        @unpack xgeo, τDIP = p
        @test xgeo == ustrip(upreferred(p.xgeo * units(p, :xgeo)))
        @test τDIP == ustrip(upreferred(p.τDIP * units(p, :τDIP)))
    end
    @testset "Prior for $spᵢ ($Dᵢ)" for (spᵢ, Dᵢ) in zip([:σ, :w′, :xgeo], [Uniform, Normal, LogNormal])
        D = prior(p, spᵢ)
        @test D isa Dᵢ
        λ2p, ∇λ2p, ∇²λ2p, p2λ = subfun(D), ∇subfun(D), ∇²subfun(D), invsubfun(D)
        @test p2λ(λ2p(0.5)) ≈ 0.5
        @test ForwardDiff.derivative(λ2p, 0.5) ≈ ∇λ2p(0.5)
        @test ForwardDiff.derivative(∇λ2p, 0.5) ≈ ∇²λ2p(0.5)
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                             Chron1.0Coupled.jl                                #
#                                                                               #
#     Illustrates the use of the Chron.jl package for eruption/deposition       #
#  age estimation and production of a stratigraphic age model.                  #
#                                                                               #
#     This file uses code cells (denoted by "## ---"). To evaluate a cell in    #
#  the Julia REPL and move to the next cell, the default shortcut in Atom is    #
#  alt-shift-enter.                                                             #
#                                                                               #
#      You may have to adjust the path below which specifies the location of    #
#  the CSV data files for each sample, depending on what you want to run.       #
#                                                                               #
#   Last modified by C. Brenhin Keller 2020-05-25                               #
#                                                                               #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## --- Load the Chron package - - - - - - - - - - - - - - - - - - - - - - - - -

    try
        using Chron
    catch
        using Pkg
        Pkg.add("Chron")
        using Chron
    end

    using Statistics, StatsBase, DelimitedFiles, SpecialFunctions
    using Plots; gr();

## --- Define sample properties - - - - - - - - - - - - - - - - - - - - - - - -

    # This example data is from Clyde et al. (2016) "Direct high-precision
    # U–Pb geochronology of the end-Cretaceous extinction and calibration of
    # Paleocene astronomical timescales" EPSL 452, 272–280.
    # doi: 10.1016/j.epsl.2016.07.041

    nSamples = 5 # The number of samples you have data for
    smpl = NewChronAgeData(nSamples)
    smpl.Name      =  ("KJ08-157", "KJ04-75", "KJ09-66",  "KJ4-72", "KJ04-70",)
    smpl.Height   .=  [     -52.0,      44.0,      54.0,      82.0,      93.0,]
    smpl.Height_sigma .= [    3.0,       1.0,       3.0,       3.0,       3.0,]
    smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
    smpl.Path = "examples/DenverUPbExampleData/" # Where are the data files?
    smpl.inputSigmaLevel = 2 # i.e., are the data files 1-sigma or 2-sigma. Integer.
    smpl.Age_Unit = "Ma" # Unit of measurement for ages and errors in the data files
    smpl.Height_Unit = "cm" # Unit of measurement for Height and Height_sigma

    # IMPORTANT: smpl.Height must increase with increasing stratigraphic height
    # -- i.e., stratigraphically younger samples must be more positive. For this
    # reason, it is convenient to represent depths below surface as negative
    # numbers.

    # For each sample in smpl.Name, we must have a csv file at smpl.Path which
    # contains each individual mineral age and uncertainty. For instance,
    # examples/DenverUPbExampleData/KJ08-157.csv contains:
    #
    #   66.12,0.14
    #   66.115,0.048
    #   66.11,0.1
    #   66.11,0.17
    #   66.096,0.056
    #   66.088,0.081
    #   66.085,0.076
    #   66.073,0.084
    #   66.07,0.11
    #   66.055,0.043
    #   66.05,0.16
    #   65.97,0.12

## --- Bootstrap pre-eruptive distribution - - - - - - - - - - - - - - - - - - -

    # Bootstrap a KDE of the pre-eruptive (or pre-depositional) mineral age
    # distribution using a KDE of stacked sample data from each data file
    BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)
    x = range(0,1,length=length(BootstrappedDistribution))
    h = plot(x, BootstrappedDistribution, label="Bootstrapped distribution",
         xlabel="Time (arbitrary units)", ylabel="Probability Density", fg_color_legend=:white)
    savefig(h, joinpath(smpl.Path,"BootstrappedDistribution.pdf"))
    display(h)

## --- Estimate the eruption age distributions for each sample  - - - - - - - -

    # Configure distribution model here
    distSteps = 5*10^5 # Number of steps to run in distribution MCMC
    distBurnin = floor(Int,distSteps/2) # Number to discard

    # Choose the form of the prior distribution to use
    # Some pre-defined possiblilities include UniformDisribution,
    # TriangularDistribution, and MeltsVolcanicZirconDistribution
    # or you can define your own as with BootstrappedDistribution above
    dist = BootstrappedDistribution

    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    # Run MCMC to estimate saturation and eruption/deposition age distributions
    @time tMinDistMetropolis(smpl,distSteps,distBurnin,dist)

    # This will save rank-order and distribution plots, and print results to a
    # csv file -- you can find them in smpl.Path

    # # (Optional) Save the sample struct for later use
    # using JLD: @save, @load
    # @save "smpl.jld" smpl

## --- Run stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - - -

    # # (Optional) Load the saved sample struct
    # @load "smpl.jld" smpl

    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    # Configure the stratigraphic Monte Carlo model
    config = NewStratAgeModelConfiguration()
    # If you in doubt, you can probably leave these parameters as-is
    config.resolution = 1.0 # Same units as sample height. Smaller is slower!
    config.bounding = 0.5 # how far away do we place runaway bounds, as a fraction of total section height
    (bottom, top) = extrema(smpl.Height)
    npoints_approx = round(Int,length(bottom:config.resolution:top) * (1 + 2*config.bounding))
    config.nsteps = 15000 # Number of steps to run in distribution MCMC
    config.burnin = 10000*npoints_approx # Number to discard
    config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    # Run the stratigraphic MCMC model
    @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)

    # # Other options:
    # # Youngest Zircon
    # for i=1:length(smpl.Name)
    #     data = readdlm(string(smpl.Path, smpl.Name[i], ".csv"),',')
    #     Iyz = argmin(data[:,1])
    #     smpl.Age[i] = minimum(data[Iyz,1])
    #     smpl.Age_sigma[i] = minimum(data[Iyz,2]/smpl.inputSigmaLevel)
    # end
    # (mdl, agedist, lldist) = StratMetropolis(smpl, config)

    # # Low-N Weighted Mean
    # for i=1:length(smpl.Name)
    #     data = readdlm(string(smpl.Path, smpl.Name[i], ".csv"),',')
    #     sI = sortperm(data[:,1])
    #     # Weighted mean of youngst 3 zircons per sample (assuming there are at least 3 zircons in sample)
    #     Ns = min(size(data,1),3)
    #     (mu, sigma) = awmean(data[sI[1:Ns],1],data[sI[1:Ns],2]./smpl.inputSigmaLevel)
    #     smpl.Age[i] = mu
    #     smpl.Age_sigma[i] = sigma
    # end
    # (mdl, agedist, lldist) = StratMetropolis(smpl, config)

## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - -

    # Plot results (mean and 95% confidence interval for both model and data)
    hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(round(Int,minimum(mdl.Height)),0.5,:blue), label="model")
    plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white) # Center line
    t = smpl.Age_Sidedness .== 0 # Two-sided constraints (plot in black)
    any(t) && plot!(hdl, smpl.Age[t], smpl.Height[t], xerror=(smpl.Age[t]-smpl.Age_025CI[t],smpl.Age_975CI[t]-smpl.Age[t]),label="data",seriestype=:scatter,color=:black)
    t = smpl.Age_Sidedness .== 1 # Minimum ages (plot in cyan)
    any(t) && plot!(hdl, smpl.Age[t], smpl.Height[t], xerror=(smpl.Age[t]-smpl.Age_025CI[t],zeros(count(t))),label="",seriestype=:scatter,color=:cyan,msc=:cyan)
    any(t) && zip(smpl.Age[t], smpl.Age[t].+nanmean(smpl.Age_sigma[t])*4, smpl.Height[t]) .|> x-> plot!([x[1],x[2]],[x[3],x[3]], arrow=true, label="", color=:cyan)
    t = smpl.Age_Sidedness .== -1 # Maximum ages (plot in orange)
    any(t) && plot!(hdl, smpl.Age[t], smpl.Height[t], xerror=(zeros(count(t)),smpl.Age_975CI[t]-smpl.Age[t]),label="",seriestype=:scatter,color=:orange,msc=:orange)
    any(t) && zip(smpl.Age[t], smpl.Age[t].-nanmean(smpl.Age_sigma[t])*4, smpl.Height[t]) .|> x-> plot!([x[1],x[2]],[x[3],x[3]], arrow=true, label="", color=:orange)
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))")
    savefig(hdl,"AgeDepthModel.pdf")
    display(hdl)

## --- Interpolate model age at a specific stratigraphic height - - - - - - - -

    # Stratigraphic height at which to interpolate
    interp_height = 0

    age_at_height = linterp1s(mdl.Height,mdl.Age,interp_height)
    age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,interp_height)
    age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,interp_height)
    print("Interpolated age at height=$interp_height: $age_at_height +$(age_at_height_max-age_at_height)/-$(age_at_height-age_at_height_min) $(smpl.Age_Unit)")

    # Optional: interpolate full age distribution
    interpolated_distribution = Array{Float64}(undef,size(agedist,2))
    for i=1:size(agedist,2)
        interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],interp_height)
    end
    hdl = histogram(interpolated_distribution, nbins=50, label="")
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit)) at height=$interp_height", ylabel="Likelihood (unnormalized)")
    savefig(hdl, "Interpolated age distribution.pdf")
    display(hdl)

## --- Calculate deposition rate binned by age  - - - - - - - - - - - - - - - -

    # Set bin width and spacing
    binwidth = round(nanrange(mdl.Age)/10,sigdigits=1) # Can also set manually, commented out below
    # binwidth = 0.01 # Same units as smpl.Age
    binoverlap = 10
    ages = collect(minimum(mdl.Age):binwidth/binoverlap:maximum(mdl.Age))
    bincenters = ages[1+Int(binoverlap/2):end-Int(binoverlap/2)]
    spacing = binoverlap

    # Calculate rates for the stratigraphy of each markov chain step
    dhdt_dist = Array{Float64}(undef, length(ages)-binoverlap, config.nsteps)
    @time for i=1:config.nsteps
        heights = linterp1(reverse(agedist[:,i]), reverse(mdl.Height), ages)
        dhdt_dist[:,i] .= abs.(heights[1:end-spacing] - heights[spacing+1:end]) ./ binwidth
    end

    # # Exact (added precision is below sampling resolution, so not useful) and very slow
    # @showprogress "Calculating dh/dt..." for i=1:config.nsteps
    #     for j=1:length(bincenters)
    #         t = (agedist[:,i] .> bincenters[j] - binwidth/2) .& (agedist[:,i] .< bincenters[j] + binwidth/2)
    #         dhdt_dist[j,i] = sum(t) * resolution / binwidth
    #     end
    # end

    # Find mean and 1-sigma (68%) CI
    dhdt = nanmean(dhdt_dist,dim=2)
    dhdt_50p = nanmedian(dhdt_dist,dim=2)
    dhdt_16p = nanpctile(dhdt_dist,15.865,dim=2) # Minus 1-sigma (15.865th percentile)
    dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)

    # Plot results
    hdl = plot(bincenters,dhdt, label="Mean", color=:black, linewidth=2)
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.4,:darkblue), linealpha=0, label="68% CI")
    plot!(hdl,bincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Accumulation Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
    savefig(hdl,"AccumulationRateModel.pdf")
    display(hdl)

## --- Multiple confidence intervals (every 10) - - - - - - - - - - - - - - - -

    dhdt_20p = nanpctile(dhdt_dist,20,dim=2)
    dhdt_80p = nanpctile(dhdt_dist,80,dim=2)
    dhdt_25p = nanpctile(dhdt_dist,25,dim=2)
    dhdt_75p = nanpctile(dhdt_dist,75,dim=2)
    dhdt_30p = nanpctile(dhdt_dist,30,dim=2)
    dhdt_70p = nanpctile(dhdt_dist,70,dim=2)
    dhdt_35p = nanpctile(dhdt_dist,35,dim=2)
    dhdt_65p = nanpctile(dhdt_dist,65,dim=2)
    dhdt_40p = nanpctile(dhdt_dist,40,dim=2)
    dhdt_60p = nanpctile(dhdt_dist,60,dim=2)
    dhdt_45p = nanpctile(dhdt_dist,45,dim=2)
    dhdt_55p = nanpctile(dhdt_dist,55,dim=2)

    # Plot results
    hdl = plot(bincenters,dhdt, label="Mean", color=:black, linewidth=2)
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_20p; reverse(dhdt_80p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_25p; reverse(dhdt_75p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_30p; reverse(dhdt_70p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_35p; reverse(dhdt_65p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_40p; reverse(dhdt_60p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_45p; reverse(dhdt_55p)], fill=(0,0.2,:darkblue), linealpha=0, label="")
    plot!(hdl,bincenters,dhdt_50p, label="Median", color=:grey, linewidth=1)
    plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Accumulation Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white)
    savefig(hdl,"AccumulationRateModelCI.pdf")
    display(hdl)

## --- Make heatmap - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    rateplotmax = 3*maximum(dhdt) # May want to adjust this -- this is just a guess
    using StatsBase: fit, Histogram
    edges = range(0, rateplotmax, length=length(ages)-spacing+1)
    dhdt_im = Array{Float64}(undef,length(ages)-spacing,length(ages)-spacing)
    for i=1:length(ages)-spacing
        dhdt_im[:,i] .= fit(Histogram, dhdt_dist[i, .~ isnan.(view(dhdt_dist,i,:))], edges, closed=:left).weights
    end

    # Rescale image to fit in UInt8 (0-255)
    using IndirectArrays: IndirectArray
    imSc = dhdt_im./nanpctile(dhdt_im[:],97.5)*256 # Rescale to include 97.5 CI (may need ot adjust)
    imSc[imSc.>255] .= 255 # Cut off to fit in Uint8

    # Apply colormap. Available colormaps include viridis, inferno, plasma, fire
    A = IndirectArray(floor.(UInt8,imSc) .+ 1, inferno)

    # Plot image
    scale=1e-5
    img = plot(bincenters,cntr(edges)*scale,A,yflip=false,xflip=false, colorbar=:right)
    plot!(img, xlabel="Age ($(smpl.Age_Unit))", ylabel="Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) * $scale, $binwidth $(smpl.Age_Unit) Bin)")
    savefig(img,"AccumulationRateModelHeatmap.pdf")
    display(img)

    # dhdt_im_log = copy(dhdt_im)
    # dhdt_im_log[dhdt_im .>0] = log10.(dhdt_im[dhdt_im .>0])
    # heatmap(bincenters,cntr(edges),dhdt_im_log, xlabel="Age ($smpl.Age_Unit)", ylabel="Rate ($smpl.Height_Unit / $smpl.Age_Unit, $binwidth $smpl.Age_Unit Bin)")

## --- Probability that a given interval of stratigraphy was deposited entirely before/after a given time

    # Stratigraphic height and absoltue age/uncert to test
    testHeight = -40.0
    testAge = 66.0
    testAge_sigma = 0.05

    # Find index of nearest model height node
    nearest = argmin((testHeight .- mdl.Height).^2)

    # Cycle through each possible age within testAge +/- 5 sigma, with resolution of 1/50 sigma
    test_ages = (testAge-5*testAge_sigma):testAge_sigma/50:(testAge+5*testAge_sigma)
    test_prob_older = Array{Float64}(undef,size(test_ages))
    # Evaluate the probability that model age is older than each test_age at the given strat level
    for i=1:length(test_ages)
        test_prob_older[i] = sum(agedist[nearest,:] .> test_ages[i]) ./ size(agedist,2)
    end

    # Normalized probability for each distance away from testAge between +5sigma and -5sigma
    prob_norm = normpdf.(testAge, testAge_sigma, test_ages) ./ sum(normpdf.(testAge, testAge_sigma, test_ages));  # SUM = 1

    # Integrate the product
    prob_older = sum(test_prob_older .* prob_norm)
    print("$(prob_older*100) % chance that $(mdl.Height[nearest]) $(smpl.Height_Unit) was deposited before $testAge +/- $testAge_sigma $(smpl.Age_Unit) Gaussian")

## --- (Optional) If your section has hiata / exposure surfaces of known duration, try this:

    # Data about hiatuses
    nHiatuses = 2 # The number of hiatuses you have data for
    hiatus = NewHiatusData(nHiatuses) # Struct to hold data
    hiatus.Height         = [-7.0, 35.0 ]
    hiatus.Height_sigma   = [ 0.0,  0.0 ]
    hiatus.Duration       = [ 0.3,  0.3 ]
    hiatus.Duration_sigma = [ 0.05, 0.05]

    # Run the model. Note the additional `hiatus` arguments
    @time (mdl, agedist, hiatusdist, lldist) = StratMetropolisDist(smpl, hiatus, config)

    # Plot results (mean and 95% confidence interval for both model and data)
    hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model")
    plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white)
    plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black)
    plot!(hdl, xlabel="Age (Ma)", ylabel="Height (cm)")

## --- (Optional) Add systematic uncertainties for U-Pb data

    # # Tracer (ET2535) uncertainty converted from per cent to relative
    # unc_tracer = 0.03/2/100
    #
    # # U-238 Decay constant and uncertainty, Myr^-1
    # lambda238 = 1.55125e-10 * 1e6
    # unc_lambda238 = 0.107/2/100 # converted from per cent to relative
    #
    # # Consider only the distribution of ages at model nodes where we have an ash bed
    # age_dist_X = Array{Float64}(undef,length(smpl.Height),size(agedist,2))
    # for i = 1:length(smpl.Height)
    #    closest_model_node = argmin(abs.(mdl.Height-smpl.Height[i]))
    #    age_dist_X[i,:] = agedist[closest_model_node,:]
    # end
    #
    # # Convert ages to 206Pb/238U ratios of the distribution
    # ratio_dist = exp.(age_dist_X.*lambda238)-1
    #
    # # Add tracer uncertainty
    # ratio_dist_tracerunc = Array{Float64}(undef,size(ratio_dist))
    # for i=1:size(ratio_dist,2)
    #     ratio_dist_tracerunc[:,i] = ratio_dist[:,i].*(1 + unc_tracer*randn())
    # end
    #
    # # Convert 206/238 ratios back to ages, in Ma
    # age_dist_XY = log.(ratio_dist_tracerunc+1)./lambda238
    #
    # # Add decay constant uncertainty
    # age_dist_XYZ = Array{Float64}(undef,size(ratio_dist))
    # for i=1:size(ratio_dist,2)
    #     age_dist_XYZ[:,i] = log.(ratio_dist_tracerunc[:,i]+1)./(lambda238.*(1 + unc_lambda238.*randn()))
    # end
    #
    # # Calculate the means and 95% confidence intervals for different levels of systematic uncertainties
    #
    # age_dist_X_mean = mean(age_dist_X,2) # Mean age
    # age_dist_X_std =  std(age_dist_X,2) # Standard deviation
    # age_dist_X_median = median(age_dist_X,2) # Median age
    # age_dist_X_025p = nanpctile(age_dist_X,2.5,dim=2) # 2.5th percentile
    # age_dist_X_975p = nanpctile(age_dist_X,97.5,dim=2) # 97.5th percentile
    #
    # age_dist_XY_mean = mean(age_dist_XY,2) # Mean age
    # age_dist_XY_std =  std(age_dist_XY,2) # Standard deviation
    # age_dist_XY_median = median(age_dist_XY,2) # Median age
    # age_dist_XY_025p = nanpctile(age_dist_XY,2.5,dim=2) # 2.5th percentile
    # age_dist_XY_975p = nanpctile(age_dist_XY,97.5,dim=2) # 97.5th percentile
    #
    # age_dist_XYZ_mean = mean(age_dist_XYZ,2) # Mean age
    # age_dist_XYZ_std =  std(age_dist_XYZ,2) # Standard deviation
    # age_dist_XYZ_median = median(age_dist_XYZ,2) # Median age
    # age_dist_XYZ_025p = nanpctile(age_dist_XYZ,2.5,dim=2) # 2.5th percentile
    # age_dist_XYZ_975p = nanpctile(age_dist_XYZ,97.5,dim=2) # 97.5th percentile
    #
    # age_X_95p = [age_dist_X_mean age_dist_X_975p-age_dist_X_mean age_dist_X_mean-age_dist_X_025p]
    # age_XY_95p = [age_dist_XY_mean age_dist_XY_975p-age_dist_XY_mean age_dist_XY_mean-age_dist_XY_025p]
    # age_XYZ_95p = [age_dist_XYZ_mean age_dist_XYZ_975p-age_dist_XYZ_mean age_dist_XYZ_mean-age_dist_XYZ_025p]

## ---


text:
module DimensionalData

# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end DimensionalData

# Standard lib
using Dates,
      LinearAlgebra,
      Random,
      Statistics,
      SparseArrays

using Base.Broadcast: Broadcasted, BroadcastStyle, DefaultArrayStyle, AbstractArrayStyle,
      Unknown

using Base: tail, OneTo, @propagate_inbounds

# Ecosystem
import Adapt,
       ArrayInterface,
       ConstructionBase,
       RecipesBase,
       Tables

using RecipesBase: @recipe

export Dimension, IndependentDim, DependentDim, XDim, YDim, ZDim, TimeDim,
       X, Y, Z, Ti, ParametricDimension, Dim, AnonDim, Coord

export Selector, At, Between, Contains, Near, Where

export Locus, Center, Start, End, AutoLocus

export Order, Ordered, Unordered, UnknownOrder, AutoOrder

export IndexOrder, ArrayOrder, Relation,
       ForwardIndex, ReverseIndex, UnorderedIndex,
       ForwardArray, ReverseArray, ForwardRelation, ReverseRelation

export Sampling, Points, Intervals

export Span, Regular, Irregular, Explicit, AutoSpan

export Mode, IndexMode, Auto, AutoMode, NoIndex

export Aligned, AbstractSampled, Sampled, AbstractCategorical, Categorical

export Unaligned, Transformed

export AbstractDimArray, DimArray, AbstractDimensionalArray, DimensionalArray

export AbstractDimTable, DimTable

export AbstractDimStack, DimStack

export DimIndices, DimKeys

export AbstractMetadata, Metadata, NoMetadata

export AbstractName, Name, NoName

export data, dims, refdims, mode, metadata, name, label, units,
       val, index, order, sampling, span, bounds, locus, <|

export dimnum, hasdim, otherdims, commondims, setdims, swapdims, sortdims,
       set, rebuild, reorder, modify, dimwise, dimwise!

export order, indexorder, arrayorder, relation

export @dim

const DD = DimensionalData

include("interface.jl")
include("mode.jl")
include("metadata.jl")
include("name.jl")
include("identify.jl")
include("dimension.jl")
include("array.jl")
include("stack.jl")
include("tables.jl")
include("selector.jl")
include("dimindices.jl")
include("indexing.jl")
include("primitives.jl")
include("broadcast.jl")
include("methods.jl")
include("matmul.jl")
include("coord.jl")
include("set.jl")
include("utils.jl")
include("plotrecipes.jl")
include("show.jl")

# For compat with old versions
const AbstractDimensionalArray = AbstractDimArray
const DimensionalArray = DimArray

const AbstractDimDataset = AbstractDimStack
const DimDataset = DimStack

function _precompile()
    precompile(DimArray, (Array{Int,1}, typeof(X())))
    precompile(DimArray, (Array{Int,2}, Tuple{typeof(X()), typeof(Y())}))
    precompile(DimArray, (Array{Float64,1}, typeof(X())))
    precompile(DimArray, (Array{Float64,2}, Tuple{typeof(X()), typeof(Y())}))
    precompile(DimArray, (Array{Float64,3}, Tuple{typeof(X()), typeof(Y()), typeof(Z())}))
end

_precompile()

end


text:
"""
   Joint{T}

   Abstract type for 3-dimensional constraint between two Body objects
"""
abstract type Joint{T,Nλ,Nb,N,Nb½} end

# joints
function joint_constraint(joint::Joint{T},
    xa::AbstractVector, qa::Quaternion, xb::AbstractVector, qb::Quaternion, η) where {T}
    return constraint_mask(joint) * displacement(joint, xa, qa, xb, qb)
end

function joint_constraint_jacobian_configuration(relative::Symbol, joint::Joint{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion, η) where {T}
    X, Q = displacement_jacobian_configuration(relative, joint, xa, qa, xb, qb, attjac=false)
    return constraint_mask(joint) * [X Q]
end

# constraints
function constraint(joint::Joint{T,Nλ,0},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion, η, μ) where {T,Nλ}
    joint_constraint(joint, xa, qa, xb, qb, η)
end

constraint(joint::Joint, pbody::Node, cbody::Node, λ, μ, timestep) = constraint(joint, next_configuration(pbody.state, timestep)..., next_configuration(cbody.state, timestep)..., λ, μ)

# constraint Jacobians
function constraint_jacobian(joint::Joint{T,Nλ,0}, η) where {T,Nλ}
    return Diagonal(REG * sones(T,Nλ))
end

function constraint_jacobian(joint::Joint{T,Nλ,Nb}, η) where {T,Nλ,Nb}
    s, γ = split_impulses(joint, η)
    c1 = [Diagonal(γ + REG * sones(T, Nb)); Diagonal(sones(Nb)); szeros(Nλ, Nb)]
    c2 = [Diagonal(s + REG * sones(T, Nb)); szeros(Nb, Nb); szeros(Nλ, Nb)]
    c3 = [szeros(Nb, Nλ); szeros(Nb, Nλ); Diagonal(REG * sones(T, Nλ))]
    return [c1 c2 c3]
end

constraint_jacobian_configuration(relative::Symbol, joint::Joint, pbody::Node, cbody::Node, λ, timestep) = constraint_jacobian_configuration(relative, joint, next_configuration(pbody.state, timestep)..., next_configuration(cbody.state, timestep)..., λ)

function constraint_jacobian_configuration(relative::Symbol, joint::Joint{T,Nλ,0},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion, η) where {T,Nλ}
    joint_constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
end

# masks
constraint_mask(::Joint{T,0}) where T = szeros(T,0,3)
constraint_mask(joint::Joint{T,1}) where T = joint.axis_mask3
constraint_mask(joint::Joint{T,2}) where T = [joint.axis_mask1; joint.axis_mask2]
constraint_mask(::Joint{T,3}) where T = SMatrix{3,3,T,9}(I)

nullspace_mask(::Joint{T,0}) where T = SMatrix{3,3,T,9}(I)
nullspace_mask(joint::Joint{T,1}) where T = [joint.axis_mask1; joint.axis_mask2]
nullspace_mask(joint::Joint{T,2}) where T = joint.axis_mask3
nullspace_mask(::Joint{T,3}) where T = szeros(T,0,3)

# impulse maps
impulse_map(relative::Symbol, joint::Joint, pbody::Node, cbody::Node, λ) = impulse_map(relative, joint, current_configuration(pbody.state)..., current_configuration(cbody.state)..., λ)

function impulse_map(relative::Symbol, joint::Joint{T,Nλ,Nb},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion,
        η) where {T,Nλ,Nb}
    return impulse_transform(relative, joint, xa, qa, xb, qb) * impulse_projector(joint)
end

function impulse_map(relative::Symbol, joint::Joint{T,Nλ,0},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    η) where {T,Nλ}
    J = constraint_jacobian_configuration(relative, joint, xa, qa, xb, qb, η)
    G = cat(Diagonal(sones(3)), LVᵀmat(relative == :parent ? qa : qb), dims=(1,2))
    return Diagonal([sones(3); 0.5 * sones(3)]) * transpose(J * G)
end

function impulse_projector(joint::Joint{T,Nλ,0}) where {T,Nλ}
    zerodimstaticadjoint(constraint_mask(joint))
end

function impulse_projector(joint::Joint{T,Nλ,Nb}) where {T,Nλ,Nb}
    zerodimstaticadjoint([szeros(Nb,3); -nullspace_mask(joint); nullspace_mask(joint); constraint_mask(joint)])
end

# inputs
function set_input!(joint::Joint, input::SVector)
    joint.input = zerodimstaticadjoint(nullspace_mask(joint)) * input
    return
end

set_input!(joint::Joint) = return

function add_input!(joint::Joint, input::SVector)
    joint.input += zerodimstaticadjoint(nullspace_mask(joint)) * input
    return
end

add_input!(joint::Joint) = return

function input_jacobian_control(relative::Symbol, joint::Joint, pbody::Node, cbody::Node)
    return input_jacobian_control(relative, joint, current_configuration(pbody.state)..., current_configuration(cbody.state)...) * zerodimstaticadjoint(nullspace_mask(joint))
end

# minimal coordinates
minimal_coordinates(joint::Joint{T,Nλ}) where {T,Nλ} = szeros(T, 3 - Nλ)

function minimal_coordinates(joint::Joint, pbody::Node, cbody::Node)
    return minimal_coordinates(joint, current_configuration(pbody.state)..., current_configuration(cbody.state)...)
end

function minimal_velocities(joint::Joint, pnode::Node, cnode::Node, timestep)
	minimal_velocities(joint, initial_configuration_velocity(pnode.state)...,
		initial_configuration_velocity(cnode.state)..., timestep)
end

################################################################################
# Utilities
################################################################################
Base.length(joint::Joint{T,Nλ}) where {T,Nλ} = Nλ
Base.zero(joint::Joint{T,Nλ}) where {T,Nλ} = szeros(T, Nλ, 6)
joint_length(joint::Joint{T,Nλ}) where {T,Nλ} = Nλ
limits_length(joint::Joint{T,Nλ,Nb}) where {T,Nλ,Nb} = Nb
impulses_length(joint::Joint{T,Nλ,Nb,N}) where {T,Nλ,Nb,N} = N

function split_impulses(joint::Joint{T,Nλ,Nb}, η) where {T,Nλ,Nb}
    s = η[SVector{Nb,Int}(1:Nb)]
    γ = η[SVector{Nb,Int}(Nb .+ (1:Nb))]
    return s, γ
end

function joint_impulse_index(joint::Joint{T,Nλ,Nb,N}, s::Int) where {T,Nλ,Nb,N}
    ind = SVector{N,Int}(s+1:s+N)
    return ind
end

input_dimension(joint::Joint{T,N}) where {T,N} = 3 - N