text:
########################
# Do not use for now
######################
#####################
# Expectation Flow circuits
# For use of algorithms depending on pairs of nodes of two circuits
#####################

"A expectation circuit node that has pair of origins of type PC and type LC"
abstract type DecoratorΔNodePair{PC<:ΔNode, LC<:ΔNode} <: ΔNode end

abstract type ExpFlowΔNode{PC, LC, F} <: DecoratorΔNodePair{PC, LC} end

const ExpFlowΔ{O} = AbstractVector{<:ExpFlowΔNode{<:O}}

struct UpExpFlow{PC, LC, F} <: ExpFlowΔNode{PC, LC, F}
    p_origin::PC
    f_origin::LC
    children::Vector{<:ExpFlowΔNode{<:PC, <:LC, <:F}}
    f::F
    fg::F
end

"""
Construct a upward expectation flow circuit from a given pair of PC and LC circuits
Note that its assuming the two circuits share the same vtree
"""
function ExpFlowΔ(pc::ProbΔ, lc::LogisticΔ, batch_size::Int, ::Type{El}) where El
    pc_type = grapheltype(pc)
    lc_type = grapheltype(lc)

    F = Array{El, 2}
    fmem = () -> zeros(1, batch_size) #Vector{El}(undef, batch_size)  #some_vector(El, batch_size) # note: fmem's return type will determine type of all UpFlows in the circuit (should be El)
    fgmem = () -> zeros(classes(lc[end]), batch_size)

    root_pc = pc[end]
    root_lc = lc[end- 1]

    cache = Dict{Pair{ΔNode, ΔNode}, ExpFlowΔNode}()
    sizehint!(cache, (length(pc) + length(lc))*4÷3)
    expFlowCircuit = Vector{ExpFlowΔNode}()

    function ExpflowTraverse(n::Prob⋁, m::Logistic⋁)
        get!(cache, Pair(n, m)) do
            children = [ ExpflowTraverse(i, j) for i in n.children for j in m.children]
            node = UpExpFlow{pc_type,lc_type, F}(n, m, children, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end
    function ExpflowTraverse(n::Prob⋀, m::Logistic⋀)
        get!(cache, Pair(n, m)) do
            children = [ ExpflowTraverse(z[1], z[2]) for z in zip(n.children, m.children) ]
            node = UpExpFlow{pc_type,lc_type, F}(n, m, children, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end
    function ExpflowTraverse(n::ProbLiteral, m::Logistic⋁)
        get!(cache, Pair(n, m)) do
            children = Vector{ExpFlowΔNode{pc_type,lc_type, F}}() # TODO
            node = UpExpFlow{pc_type,lc_type, F}(n, m, children, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end
    function ExpflowTraverse(n::ProbLiteral, m::LogisticLiteral)
        get!(cache, Pair(n, m)) do
            children = Vector{ExpFlowΔNode{pc_type,lc_type, F}}() # TODO
            node = UpExpFlow{pc_type,lc_type, F}(n, m, children, fmem(), fgmem())
            push!(expFlowCircuit, node)
            return node
        end
    end

    ExpflowTraverse(root_pc, root_lc)
    expFlowCircuit
end

function exp_pass_up(pc::ProbΔ, lc::LogisticΔ, data::XData{E}) where{E <: eltype(F)} where{PC, LC, F}
    expFlowCircuit = ExpFlowΔ(pc, lc, num_examples(data), Float64);
    for n in expFlowCircuit
        exp_pass_up_node(n, data)
    end
    expFlowCircuit
end

function exp_pass_up(fc::ExpFlowΔ, data::XData{E}) where{E <: eltype(F)} where{PC, LC, F}
    #TODO write resize_flows similar to flow_circuits
    #     and give as input the expFlowCircuit instead
    #expFlowCircuit = ExpFlowΔ(pc, lc, num_examples(data), Float64);
    for n in fc
        exp_pass_up_node(n, data)
    end
end

function exp_pass_up_node(node::ExpFlowΔNode{PC,LC,F}, data::XData{E}) where{E <: eltype(F)} where{PC, LC, F}
    pType = typeof(node.p_origin)
    fType = typeof(node.f_origin)

    if node.p_origin isa Prob⋁ && node.f_origin isa Logistic⋁
        #todo this ordering might be different than the ExpFlowΔNode children
        pthetas = [exp(node.p_origin.log_thetas[i])
                    for i in 1:length(node.p_origin.children) for j in 1:length(node.f_origin.children)]
        fthetas = [node.f_origin.thetas[j,:] # only taking the first class for now
            for i in 1:length(node.p_origin.children) for j in 1:length(node.f_origin.children)]

        node.f .= 0.0
        node.fg .= 0.0
        for z = 1:length(node.children)
            node.f  .+= pthetas[z] .* node.children[z].f
            node.fg .+= (pthetas[z] .* fthetas[z]) .* node.children[z].f
            node.fg .+= pthetas[z] .* node.children[z].fg
        end
    elseif node.p_origin isa Prob⋀ && node.f_origin isa Logistic⋀
        node.f .= node.children[1].f .* node.children[2].f # assume 2 children
        node.fg .= (node.children[1].f .* node.children[2].fg) .+
                   (node.children[2].f .* node.children[1].fg)

    elseif node.p_origin isa ProbLiteral
        if node.f_origin isa Logistic⋁
            m = node.f_origin.children[1]
        elseif node.f_origin isa LogisticLiteral
            m = node.f_origin
        else
            error("Invalid Types of pairs {$pType} - {$fType}")
        end

        var = lit2var(literal(m))
        X = feature_matrix(data)
        if positive(node.p_origin) && positive(m)
            node.f[:, X[:, var] .!= 0 ] .= 1.0 # positive and missing observations
            node.f[:, X[:, var] .== 0 ] .= 0.0
        elseif negative(node.p_origin) && negative(m)
            node.f[:, X[:, var] .!= 1 ] .= 1.0 # negative and missing observations
            node.f[:, X[:, var] .== 1 ] .= 0.0
        else
            node.f .= 0.0
        end

        if node.f_origin isa Logistic⋁
            node.fg .= node.f .* transpose(node.f_origin.thetas)
        else
            node.fg .= 0.0
        end

    else
        error("Invalid Types of pairs {$pType} - {$fType}")
    end

end


text:
"""
    BenchmarkCI

* [`BenchmarkCI.judge`](@ref): Run `benchmarkpkg` on `target` and `baseline`.
* [`BenchmarkCI.postjudge`](@ref): Post judgement as comment.
* [`BenchmarkCI.displayjudgement`](@ref): Print result of `BenchmarkCI.judge`.
* [`BenchmarkCI.pushresult`](@ref): Push benchmark result to a repository.
"""
module BenchmarkCI

import BenchmarkTools
import CpuId
import Dates
import GitHub
import JSON
import LinearAlgebra
import Markdown
import Tar
using Base64: base64decode
using Logging: ConsoleLogger
using Pkg: TOML
using PkgBenchmark:
    BenchmarkConfig,
    BenchmarkJudgement,
    BenchmarkResults,
    PkgBenchmark,
    baseline_result,
    export_markdown,
    target_result
using Setfield: @set, @set!
using UnPack: @unpack
using Zstd_jll: zstdmt

if VERSION < v"1.2-"
    mktempdir(f; _...) = Base.mktempdir(f)
end

include("runtimeinfo.jl")
include("gitutils.jl")

function versionof(pkg::Module)
    filepath = pathof(pkg)
    if filepath === nothing
        @error "`pathof($pkg)` does not return a path"
        return nothing
    end
    projpath = joinpath(dirname(dirname(abspath(filepath))), "Project.toml")
    if !isfile(projpath)
        @error "`$projpath` does not exist"
        return nothing
    end
    try
        proj = TOML.parsefile(projpath)
        return VersionNumber(proj["version"])
    catch err
        @error(
            "Failed to read version from `$projpath`",
            exception = (err, catch_backtrace())
        )
        return nothing
    end
end

const BENCHMARK_CI_VERSION = versionof(@__MODULE__)

macro versionof(pkg::Symbol)
    QuoteNode(versionof(getproperty(__module__, pkg)))
end

Base.@kwdef struct CIResult
    judgement::BenchmarkJudgement
    title::String = "Benchmark result"
end

const DEFAULT_WORKSPACE = ".benchmarkci"

function getnested(d, keys...)
    for k in keys
        d = get(d, k, nothing)
        d === nothing && return nothing
    end
    return d
end

function setnested!(d0, v, keys...)
    d = d0
    for k in keys[1:end-1]
        d = get!(d, k, empty(d0))
    end
    d[keys[end]] = v
    return d
end

mirrornested!(dest, src, keys...) = setnested!(dest, getnested(src, keys...), keys...)

is_in_ci(ENV = ENV) =
    lowercase(get(ENV, "CI", "false")) == "true" || haskey(ENV, "GITHUB_EVENT_PATH")

function find_project_toml(project)
    dir = isdir(project) ? project : dirname(project)
    candidates = joinpath.(dir, ("JuliaProject.toml", "Project.toml"))
    i = findfirst(isfile, candidates)
    i === nothing && return nothing
    return candidates[i]
end

function find_manifest_toml(project)
    dir = isdir(project) ? project : dirname(project)
    candidates = joinpath.(dir, ("JuliaManifest.toml", "Manifest.toml"))
    i = findfirst(isfile, candidates)
    i === nothing && return nothing
    return candidates[i]
end

function benchmarkci_info()
    return Dict(
        :versions => Dict(
            :BenchmarkCI => string(BENCHMARK_CI_VERSION),
            :PkgBenchmark => string(@versionof(PkgBenchmark)),
            :BenchmarkTools => string(@versionof(BenchmarkTools)),
        ),
        # Increment this when tweaking result format:
        :format_version => -4611686018427387904 + 3,
    )
end

function julia_info(cmd::Cmd)
    # Using `load_path_setup_code` requires that `cmd` and this
    # `julia` is compatible enough s.t. same JSON.jl version works.
    code = """
    $(Base.load_path_setup_code())
    let JSON = Base.require(Base.PkgId(
            Base.UUID("682c06a0-de6a-54ab-a142-c8b1cf79cde6"),
            "JSON",
        ))
        JSON.print(
            Dict(
                :VERSION => string(VERSION),
                :libm_name => Base.libm_name,
                :libllvm_version => string(Base.libllvm_version),
                :Sys => Dict(
                    :WORD_SIZE => Sys.WORD_SIZE,
                    :JIT => Sys.JIT,
                    :CPU_NAME => Sys.CPU_NAME,
                ),
            )
        )
    end
    """
    return open(JSON.parse, `$cmd --startup-file=no -e $code`)
end

function metadata_from(; target, baseline, pkgdir, script, project)
    return Dict(
        :target => target,
        :target_julia_info => julia_info(target.juliacmd),
        :baseline => baseline,
        :baseline_julia_info => (
            if target.juliacmd.exec[1] == baseline.juliacmd.exec[1]
                nothing
            else
                julia_info(baseline.juliacmd)
            end
        ),
        :cpu_info => cpu_info(),
        :pkgdir => pkgdir,
        :script => script,
        :project => project,
        :BenchmarkCI => benchmarkci_info(),
    )
end

function cpu_info()
    return Dict(
        :cpubrand => CpuId.cpubrand(),
        :cachesize => CpuId.cachesize(),
        :cpumodel => CpuId.cpumodel(),
        :address_size => CpuId.address_size(),
        :physical_address_size => CpuId.physical_address_size(),
        :cachelinesize => CpuId.cachelinesize(),
        :cpucores => CpuId.cpucores(),
        :cputhreads => CpuId.cputhreads(),
        :cpu_base_frequency => CpuId.cpu_base_frequency(),
        :cpu_max_frequency => CpuId.cpu_max_frequency(),
        :cpu_bus_frequency => CpuId.cpu_bus_frequency(),
        :simdbytes => CpuId.simdbytes(),
        :cpufeature => Dict(
            :TSC => CpuId.cpufeature(CpuId.TSC),
            :TSCINV => CpuId.cpufeature(CpuId.TSCINV),
            :PDCM => CpuId.cpufeature(CpuId.PDCM),
            :PDCM => CpuId.cpufeature(CpuId.PDCM),
        ),
        :perf_revision => CpuId.perf_revision(),
        :perf_fix_counters => CpuId.perf_fix_counters(),
        :perf_fix_bits => CpuId.perf_fix_bits(),
        :perf_gen_counters => CpuId.perf_gen_counters(),
        :perf_gen_bits => CpuId.perf_gen_bits(),
    )
end

function generate_script(default_script, project, should_resolve)
    default_script = abspath(default_script)
    project = abspath(project)
    """
    let Pkg = Base.require(Base.PkgId(
            Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"),
            "Pkg",
        ))
        Pkg.activate($(repr(project)))
        $(repr(should_resolve)) && Pkg.resolve()
        Pkg.instantiate()
    end
    include($(repr(default_script)))
    """
end

ensure_origin(::Nothing) = nothing
ensure_origin(config::BenchmarkConfig) = ensure_origin(config.id)
function ensure_origin(committish)
    if startswith(committish, "origin/")
        _, remote_branch = split(committish, "/"; limit = 2)
        cmd = `git fetch origin "+refs/heads/$remote_branch:refs/remotes/origin/$remote_branch"`
        @debug "Fetching $committish..." cmd
        run(cmd)
    end
end

function maybe_with_merged_project(f, project, pkgdir)
    project = abspath(project)
    if find_manifest_toml(project) !== nothing
        @info "Using existing manifest file."
        return f(project, false)  # should_resolve = false
    else
        if isfile(project)
            file = project
        else
            candidates = joinpath.(project, ("JuliaProject.toml", "Project.toml"))
            i = findfirst(isfile, candidates)
            if i === nothing
                error("One of the following files must exist:\n", join(candidates, "\n"))
            end
            file = candidates[i]
        end
        return mktempdir(prefix = "BenchmarkCI_jl_") do tmp
            tmpproject = joinpath(tmp, "Project.toml")
            cp(file, tmpproject)
            code = """
            using Pkg
            Pkg.develop(Pkg.PackageSpec(path = $(repr(pkgdir))))
            """
            run(setenv(
                `$(Base.julia_cmd()) --startup-file=no --project=$tmpproject -e $code`,
                "JULIA_LOAD_PATH" => "@:@stdlib",
            ))
            @info "Using temporary project `$tmp`."
            f(tmpproject, true)  # should_resolve = true
        end
    end
end

function format_period(seconds::Real)
    seconds < 60 && return string(floor(Int, seconds), " seconds")
    minutes = floor(Int, seconds / 60)
    return string(minutes, " minutes ", floor(Int, seconds - 60 * minutes), " seconds")
end

"""
    judge()

Run `benchmarkpkg` on `target` and `baseline`.

# Keyword Arguments
- `target :: Union{Nothing, AbstractString, BenchmarkConfig} = nothing`:
  Benchmark target configuration.  Default to the checked out working tree.  A
  git commitish can be passed as a string.  `PkgBenchmark.BenchmarkConfig`
  can be used for more detailed control (see PkgBenchmark.jl manual).
- `baseline :: Union{AbstractString, BenchmarkConfig} = "origin/master"`:
  Benchmark baseline configuration.  See `target.`
- `pkgdir :: AbstractString = pwd()`: Package root directory.
- `script :: AbstractString = "\$pkgdir/benchmark/benchmarks.jl"`: The script
  that defines the `SUITE` global variable/constant.
- `project :: AbstractString = dirname(script)`: The project used to define
  and run benchmarks.
- `postprocess`, `retune`, `verbose`, `logger_factory`: Passed to
  `PkgBenchmark.benchmarkpkg`.
"""
judge(; target = nothing, baseline = "origin/master", kwargs...) =
    judge(target, baseline; kwargs...)

function judge(
    target,
    baseline = "origin/master";
    workspace = DEFAULT_WORKSPACE,
    pkgdir = pwd(),
    script = joinpath(pkgdir, "benchmark", "benchmarks.jl"),
    project = dirname(script),
    logger_factory = is_in_ci() ? ConsoleLogger : nothing,
    kwargs...,
)
    target = BenchmarkConfig(target)
    if !(baseline isa BenchmarkConfig)
        baseline = @set target.id = baseline
    end

    # Make sure `origin/master` etc. exists:
    ensure_origin(target)
    ensure_origin(baseline)

    mkpath(workspace)
    script_wrapper = abspath(joinpath(workspace, "benchmarks_wrapper.jl"))

    let metadata = metadata_from(;
            target = target,
            baseline = baseline,
            pkgdir = pkgdir,
            script = script,
            project = project,
        )
        open(joinpath(workspace, "metadata.json"); write = true) do io
            JSON.print(io, metadata)
        end
    end
    let env = Base.active_project(),
        projecttoml = find_project_toml(env),
        manifesttoml = find_manifest_toml(env),
        dest = joinpath(workspace, "runnerenv")
        mkpath(dest)
        cp(projecttoml, joinpath(dest, "Project.toml"); force = true)
        cp(manifesttoml , joinpath(dest, "Manifest.toml"); force = true)
    end

    maybe_with_merged_project(project, pkgdir) do tmpproject, should_resolve
        tmpprojecttoml = find_project_toml(tmpproject)
        cp(tmpprojecttoml, joinpath(workspace, "Project.toml"); force = true)
        tmpmanifest = something(find_manifest_toml(tmpproject))
        cp(tmpmanifest , joinpath(workspace, "Manifest.toml"); force = true)
        write(script_wrapper, generate_script(script, tmpproject, should_resolve))
        _judge(;
            target = target,
            baseline = baseline,
            workspace = workspace,
            pkgdir = pkgdir,
            benchmarkdir = dirname(script),
            benchmarkpkg_kwargs = (;
                kwargs...,
                logger_factory = logger_factory,
                script = script_wrapper,
            ),
        )
    end
end

function noisily(f, yes::Bool = is_in_ci(); interval = 60 * 5)
    if yes
        t0 = time_ns()
        timer = Timer(interval; interval = interval) do _
            tstr = format_period(floor(Int, (time_ns() - t0) / 1e9))
            @info "$tstr passed.  Still running `judge`..."
        end
        try
            f()
        finally
            close(timer)
        end
    else
        f()
    end
end

function _judge(; target, baseline, workspace, pkgdir, benchmarkdir, benchmarkpkg_kwargs)

    target_git_info = git_commit_info(something(target.id, "HEAD"), benchmarkdir)
    baseline_git_info = git_commit_info(something(baseline.id, "HEAD"), benchmarkdir)
    noisily() do
        start_date = Dates.now()
        time_target = @elapsed group_target = PkgBenchmark.benchmarkpkg(
            pkgdir,
            target;
            resultfile = joinpath(workspace, "result-target.json"),
            benchmarkpkg_kwargs...,
        )
        mid_date = Dates.now()
        @debug("`git status`", output = Text(read(`git status`, String)))
        @debug("`git diff`", output = Text(read(`git diff`, String)))
        time_baseline = @elapsed group_baseline = PkgBenchmark.benchmarkpkg(
            pkgdir,
            baseline;
            resultfile = joinpath(workspace, "result-baseline.json"),
            benchmarkpkg_kwargs...,
        )
        end_date = Dates.now()
        @info """
        Finish running benchmarks.
        * Target: $(format_period(time_target))
        * Baseline: $(format_period(time_baseline))
        """
        let run_info = Dict(
                :target => merge!(
                    Dict(
                        :start_date => start_date,
                        :end_date => mid_date,
                        :elapsed_time => time_target,
                    ),
                    target_git_info,
                ),
                :baseline => merge!(
                    Dict(
                        :start_date => mid_date,
                        :end_date => end_date,
                        :elapsed_time => time_baseline,
                    ),
                    baseline_git_info,
                ),
            )
            open(joinpath(workspace, "run_info.json"); write = true) do io
                JSON.print(io, run_info)
            end
        end
        judgement = PkgBenchmark.judge(group_target, group_baseline)
        if is_in_ci()
            display(judgement)
        end
        return judgement
    end
end

function git_commit_info(rev::AbstractString, benchmarkdir::AbstractString)
    git = `git --no-pager`
    benchmarkdir = relpath(benchmarkdir)
    git_tree_sha1_benchmark = try
        strip(read(`$git rev-parse $rev:$benchmarkdir`, String))
    catch err
        nothing
    end
    git_tree_sha1_src = try
        strip(read(`$git rev-parse $rev:src`, String))
    catch err
        nothing
    end
    info = Dict(
        :git_tree_sha1 => strip(read(`$git rev-parse $rev:`, String)),
        :git_tree_sha1_benchmark => git_tree_sha1_benchmark,
        :git_tree_sha1_src => git_tree_sha1_src,
        :git_commit_sha1 => strip(read(`$git rev-parse $rev`, String)),
        :git_commit_date =>
            strip(read(`$git show --no-patch --format=%cI $rev`, String)),
    )
end

function _loadjudge(workspace)
    group_target = PkgBenchmark.readresults(joinpath(workspace, "result-target.json"))
    group_baseline = PkgBenchmark.readresults(joinpath(workspace, "result-baseline.json"))
    return PkgBenchmark.judge(group_target, group_baseline)
end

# Used only for testing:
_loadciresult(workspace::AbstractString = DEFAULT_WORKSPACE) =
    CIResult(judgement = _loadjudge(workspace))

"""
    postjudge(; title = "Benchmark result")

Post judgement as a comment in the corresponding pull request.
"""
postjudge(workspace::AbstractString = DEFAULT_WORKSPACE; kwargs...) =
    postjudge(_loadjudge(workspace); kwargs...)

function postjudge(judgement::BenchmarkJudgement; title = "Benchmark result")
    event_path = get(ENV, "GITHUB_EVENT_PATH", nothing)
    if event_path !== nothing
        post_judge_github(event_path, CIResult(judgement = judgement, title = title))
        return
    end
    displayjudgement(judgement)
end

function printcommentmd(io, ciresult)
    println(io, "")
    println(io, "", ciresult.title, "")
    println(io)
    printresultmd(io, ciresult)
    println(io)
    println(io)
    println(io, "")
end

function printresultmd(io, ciresult)
    judgement = ciresult.judgement
    println(io, "# Judge result")
    export_markdown(io, judgement)
    println(io)
    println(io)
    println(io, "---")
    println(io, "# Target result")
    export_markdown(io, target_result(judgement))
    println(io)
    println(io)
    println(io, "---")
    println(io, "# Baseline result")
    export_markdown(io, baseline_result(judgement))
    println(io)
    println(io)
    println(io, "---")
    println(io, "# Runtime information")
    show(io, MIME"text/markdown"(), runtimeinfo())
    println(io)
    md = try
        sprint() do buffer
            show(buffer, MIME"text/markdown"(), CpuId.cpuinfo())
        end
    catch err
        @error(
            """`show(_, "text/markdown", CpuId.cpuinfo())` failed""",
            exception = (err, catch_backtrace())
        )
        nothing
    end
    md === nothing || println(io, md)
end

function printcommentjson(io, ciresult)
    comment = sprint() do io
        printcommentmd(io, ciresult)
    end
    # https://developer.github.com/v3/issues/comments/#create-a-comment
    JSON.print(io, Dict("body" => comment::AbstractString))
end

function error_on_missing_github_token()
    error("""
    Environment variable `GITHUB_TOKEN` is not set.  The workflow file
    must contain configuration such as:

        - name: Post result
          run: julia -e "using BenchmarkCI; BenchmarkCI.postjudge()"
          env:
            GITHUB_TOKEN: \${{ secrets.GITHUB_TOKEN }}
    """)
end

function post_judge_github(event_path, ciresult)
    event = JSON.parsefile(event_path)
    url = event["pull_request"]["comments_url"]
    # https://developer.github.com/v3/activity/events/types/#pullrequestevent
    @debug "Posting to: $url"

    GITHUB_TOKEN = get(ENV, "GITHUB_TOKEN", nothing)
    GITHUB_TOKEN === nothing && error_on_missing_github_token()

    cmd = ```
    curl
    --include
    --request POST
    $url
    -H "Content-Type: application/json"
    -H "Authorization: token $GITHUB_TOKEN"
    --data @-
    ```

    response = sprint() do stdout
        open(pipeline(cmd, stdout = stdout, stderr = stderr), write = true) do io
            printcommentjson(io, ciresult)
        end
    end
    @debug "Response from GitHub" Text(response)
    @info "Comment posted."
end

"""
    pushresult(; url, branch, sshkey, title)

Push benchmark result to `branch` in `url`.

# Keyword Arguments
- `url::Union{AbstractString,Nothing} = nothing`: Repository URL.
- `branch::AbstractString = "benchmark-results"`: Branch where the
  results are pushed.
- `sshkey::Union{AbstractString,Nothing} = nothing`: Documenter.jl-style
  SSH private key (base64-encoded private key).  Default to `ENV["SSH_KEY"]`.
- `title::AbstractString = "Benchmark result"`: The title to be used in
  benchmark report.
"""
function pushresult(;
    url::Union{AbstractString,Nothing} = nothing,
    branch::AbstractString = "benchmark-results",
    sshkey::Union{AbstractString,Nothing} = nothing,
    workspace::AbstractString = DEFAULT_WORKSPACE,
    title::AbstractString = "Benchmark result",
)
    workspace = abspath(workspace)
    default_url = nothing
    repo = nothing
    commit_info = github_commit_info()
    if haskey(ENV, "GITHUB_TOKEN")
        sha = github_sha()
        auth = GitHub.authenticate(ENV["GITHUB_TOKEN"])
        repo = GitHub.repo(ENV["GITHUB_REPOSITORY"]; auth = auth)
        default_url = "git@github.com:$(repo.full_name).git"
    end
    if url === nothing
        default_url === nothing && error("cannot auto detect url")
        url = default_url
    end
    if sshkey === nothing
        let key = get(ENV, "SSH_KEY", "")
            if key == ""
                @info "Empty `SSH_KEY`. Skipping `pushresult`."
                displayjudgement(workspace)
                return
            else
                sshkey = String(base64decode(key))
            end
        end
    end

    judgement = _loadjudge(workspace)
    if commit_info !== nothing
        open(joinpath(workspace, "commit_info.json"); write = true) do io
            JSON.print(io, commit_info)
        end
    end
    local datadir
    GitUtils.updating(
        url,
        branch;
        sshkey = sshkey,
        commit_message = "Add: $title",
    ) do ctx
        datadir = Dates.format(Dates.now(), joinpath("yyyy", "mm", "dd", "HHMMSS"))
        mkpath(datadir)
        compress_tar(joinpath(datadir, "result.tar.zst"), workspace)
        open(joinpath(datadir, "result.md"); write = true) do io
            println(io, "# ", title)
            println(io)
            printinfomd(io, commit_info)
            println(io)
            printresultmd(io, CIResult(title = title, judgement = judgement))
        end
    end
    if repo !== nothing
        status_params = Dict(
            "state" => "success",
            "context" => "benchmarkci/$title",
            "description" => "Benchmarks complete!",
            "target_url" => as_target_url(url, branch, datadir),
        )
        @info "Creating status" sha params = status_params
        GitHub.create_status(repo, sha; auth = auth, params = status_params)
    end
    return
end

# From `post_status` in Documenter.jl
function github_sha()
    if get(ENV, "GITHUB_EVENT_NAME", nothing) == "pull_request"
        event_path = get(ENV, "GITHUB_EVENT_PATH", nothing)
        event_path === nothing && return
        event = JSON.parsefile(event_path)
        if haskey(event, "pull_request") &&
            haskey(event["pull_request"], "head") &&
            haskey(event["pull_request"]["head"], "sha")
            return event["pull_request"]["head"]["sha"]
        end
    elseif get(ENV, "GITHUB_EVENT_NAME", nothing) == "push"
        return ENV["GITHUB_SHA"]
    end
end

function as_https(url::AbstractString)
    if (m = match(r"^git@github.com:([^/]+)/(.*?)(\.git)?$", url)) !== nothing
        # handle `git@github.com:$USER/$REPO.git`
        return "https://github.com/$(m[1])/$(m[2])"
    end
    return nothing
end

function as_target_url(url::AbstractString, branch::AbstractString, datadir::AbstractString)
    https = as_https(url)
    https === nothing && return nothing
    return "$https/blob/$branch/$datadir/result.md"
end

struct GitHubCommitInfo
    repo_url::Union{Nothing,String}
    pull_request::Dict{String,Any}
end

GitHubCommitInfo() = GitHubCommitInfo(nothing, Dict{String,Any}())

function github_commit_info()
    commit_info::GitHubCommitInfo = GitHubCommitInfo()

    event_path = get(ENV, "GITHUB_EVENT_PATH", nothing)
    event_path === nothing && return commit_info
    event = JSON.parsefile(event_path)

    # https://docs.github.com/en/developers/webhooks-and-events/webhook-events-and-payloads#pull_request
    if get(ENV, "GITHUB_EVENT_NAME", nothing) == "pull_request"
        if (pull_request = get(event, "pull_request", nothing)) !== nothing
            for keys in Tuple[
                ("html_url",),
                ("title",),
                ("user", "login"),
                ("head", "sha"),
                ("head", "ref"),
                ("head", "repo", "html_url"),
                ("base", "sha"),
                ("base", "ref"),
            ]
                mirrornested!(commit_info.pull_request, pull_request, keys...)
            end
        end
    end

    @set! commit_info.repo_url = getnested(event, "repository", "html_url")

    return commit_info
end

function printinfomd(io, commit_info::GitHubCommitInfo)
    @unpack repo_url, pull_request = commit_info
    pull_request_url = getnested(pull_request, "html_url")
    pull_request_title = getnested(pull_request, "title")
    pull_request_commit_sha = getnested(pull_request, "head", "sha")
    target_commit_url = nothing

    if repo_url !== nothing && pull_request_commit_sha !== nothing
        target_commit_url = "$repo_url/commit/$pull_request_commit_sha"
    end

    if target_commit_url !== nothing
        println(
            io,
            "* Pull request commit:",
            " [`$pull_request_commit_sha`]($target_commit_url)",
        )
    end
    if pull_request_url !== nothing
        print(io, "* Pull request: <$pull_request_url>")
        if pull_request_title !== nothing
            print(io, " ($pull_request_title)")
        end
        println(io)
    end
end

function compress_tar(dest, src)
    proc = zstdmt() do cmd
        open(`$cmd -f - -o $dest`; write = true)
    end
    Tar.create(src, proc)
    close(proc)
    wait(proc)
end

function decompress_tar(dest, src)
    proc = zstdmt() do cmd
        open(pipeline(`$cmd -d`, stdin = src); read = true)
    end
    Tar.extract(proc, dest)
    close(proc)
    wait(proc)
end

"""
    displayjudgement()

Print result of `BenchmarkCI.judge`.
"""
displayjudgement(workspace::AbstractString = DEFAULT_WORKSPACE) =
    displayjudgement(_loadjudge(workspace))

function displayjudgement(judgement::BenchmarkJudgement)
    io = IOBuffer()
    printresultmd(io, CIResult(judgement = judgement))
    seekstart(io)
    display(Markdown.parse(io))
    display(Text("\n"))
end

runall(args...; kwargs...) = postjudge(judge(args...; kwargs...))

end # module

ext: jl


text:
module ThermalMemoryEraser
using Reexport

include("force/Forces.jl")
@reexport using .Force

include("update_methods.jl")
@reexport using .UpdateMethods

include("Trajectories.jl")
@reexport using .Trajectories

include("metrics.jl")
@reexport using .Metrics

include("configuration.jl")
@reexport using .Configurations

include("ensemble.jl")
@reexport using .Ensembles

end # module


text:
import Term: cleantext, chars, textlen, split_lines, replace_multi

cleanstring(x) = replace(string(x), "\r\n" => "\n")
cleansprint(fn, x) = replace(sprint(fn, x), "\r\n" => "\n")

"""
Write one or multiple strings to a text file so that they can
later be used for testing.
"""
tofile(content::Vector, filepath) =
    open(filepath, "w") do io
        map(ln -> println(io, ln), content)
        close(io)
    end

tofile(content::String, filepath) =
    open(filepath, "w") do io
        print(io, content)
        close(io)
    end

"""
Load a string representation of a renderable from file
and clean it up
"""
fromfile(filepath) = replace_multi(read(filepath, String), "\\n" => "\n", "\\e" => "\e")

fromfilelines(filepath) = lines = readlines(filepath)

same_widths(text::String) = length(unique(textlen.(split_lines(text)))) == 1

function check_widths(text, width)
    for line in split_lines(text)
        @test textlen(line) <= width
    end
end

"""
Extensively test a panel making sure it has
the right size and Measure
"""
function testpanel(p, w, h)
    # isnothing(h) || println(vLine(h) * p)
    # check all lines have the same length
    _p = string(p)

    dw = displaysize(stdout)[2]
    if isnothing(w) || w > dw
        return nothing
    else
        widths = textwidth.(cleantext.(split(_p, "\n")))
    end

    # println(p, p.measure, widths)
    @test length(unique(widths)) == 1

    # check it has the right measure
    if !isnothing(w)
        @test p.measure.w == w
        @test textlen(cleantext(p.segments[1].text)) == w
        @test length(chars(cleantext(p.segments[1].text))) == w
    end

    if !isnothing(h)
        @test p.measure.h == h
        @test length(p.segments) == h
    end
end

function testtree(p, w, h)
    # check all lines have the same length
    _p = string(p)

    dw = displaysize(stdout)[2]
    if isnothing(w) || w > dw
        return nothing
    else
        widths = textwidth.(cleantext.(split(_p, "\n")))
    end

    # println(p, p.measure, widths)
    @test length(unique(widths)) == 1

    # check it has the right measure
    if !isnothing(w)
        @test p.measure.w == w
    end

    if !isnothing(h)
        @test p.measure.h == h
        @test length(p.segments) == h
    end
end

nlines(x) = length(split(x, "\n"))
lw(x) = max(length.(split(x, "\n"))...)

"""
Include but with a timer
"""

macro timeit_include(path::AbstractString)
    return :(@timeit TIMEROUTPUT $path include($path))
end

# ------------------------------- test no throw ------------------------------ #

struct NoException <: Exception end

macro test_nothrow(ex)
    return esc(:(@test_throws NoException ($(ex); throw(NoException()))))
end


text:
__precompile__(false)

eye(T,n,m) = Matrix{T}(LinearAlgebra.I,n,m)

Mat2(T) = Mat2x2(T) = zeros(T,2,2)
Mat3(T) = Mat3x3(T) = zeros(T,3,4)
Mat4(T) = Mat4x4(T) = zeros(T,4,4)

const Mat2f = const Mat2x2f = SMatrix{2,2,Float32}
const Mat3f = const Mat3x3f = SMatrix{3,3,Float32}
const Mat4f = const Mat4x4f = SMatrix{4,4,Float32}

const zerosMat4x4f = Mat4(Float32)
const eyeMat4x4f = eye(Float32,4,4)

export eye
export Mat2
export Mat3
export Mat4
export Mat2x2
export Mat2x2f
export Mat3x3
export Mat3x3f
export Mat4x4
export Mat4x4f
export zerosMat4x4f
export eyeMat4x4f

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
# As the submarine starts booting up things like the Retro Encabulator, you
# realize that maybe you don't need all these submarine features after all.
#
# What is the smallest model number accepted by MONAD?

# read input
program = map(eachline("input.txt")) do l
  (opcode, operands...) = split(l, ' ')
  op1 = operands[1][1]
  op2 = nothing
  if length(operands) > 1
    op2 = operands[2]
    if '0' <= op2[1] <= '9' || op2[1] == '-'
      op2 = parse(Int, op2)
    else
      op2 = op2[1]
    end
  end
  (opcode, op1, op2)
end

# location of all inp instructions
inpip = accumulate((i,_) -> findnext(((o,_,_),) -> o == "inp", program, i+1), 1:14)

# look at part 1 for an explanation
result = Vector{Char}(undef, 14)
stack = Vector{NTuple{2,Int}}()
sizehint!(stack, 7)
for i in 1:14
  ip = inpip[i]
  (bbb,ccc) = (program[ip+5][3], program[ip+15][3])
  if bbb > 0
    push!(stack, (i, ccc))
  else
    (j,delta) = pop!(stack)
    delta += bbb
    result[i] = max(1, 1 + delta) + '0'
    result[j] = max(1, 1 - delta) + '0'
  end
end

println("Model: ", String(result))


text:
import CacheServers: iscached, update!, pull, clear!

struct Scalar{T}
    val::T
end

struct Grad{T}
    val::T
end

struct Param{T}
    val::T
end

mutable struct Variable{T}
    val::T
    grad::T
end

Variable(val::T) where T = Variable(val, T(0))

iscached(val) = true

function update!(var::Variable{T}, grad::Grad{T}) where {T <: Number}
    var.grad = grad.val
    var
end

function update!(var::Variable{T}, param::Param{T}) where {T <: Number}
    var.val = param.val
    var
end

pull(var::Variable, ::Type{Grad}) = var.grad
pull(var::Variable, ::Type{Param}) = var.val
clear!(var::Variable{<:Number}) = var.grad = 0


text:
function show{FSA <: FixedVector}(io::IO, a::Vector{FSA})
	println(io, length(a), "-element", typeof(a), ":")
	for elem in a
		print(io, " ")
		showcompact(io, elem)
		println(io)
	end
end

immutable MaxFunctor <: Functor{2} end
immutable MinFunctor <: Functor{2} end
immutable ExtremaFun <: Functor{2} end

@compat (::MaxFunctor)(a, b) = max(a, b)
@compat (::MinFunctor)(a, b) = min(a, b)
@compat (::ExtremaFun)(reducevalue, a) = min(reducevalue[1], a), max(reducevalue[2], a)

minimum{T <: FixedArray}(a::Vector{T}) = reduce(MinFunctor(), a)
maximum{T <: FixedArray}(a::Vector{T}) = reduce(MaxFunctor(), a)
extrema{T <: FixedArray}(a::AbstractVector{T}) = reduce(ExtremaFun(), a)

function isapprox{FSA <: FixedArray, A <: FixedArray}(x::FSA, y::A; rtol::Real=Base.rtoldefault(eltype(x),eltype(y)), atol::Real=0, norm::Function=vecnorm)
    # Same behaviour as AbstractArray: julia/base/linalg/generic.jl
	d = norm(x - y)
    return isfinite(d) ? d <= atol + rtol*max(norm(x), norm(y)) : x == y
end
function isapprox{FSA <: FixedArray, A <: Array}(x::FSA, y::A; rtol::Real=Base.rtoldefault(eltype(x),eltype(y)), atol::Real=0, norm::Function=vecnorm)
    d = norm(Array(x) - y) # operations between Arrays and FixedArrays not defined...
    return isfinite(d) ? d <= atol + rtol*max(norm(x), norm(y)) : x == y
end

(.+){T<:FixedArray, ND}(a::Array{T, ND}, x::T) = T[elem .+ x for elem in a]
(+){T<:FixedArray, ND}(a::Array{T, ND}, x::T) = a .+ x


text:
module AlgebraHelpersTest

using Test
using ReactiveMP
using Random

using LinearAlgebra

@testset "Helpers" begin

    @testset "negate_inplace!" begin
        rng = MersenneTwister(1234)

        for size in 2:4
            A = rand(rng, size, size)
            B = copy(A)
            @test B == A
            @test ReactiveMP.negate_inplace!(B) == -A
            @test B == -A

            a = rand(rng)
            b = copy(a)
            @test ReactiveMP.negate_inplace!(b) == -a
            @test b == a # numbers stay unchanged
        end
    end

    @testset "negate_inplace!" begin
        rng = MersenneTwister(1234)

        for size in 2:4
            A = rand(rng, size, size)
            B = copy(A)
            c = randn(rng)
            @test B == A
            @test ReactiveMP.mul_inplace!(c, B) == c * A
            @test B == c * A

            a = rand(rng)
            b = copy(a)
            c = randn(rng)
            @test ReactiveMP.mul_inplace!(c, b) == c * a
            @test b == a # numbers stay unchanged
        end
    end

    @testset "rank1update" begin
        rng = MersenneTwister(1234)

        for size in 2:4, T1 in (Float32, Float64), T2 in (Float32, Float64), T3 in (Float32, Float64)
            A = rand(rng, T1, size, size)
            x = rand(rng, T2, size)
            @test ReactiveMP.rank1update(A, x) == (A + x*x')
            y = rand(rng, T3, size)
            @test ReactiveMP.rank1update(A, x, y) == (A + x*y')
        end

    end

    @testset "mul_trace" begin
        rng = MersenneTwister(1234)

        for size in 2:4, T1 in (Float32, Float64), T2 in (Float32, Float64)
            A = rand(rng, T1, size, size)
            B = rand(rng, T2, size, size)
            @test ReactiveMP.mul_trace(A, B) ≈ tr(A * B)

            a = rand(rng, T1)
            b = rand(rng, T2)
            @test ReactiveMP.mul_trace(a, b) ≈ a * b
        end

    end

end

end


text:
using Libdl: dlext

# Platform-independent library name
libcptr_path() = joinpath(".", "libcptr."*dlext)

# Set array length to 4
n = Cint(4)

# Allocate object in Fortran
cptr = ccall((:allocatestuff, libcptr_path()), Ptr{Cvoid}, (Ref{Cint},), n)

# Print some information about it
ccall((:enquirestuff, libcptr_path()), Cvoid, (Ref{Ptr{Cvoid}},), cptr)

# Release memory
ccall((:clearstuff, libcptr_path()), Cvoid, (Ref{Ptr{Cvoid}},), cptr)

# Print pointer value
println(" ", Int(cptr))

# Pass string from Julia to Fortran
ccall((:print_string, libcptr_path()), Cvoid, (Cstring,), "wololo😊")

# Pass array from Julia to Fortran
values = Cint[1, 1, 2, 3, 5, 8, 13]
ccall((:printintarray, libcptr_path()), Cvoid, (Ref{Cint}, Ref{Cint},), values, length(values))

max_stars_repo_name: sloede/calling-fortran-from-julia

max_issues_repo_name: sloede/calling-fortran-from-julia

max_forks_repo_name: sloede/calling-fortran-from-julia


text:
module test

immutable Test
    sex::ASCIIString
    age::Int64
    function Student(s::ASCIIString,a::Int64)
        @assert a>=18
        new(s,a)
    end
    function get_sex(st::Student)
        st.sex
    end
end

@show st = Student("FM",19)
@show st.sex = "M"
@show get_sex(st)

end


text:
module QueryableBackend

import IteratorInterfaceExtensions, TableTraits, QueryOperators

include("queryable/queryable.jl")
include("queryable/queryable_filter.jl")
include("queryable/queryable_map.jl")

include("source_queryable.jl")

end


text:
#############################################################################
#############################################################################
#
# A script that runs all unit tests in the project.
#
#############################################################################
#############################################################################

include("../poly_factorization_project.jl")

####
# Execute unit tests for integers
###
include("integers_test.jl")
test_euclid_ints()
test_ext_euclid_ints()

####
# Execute unit tests for polynomials
####
include("polynomials_test.jl")
prod_test_poly()
prod_derivative_test_poly()
ext_euclid_test_poly()
division_test_poly()

####
# Execute unit tests for polynomial factorization
####
include("factorization_test.jl")
factor_test_poly()


text:
"""
Returns an array of data filtered according to absorbing state being the index of the array.
"""
function absorb_filter_data(data, maph::MAPHDist)
    p, q = model_size(maph)
    filter_data = []

    for i = 1:q
        temp_data = filter(data) do obs first(obs).a == i end
        push!(filter_data, temp_data)
    end

    #QQQQ this whole function can be just this one line ??
    # [filter((obs)->first(obs).a == i, data) for i in 1:q]
    return filter_data
end

"""
Returns an array of data filtered according to the sojourn time being less than the index of the time array.
"""
function time_filter_data(data, n::Int64)

    times = ((d)->first(d).y).(data)
    min_time = quantile(times, 0.1)
    max_time = quantile(times, 0.9)
    @show min_time, max_time
    time_vec = Array(LinRange(min_time, max_time, n)) #QQQQ - cleanup later

    filter_data = []

    for i = 1:(n-1)
        temp_data = filter(data) do obs first(obs).y ≥ time_vec[i] && first(obs).y < time_vec[i+1] end
        push!(filter_data,temp_data)
        @show length(temp_data)
    end

    return filter_data, time_vec
end

"""
Test the sufficient stats with full trace data
"""
function full_trace_sufficient_stats_test(;N=10^4)
    Λ₄, λ45, λ54, Λ₅ = 5, 2, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]

    maph = MAPHDist([0.5,0.5]',T_example, T0_example)

    p,q = model_size(maph)
    test_stats = MAPHSufficientStats(maph)

    for _ in 1:N
        times, states = rand(maph, full_trace = true)
        ss = sufficient_stat_from_trajectory(maph, times, states)
        test_stats.N += ss.N
        test_stats.Z += ss.Z
        test_stats.B += ss.B
    end

    test_stats.B = test_stats.B/sum(test_stats.B)
    test_stats.Z = test_stats.Z/sum(test_stats.Z)
    test_stats.N = test_stats.N/sum(test_stats.N)

    test_α = test_stats.B

    @show test_stats.Z

    computered_intensity = (test_stats.N.*(1 ./ test_stats.Z))

    @show computered_intensity

    # @show(stats.N,test_stats.N)
    return true
end

"""
QQQQ
"""
function sufficient_stats_test(; sim_runs = 10^6)

    #Set heuristically the number of times steps to be the sqrt of the number of runs
    timesteps = round(Int,sqrt(sim_runs))

    Λ₄, λ45, λ54, Λ₅ = 5, 2, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]
    initial_dist = [0.5,0.5]

    maph = MAPHDist(initial_dist', T_example, T0_example)

    p,q = model_size(maph)

    data = []
    full_trace =[]

    @showprogress "Simulating data" for i in 1:sim_runs
        times, states = rand(maph, full_trace = true)
        push!(full_trace, (times,states))
        push!(data, (observation_from_full_traj(times,states),i))
    end

    absorb = absorb_filter_data(data, maph)
    absorb_data = vcat(absorb[1],absorb[2],absorb[3])

    #QQQQ this is for now just looking at trajectory observed in first state
    time_bin, time_vec = time_filter_data(absorb_data, timesteps)

    Z_errors = Float64[]
    N_errors = Float64[]

    #loop over all bins
    # println("\n start initialization...")
    # @showprogress "Checking sufficient stats"

    for i in 1:length(time_bin)
        ss_i = MAPHSufficientStats[]
        for trace in full_trace[last.(time_bin[i])]
            ss = sufficient_stat_from_trajectory(maph, trace[1], trace[2])
            push!(ss_i, ss)
        end

        if !isempty(ss_i)
            mean_observed_ss = mean(ss_i)
            time_slice = first(data[last(last.(time_bin[i]))])[1]
            obs = first(data[last(last.(time_bin[i]))])
            @show obs
            computed_ss = sufficient_stats(obs, maph)
            # @show computed_ss.N, mean_observed_ss.N

            errs_N = (mean_observed_ss.N - computed_ss.N)#./ computed_ss
            push!(N_errors, norm(mean_observed_ss.N - computed_ss.N)/time_slice)

            # computed_ss.B ≈ initial_dist || return false

            # @show computed_ss.Z, mean_observed_ss.Z
            push!(Z_errors, norm(mean_observed_ss.Z - computed_ss.Z)/time_slice)

            # @show obs
            # sufficient_stats()
        end
    end

    return N_errors, Z_errors, time_vec

    # ab1 = filter(x->x.a==1,first.(data))
    # ab2 = filter(x->x.a==2,first.(data))
    # ab3 = filter(x->x.a==3,first.(data))

    # probs = [length(ab1),length(ab2),length(ab3)]./length(data)
    # means = [mean(first.(ab1)),mean(first.(ab2)),mean(first.(ab3))]
    # scvs = [scv(first.(ab1)),scv(first.(ab2)),scv(first.(ab3))]

    # dist = MAPHDist(10,probs,means,scvs)
    # println("finish initialization")

    # println("starting simulations")
    # fit!(dist,first.(data))
    # println("\nfinished simulations")

    return true

end

"""
QQQQ - An example with a deterministic path....
"""
function sufficient_stats_test3()

    ϵ = 0.00000001

    T_example = [-(1.0+4ϵ) 1.0-4ϵ ϵ
                ϵ -(1.0 + 4ϵ) 1.0-4ϵ
                ϵ ϵ -(1 + 4ϵ) ]

    T0_example = [ϵ ϵ ϵ;
                  ϵ ϵ ϵ;
                  1-4ϵ ϵ ϵ]

    T_example2 = [-(1.0+3ϵ) 1.0 ϵ
                 ϵ -(1.0+3ϵ) 1.0
                 ϵ ϵ -(1+3ϵ)]
    T0_example2 = [ϵ ϵ;
                   ϵ ϵ;
                   1-3ϵ ϵ]

    display(sum(T0_example2[:,1]))
    maph = MAPHDist([0.5,0.5, 0.0]',T_example, T0_example)

    obs = (y=100.0, a=2)
    sufficient_stats(obs, maph)
end

"""
QQQQ - this is probably temporary.
"""
function analyze_ss_with_plots()

    Λ₄, λ45, λ54, Λ₅ = 5, 2, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]
    initial_dist = [0.5,0.5]

    maph = MAPHDist(initial_dist', T_example, T0_example)

    SingleObs = NamedTuple{(:y, :a), Tuple{Float64, Int64}}

    for y in 0.5:0.5
        obs = (y=y, a = 1)
        ss = sufficient_stats(obs,maph)
        ss.N[1,1] #this means the expected number of jumps from the first transient state
                        #to the first absorbing state (a=1)
        # @show ss.N
        return ss.N
    end

end


text:
@testset "MJCore" begin

    @testset "mjmujoco" begin
        include("mjmujoco.jl")
    end

    @testset "cglobals" begin
        include("cglobals.jl")
    end

    @testset "Uninitialized constructor" begin
        skipped = (mjModel, mjData)
        toobig = (mjUI, mjuiDef, mjvFigure)
        uninit_types = Iterators.filter(T->!(T in skipped), MuJoCo.union_types(MJCore.MJTypes))
        @testset "$T" for T in uninit_types
            if T in toobig
                @test_skip test_uninit(T)
            else
                @test test_uninit(T)
            end
        end
    end
end


text:
function parse_bnf(bnf::IO, syntax)

	# using a flexible grammar incorporating common variants of BNF and EBNF so as to avoid need for user to specify the variant

	# TODO for antlr4:
	# handling escape sequences in literals and some in regular expressions, e.g. // and /'
	# regular expressions which are not bracketed expressions
	# multi-line comments
	# directives (-->skip)
	if !(syntax in [:ebnf, :antlr4])
		error("specify :ebnf or :antlr4 as the syntax")
	end

	# whitespace and comments
	whitespace = p"\s"
	if syntax == :antlr4
		comment = P"//.*$" | P"/\*.*\*/"
	else
		comment = P"\(\*.*\*\)"
	end
	ws = Drop(Star(whitespace | comment))

	# variables
	identifier = p"\w+"
	# allow variable names with and without <...> delimiter:
	variablename = identifier | (E"<" + identifier + E">")
	# allow ::=, :, and = as symbol separating lhs and rhs of each rule
	definesymbol = E"::=" | E":" | E"="
	variabledef = variablename + ws + definesymbol								|> ns->ASTNode(:variabledef, ASTNode[], Dict{Symbol,Any}(:name=>ns[1]))
	variableref = variablename													|> ns->ASTNode(:variableref, ASTNode[], Dict{Symbol,Any}(:name=>ns[1]))

	# terminals
	# allow both "..." and '...' as delimiters for literals
	terminal = (E"'" + p"[^']*" + E"'") | (E"\"" + p"[^\"]*" + E"\"") 			|> ns->ASTNode(:literal, ASTNode[], Dict{Symbol,Any}(:value=>ns[1]))

	# regular expression (using [] syntax meaning optionality elsewhere) - we retain square brackets
	if syntax == :antlr4
		# todo - proper handling of escaped ] etc.
		regexp = p"\[.+\]"														|> ns->ASTNode(:regexp, ASTNode[], Dict{Symbol,Any}(:value=>ns[1]))
	end

	# symbols
	if syntax == :antlr4
		symbol = terminal |  regexp | variableref
	else
		symbol = terminal |  variableref
	end

	# bracketed expressions
	# (...) grouping
	# [...] optional
	# {...} repeat
	expr = Delayed()
	groupexpr = E"(" + ws + expr + ws + E")"
	if syntax != :antlr4
		# for anltr, square brackets indiciate regular expressions, so omit this
		optionalexpr = E"[" + ws + expr + ws + E"]"								|> ns->ASTNode(:optional, convert(Array{ASTNode}, ns))
	end
	repeatexpr = E"{" + ws + expr + ws + E"}"									|> ns->ASTNode(:quantifier, convert(Array{ASTNode}, ns), Dict{Symbol,Any}(:min=>0))

	# operators
	# precedence (highest first)
	# ~ negation (Antlr)
	# * (repetition)
	# , (concatenate)
	# | (choice)
	if syntax == :antlr4
		negatableexpr = symbol | groupexpr		# TODO negation don't make sense for optional and repeat exprs? - in any case, these aren't used by ANTLR
		negexpr = E"~" + ws + negatableexpr										|> ns->ASTNode(:negation, convert(Array{ASTNode}, ns))
	end

	if syntax == :antlr4
		quantifiableexpr = symbol | groupexpr | negexpr		# TODO quantifiers don't make sense for optional and repeat exprs?
	else
		quantifiableexpr = symbol | groupexpr	# TODO quantifiers don't make sense for optional and repeat exprs?
	end
	quesexpr = quantifiableexpr + ws + E"?"										|> ns->ASTNode(:optional, convert(Array{ASTNode}, ns))
	starexpr = quantifiableexpr + ws + E"*"										|> ns->ASTNode(:quantifier, convert(Array{ASTNode}, ns), Dict{Symbol,Any}(:min=>0))
	plusexpr = quantifiableexpr + ws + E"+"										|> ns->ASTNode(:quantifier, convert(Array{ASTNode}, ns), Dict{Symbol,Any}(:min=>1))

	if syntax == :antlr4
		primaryexpr = symbol | groupexpr | repeatexpr | starexpr | plusexpr | quesexpr | negexpr
	else
		primaryexpr = symbol | groupexpr | optionalexpr | repeatexpr | starexpr | plusexpr | quesexpr
	end

	concatsymbol = E"," | P"\s+"
	concatexpr = Star(primaryexpr + ws + concatsymbol + ws) + primaryexpr		|> ns->ASTNode(:concatenation, convert(Array{ASTNode}, ns))

	altsymbol = E"|"
	altexpr = Star(concatexpr + ws + altsymbol + ws) + concatexpr				|> ns->ASTNode(:alternation, convert(Array{ASTNode}, ns))

	expr.matcher = altexpr

	# rules
	# allow ;, ., or end-of-line to terminate rule
	ruleterminatorsym = E";" | E"." | P"$"
	rule = P"^" + ws + variabledef + ws + expr + ws + ruleterminatorsym			|> ns->begin ns[1].children=ASTNode[ns[2],]; ns[1] end

	# grammar
	grammar = ws + Star(rule + ws) + rule + ws 									|> ns->ASTNode(:bnf, convert(Array{ASTNode}, ns))

	astarray = parse_try(bnf, Try(grammar + Eos()))
 	# note: Try not only enables (limited?) backtracking but also
 	# permits use of parse_try to process whole file (parse_one cannot be used
 	# in the same way in this particular idiom)

	if isempty(astarray)
		error("No production rules found")
	end

	astarray[1]

end


text:
using Revise
using Traceur
using BenchmarkTools
using ZonalFlow

domain  = Domain(extent=(2π,2π),res=(4,4));
coeffs  = Coefficients(Ω=2π,θ=0.0,μ=0.01,ν=1.0,ν₄=1.0,linear=false);
forcing = Stochastic(kf=3,dk=1,ε=0.01);
prob    = BetaPlane(domain,coeffs,forcing);

tspan   = (0.0,1000.0);
tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            save_start=false
           );

# eq = NL()
# eq = CE2()
# eq = GQL(2)
eq = GCE2(2)

u0 = get_de_ic(prob,eq);
p  = get_de_params(prob,eq)

du = similar(u0);
@benchmark f!($du,$u0,$p,$tspan)

_prob,_alg = get_de_probalg(prob,eq,u0,tspan,p);
@benchmark solve($_prob,$_alg;$tsargs...)


text:
function steadystate!(m::HetDSGE;
                      βlo::Float64 = 0.5*exp(m[:γ])/(1 + m[:r]),
                      βhi::Float64 = exp(m[:γ])/(1 + m[:r]),
                      excess::Float64 = 5000.,
                      tol::Float64 = 1e-4,
                      maxit::Int64 = 20,
                      βband::Float64 = 1e-2)

    # Load settings
    nx = get_setting(m, :nx)
    ns = get_setting(m, :ns)

    xlo = get_setting(m, :xlo)
    xhi = get_setting(m, :xhi)
    zlo = get_setting(m, :zlo)
    zhi = get_setting(m, :zhi)

    # Load parameters
    R = 1 + m[:r].value
    H = m[:H].value
    η = m[:η].value
    γ = m[:γ].value
    ω = m[:ωstar].value
    T = m[:Tstar].value

    # Load grids
    xgrid = m.grids[:xgrid].points
    sgrid = m.grids[:sgrid].points
    xswts = m.grids[:weights_total]
    f     = m.grids[:fgrid]

    βlo_temp = βlo
    βhi_temp = βhi

    counter = 1
    n  = ns*nx
    c  = zeros(n)
    bp = zeros(n)
    KF = zeros(n,n)
    μ  = zeros(n)

    # Initial guess
    β   = 1.0
    Win = Vector{Float64}(undef, n)
    Win_guess = ones(n)

    if get_setting(m, :use_last_βstar) && !isnan(m[:βstar].value)
        βlo = βhi = m[:βstar].value
    elseif !isnan(m[:βstar].value)
        βlo_temp = m[:βstar].value - βband
        βhi_temp = m[:βstar].value + βband

        c, bp, Win, KF = policy_hetdsge(nx, ns, βlo_temp, R, ω, H, η, T, γ,
                                        zhi, zlo, xgrid, sgrid, xswts, Win_guess, f)
        excess_lo, μ = compute_excess(xswts, KF, bp)

        if excess_lo < 0 && abs(excess_lo) > tol
            βlo = βlo_temp
            c, bp, Win, KF = policy_hetdsge(nx, ns, βhi_temp, R, ω, H, η, T, γ, zhi,
                                            zlo, xgrid, sgrid, xswts, Win_guess, f)
            excess_hi, μ = compute_excess(xswts, KF, bp)

            if excess_hi > 0
                βhi = βhi_temp
            end
        elseif excess_lo >= 0
            βhi = βlo_temp
        end
    end

    while abs(excess) > tol && counter < maxit # clearing markets
        β = (βlo + βhi) / 2.0
        c, bp, Win, KF = policy_hetdsge(nx, ns, β, R, ω, H, η, T, γ, zhi, zlo, xgrid, sgrid,
                                        xswts, Win_guess, f)
        excess, μ = compute_excess(xswts, KF, bp)

        # bisection
        if excess > 0
            βhi = β
        elseif excess < 0
            βlo = β
        end
        counter += 1
        if get_setting(m, :use_last_βstar) && !isnan(m[:βstar].value)
            break
        end
    end
    m[:lstar] = Win
    m[:cstar] = c
    m[:μstar] = μ
    m[:βstar] = β
    m[:β_save] = β

    #return Win, c, μ, β
    nothing
end

@inline function compute_excess(xswts::Vector{T}, KF::Matrix{T}, bp::Vector{T}) where T<:Float64
    LPMKF = xswts[1] * KF
    # Find eigenvalue closest to 1
    (D,V) = (eigen(LPMKF)...,)
    max_D = argmax(abs.(D))
    D = D[max_D]

    if abs(D - 1) > 2e-1 # that's the tolerance we are allowing
        @warn "Your eigenvalue is too far from 1, something is wrong."
    end
    # Pick eigenvector associated w/ largest eigenvalue and moving it back to values
    μ = real(V[:,max_D])
    μ = μ ./ dot(xswts, μ) # Scale of eigenvectors not determinate: rescale to integrate to 1
    excess = dot(xswts, (μ .* bp)) # compute excess supply of savings, which is a fn of w

    return excess, μ
end

function policy_hetdsge(nx::Int, ns::Int, β::AbstractFloat, R::AbstractFloat,
                        ω::AbstractFloat, H::AbstractFloat,
                        η::AbstractFloat, T::AbstractFloat,
                        γ::AbstractFloat, zhi::AbstractFloat,
                        zlo::AbstractFloat,
                        xgrid::Vector{Float64}, sgrid::Vector{Float64},
                        xswts::Vector{Float64}, Win::Vector{Float64},
                        f::Array{Float64,2}, damp::Float64 = 0.5, dist::Float64 = 1.,
                        tol::Float64 = 1e-4, maxit::Int64 = 500)
    n    = nx*ns
    c    = zeros(n)                # consumption
    bp   = Vector{Float64}(undef, n)      # savings
    Wout = Vector{Float64}(undef, length(Win))
    counter = 1
    qfunction_hetdsge(x::Float64) = mollifier_hetdsge(x, zhi, zlo) #/sumz

    while dist>tol && counter<maxit # for debugging
        # compute c(w) given guess for Win = β*R*E[u'(c_{t+1})]
        for iss in 1:ns
            for ia in 1:nx
                c[nx*(iss-1)+ia] = min(1/Win[nx*(iss-1)+ia],xgrid[ia]+η)
            end
        end
        bp = repeat(xgrid, ns) - c  # compute bp(w) given guess for Win
        Wout = parameterized_expectations_hetdsge(nx, ns, β, R, ω, H, T, γ,
                                                  qfunction_hetdsge, xgrid, sgrid, xswts, c, bp, f)

        dist = maximum(abs.(Wout - Win))
        Win  = damp*Wout + (1.0-damp) * Win
        counter += 1
    end
    if counter == maxit
        @warn "Euler iteration did not converge"
    end
    tr = kolmogorov_fwd_hetdsge(nx, ns, ω, H, T, R, γ, qfunction_hetdsge, xgrid, sgrid, bp, f)
    return c, bp, Wout, tr
end

@inline function parameterized_expectations_hetdsge(nx::Int, ns::Int, β::AbstractFloat,
                                                    R::AbstractFloat, ω::AbstractFloat,
                                                    H::AbstractFloat, T::AbstractFloat,
                                                    γ::AbstractFloat,
                                                    qfunc::Function,
                                                    xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                                    xswts::Vector{Float64}, c::Vector{Float64},
                                                    bp::Vector{Float64}, f::Array{Float64,2})
    l_out = zeros(nx*ns)
    for iss=1:ns
        for ia=1:nx
            sumn = 0.0
            for isp=1:ns
                for iap=1:nx
                    sumn += (xswts[nx*(isp-1)+iap]/c[nx*(isp-1)+iap]) *
                        qfunc((xgrid[iap] - R*(exp(-γ))*bp[nx*(iss-1)+ia] - T) /
                              (ω*H*sgrid[isp])) * f[iss,isp] ./ sgrid[isp]
                end
            end
            l_out[nx*(iss-1)+ia] = (β*R*(exp(-γ))/ω*H)*sumn
        end
    end
    return l_out
end

@inline function kolmogorov_fwd_hetdsge(nx::Int, ns::Int, ω::AbstractFloat,
                                H::AbstractFloat, T::AbstractFloat,
                                R::AbstractFloat, γ::AbstractFloat,
                                qfunc::Function,
                                xgrid::Vector{Float64},
                                sgrid::Vector{Float64}, bp::Vector{Float64}, f::Array{Float64,2})
    tr = zeros(nx*ns,nx*ns)
    for iss=1:ns
        for ia=1:nx
            for isp=1:ns
                for iap=1:nx
                    tr[nx*(isp-1)+iap, nx*(iss-1)+ia] = qfunc((xgrid[iap] - R*(exp(-γ))*bp[nx*(iss-1)+ia] - T)/(ω*H*sgrid[isp])) * f[iss,isp] ./ (ω * H * sgrid[isp])
                end
            end
        end
    end
    return tr
end

@inline function mollifier_hetdsge(z::AbstractFloat,ehi::AbstractFloat,elo::AbstractFloat)
    # mollifier function
    In = 0.443993816237631
    if z<ehi && z>elo
        temp = -1.0 + 2.0 * (z - elo) / (ehi - elo)
        return (2.0 / (ehi - elo)) * exp(-1.0 / (1.0 - temp^2)) / In
    end
    return 0.0
end

@inline function dmollifier_hetdsge(x::AbstractFloat, ehi::AbstractFloat, elo::AbstractFloat)
    In = 0.443993816237631
    if x<ehi && x>elo
        temp = (-1.0 + 2.0*(x-elo)/(ehi-elo))
        out  = -(2*temp./((1 - temp.^2).^2)).*(2/(ehi-elo)).*mollifier_hetdsge(x, ehi, elo)
    else
        out = 0.0
    end
    return out
end


text:
# --- Models --- #

"""
    AbstractModel

Supertype for all models.
"""
abstract type AbstractModel end

"""
    NoStateModel <: AbstractModel

Supertype for all models which contain no state variables.
"""
abstract type NoStateModel <: AbstractModel end

"""
    InterfaceModel <: AbstractModel

Supertype for all models which are used to extend 2D models to 3D models.
"""
abstract type InterfaceModel <: AbstractModel end

"""
    ResidualModel <: AbstractModel

Supertype for all models whose state variables may be defined by solving a set
of residual equations.
"""
abstract type ResidualModel <: AbstractModel end

"""
    UnsteadyModel <: AbstractModel

Supertype for all models which represent dynamic systems of equations.
"""
abstract type UnsteadyModel <: AbstractModel end

# --- Inplaceness Trait --- #

abstract type InPlaceness end
struct InPlace <: InPlaceness end
struct OutOfPlace <: InPlaceness end

"""
    inplaceness(::Type{T})

Return `InPlace()` if functions associated with model `T` are in-place
or `OutOfPlace()` if functions associated with model `T` are out-of-place.
"""
inplaceness(::Type{T}) where T

# default definition for models with no state variables
inplaceness(::Type{T}) where T<:NoStateModel = OutOfPlace()

# default definition for coupled models
function inplaceness(::Type{T}) where T <: NTuple{N,AbstractModel} where N
    model_types = (T.parameters...,)
    if isinplace(coupling_inplaceness(model_types...)) || any(isinplace.(model_types))
        return InPlace()
    else
        return OutOfPlace()
    end
end

"""
    coupling_inplaceness(::Type{T1}, ::Type{T2}, ..., ::Type{TN})

Return `InPlace()` if the functions associated with the coupling function
for coupled models `T1`, `T2`, ... `TN` are in-place or `OutOfPlace()`
if the functions associated with the coupling function for coupled models `T1`,
`T2`, ... `TN` are out-of-place.
"""
coupling_inplaceness(::Vararg{Type,N}) where N

isinplace(model::T) where T = isinplace(inplaceness(T))
isinplace(::Type{T}) where T = isinplace(inplaceness(T))
isinplace(::InPlace) = true
isinplace(::OutOfPlace) = false

# --- Matrix Type Trait --- #

abstract type MatrixType end
struct Empty <: MatrixType end
struct Zeros <: MatrixType end
struct Identity <: MatrixType end
struct Constant <: MatrixType end
struct Invariant <: MatrixType end
struct Linear <: MatrixType end
struct Nonlinear <: MatrixType end

"""
    rate_jacobian_type(::Type{T})

Return
 - `Empty()`, if the jacobian of the residual expression associated with model `T`
    with respect to the state rates is an empty matrix
 - `Zeros()`, if the jacobian of the residual expression associated with model `T`
    with respect to the state rates is a zero matrix
 - `Identity()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state rates is the identity matrix
 - `Invariant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state rates is independent of the state rates,
    state variables, inputs, parameters, and time.
 - `Constant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state rates is independent of the state rates,
    state variables, inputs, and time.
 - `Linear()`, if the jacobian of the residual expression associated with model `T`
    with respect to the state rates may vary with respect to time, and is
    linear with respect to the state rates.
 - `Nonlinear()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state rates may vary with respect to time,
    and is nonlinear with respect to the state rates.

If no method is defined for the specified type, return `Nonlinear()`.
"""
rate_jacobian_type(::Type{T}) where T = Nonlinear()

# models with no state variables have no rate jacobian
rate_jacobian_type(::Type{T}) where T <: NoStateModel = Empty()

# definition for combinations of models
function rate_jacobian_type(::Type{T}) where T<:NTuple{N,AbstractModel} where N
    model_types = (T.parameters...,)
    if isempty(coupling_rate_jacobian_type(model_types...)) &&
        all(isempty.(rate_jacobian_type.(model_types)))
        return Empty()
    elseif iszero(coupling_rate_jacobian_type(model_types...)) &&
        all(iszero.(rate_jacobian_type.(model_types)))
        return Zeros()
    elseif iszero(coupling_rate_jacobian_type(model_types...)) &&
        all(isidentity.(rate_jacobian_type.(model_types)))
        return Identity()
    elseif isinvariant(coupling_rate_jacobian_type(model_types...)) &&
        all(isinvariant.(input_jacobian_type.(model_types))) &&
        all(isinvariant.(rate_jacobian_type.(model_types)))
        return Invariant()
    elseif isconstant(coupling_rate_jacobian_type(model_types...)) &&
        all(isconstant.(input_jacobian_type.(model_types))) &&
        all(isconstant.(rate_jacobian_type.(model_types)))
        return Constant()
    else
        return Linear()
    end
end

"""
    state_jacobian_type(::Type{T})

Return
 - `Empty()`, if the jacobian of the residual expression associated with model `T`
    with respect to the state variables is an empty matrix
 - `Zeros()`, if the jacobian of the residual expression associated with model `T`
    with respect to the state variables is a zero matrix
 - `Identity()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state variables is the identity matrix
 - `Invariant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state variables is independent of the state rates,
    state variables, inputs, parameters, and time.
 - `Constant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state variables is independent of the state rates,
    state variables, inputs, and time.
 - `Linear()`, if the jacobian of the residual expression associated with model `T`
    with respect to the state variables may vary with respect to time, and is
    linear with respect to the states.
 - `Nonlinear()`, if the jacobian of the residual expression associated with
    model `T` with respect to the state variables may vary with respect to time,
    and is nonlinear with respect to the states.

If no method is defined for the specified type, return `Nonlinear()`.
"""
state_jacobian_type(::Type{T}) where T = Nonlinear()

# models with no state variables have no state jacobian
state_jacobian_type(::Type{T}) where T <: NoStateModel = Empty()

# definition for combinations of models
function state_jacobian_type(::Type{T}) where T<:NTuple{N,AbstractModel} where N
    model_types = (T.parameters...,)
    if isempty(coupling_state_jacobian_type(model_types...)) &&
        all(isempty.(state_jacobian_type.(model_types)))
        return Empty()
    elseif iszero(coupling_state_jacobian_type(model_types...)) &&
        all(iszero.(state_jacobian_type.(model_types)))
        return Zeros()
    elseif iszero(coupling_state_jacobian_type(model_types...)) &&
        all(isidentity.(state_jacobian_type.(model_types)))
        return Identity()
    elseif isinvariant(coupling_state_jacobian_type(model_types...)) &&
        all(isinvariant.(input_jacobian_type.(model_types))) &&
        all(isinvariant.(state_jacobian_type.(model_types)))
        return Constant()
    elseif isconstant(coupling_state_jacobian_type(model_types...)) &&
        all(isconstant.(input_jacobian_type.(model_types))) &&
        all(isconstant.(state_jacobian_type.(model_types)))
        return Constant()
    elseif islinear(coupling_state_jacobian_type(model_types...)) &&
        all(islinear.(input_jacobian_type.(model_types))) &&
        all(islinear.(state_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

"""
    input_jacobian_type(::Type{T})

Return
 - `Empty()`, if the jacobian of the residual expression associated with model `T`
    with respect to the inputs is an empty matrix
 - `Zeros()`, if the jacobian of the residual expression associated with model `T`
    with respect to the inputs is a zero matrix
 - `Identity()`, if the jacobian of the residual expression associated with
    model `T` with respect to the inputs is the identity matrix
 - `Invariant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the inputs is independent of the state rates,
    state variables, inputs, parameters, and time.
 - `Constant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the inputs is independent of the state rates,
    state variables, inputs, and time.
 - `Linear()`, if the jacobian of the residual expression associated with model `T`
    with respect to the inputs may vary with respect to time, and is
    linear with respect to the inputs.
 - `Nonlinear()`, if the jacobian of the residual expression associated with
    model `T` with respect to the inputs may vary with respect to time,
    and is nonlinear with respect to the inputs.

If no method is defined for the specified type, return `Nonlinear()`.
"""
input_jacobian_type(::Type{T}) where T = Nonlinear()

# models with no state variables have no input jacobian
input_jacobian_type(::Type{T}) where T<:NoStateModel = Empty()

# definition for combinations of models
function input_jacobian_type(::Type{T}) where T<:NTuple{N,AbstractModel} where N
    model_types = (T.parameters...,)
    if all(isempty.(input_jacobian_type.(model_types)))
        return Empty()
    elseif all(iszero.(input_jacobian_type.(model_types)))
        return Zeros()
    elseif all(isidentity.(input_jacobian_type.(model_types)))
        return Identity()
    elseif all(isinvariant.(input_jacobian_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(input_jacobian_type.(model_types)))
        return Constant()
    elseif all(islinear.(input_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

"""
    parameter_jacobian_type(::Type{T})

Return
 - `Empty()`, if the jacobian of the residual expression associated with model `T`
    with respect to the parameters is an empty matrix
 - `Zeros()`, if the jacobian of the residual expression associated with model `T`
    with respect to the parameters is a zero matrix
 - `Identity()`, if the jacobian of the residual expression associated with
    model `T` with respect to the parameters is the identity matrix
 - `Invariant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the parameters is independent of the state rates,
    state variables, inputs, parameters, and time.
 - `Constant()`, if the jacobian of the residual expression associated with
    model `T` with respect to the parameters is independent of the state rates,
    state variables, inputs, and time.
 - `Linear()`, if the jacobian of the residual expression associated with model `T`
    with respect to the parameters may vary with respect to time, and is
    linear with respect to the parameters.
 - `Nonlinear()`, if the jacobian of the residual expression associated with
    model `T` with respect to the parameters may vary with respect to time,
    and is nonlinear with respect to the parameters.

If no method is defined for the specified type, return `Nonlinear()`.
"""
parameter_jacobian_type(::Type{T}) where T = Nonlinear()

# models with no state variables have no input jacobian
parameter_jacobian_type(::Type{T}) where T <: NoStateModel = Empty()

# definition for combinations of models
function parameter_jacobian_type(::Type{T}) where T<:NTuple{N,AbstractModel} where N
    model_types = (T.parameters...,)
    if all(isempty.(parameter_jacobian_type.(model_types)))
        return Empty()
    elseif all(iszero.(parameter_jacobian_type.(model_types)))
        return Zeros()
    elseif all(isidentity.(parameter_jacobian_type.(model_types)))
        return Identity()
    elseif all(isinvariant.(parameter_jacobian_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(parameter_jacobian_type.(model_types)))
        return Constant()
    elseif all(islinear.(parameter_jacobian_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

"""
    time_gradient_type(::Type{T})

Return
 - `Empty()`, if the derivative of the residual expression associated with model `T`
    with respect to time is an empty matrix
 - `Zeros()`, if the derivative of the residual expression associated with model `T`
    with respect to time is a zero matrix
 - `Invariant()`, if the derivative of the residual expression associated with
    model `T` with respect to time is independent of the state rates,
    state variables, inputs, parameters, and time.
 - `Constant()`, if the derivative of the residual expression associated with
    model `T` with respect to time is independent of the state rates,
    state variables, inputs, and time.
 - `Linear()`, if the derivative of the residual expression associated with model `T`
    with respect to time may vary with respect to time, and is
    linear with respect to time.
 - `Nonlinear()`, if the derivative of the residual expression associated with
    model `T` with respect to time may vary with respect to time,
    and is nonlinear with respect to time.

If no method is defined for the specified type, return `Nonlinear()`.
"""
time_gradient_type(::Type{T}) where T = Nonlinear()

# models with no state variables have no input jacobian
time_gradient_type(::Type{T}) where T <: NoStateModel = Empty()

# definition for combinations of models
function time_gradient_type(::Type{T}) where T<:NTuple{N,AbstractModel} where N
    model_types = (T.parameters...,)
    if all(isempty.(time_gradient_type.(model_types)))
        return Empty()
    elseif all(iszero.(time_gradient_type.(model_types)))
        return Zeros()
    elseif all(isinvariant.(time_gradient_type.(model_types)))
        return Invariant()
    elseif all(isconstant.(time_gradient_type.(model_types)))
        return Constant()
    elseif all(islinear.(time_gradient_type.(model_types)))
        return Linear()
    else
        return Nonlinear()
    end
end

coupling_rate_jacobian_type(model_types...) = Nonlinear()
coupling_state_jacobian_type(model_types...) = Nonlinear()
coupling_parameter_jacobian_type(model_types...) = Nonlinear()
coupling_time_gradient_type(model_types...) = Nonlinear()

isempty(::MatrixType) = false
isempty(::Empty) = true

iszero(::MatrixType) = false
iszero(::Zeros) = true

isidentity(::MatrixType) = false
isidentity(::Identity) = true

isinvariant(::MatrixType) = false
isinvariant(::Empty) = true
isinvariant(::Zeros) = true
isinvariant(::Identity) = true

isconstant(::MatrixType) = false
isconstant(::Empty) = true
isconstant(::Zeros) = true
isconstant(::Identity) = true
isconstant(::Invariant) = true
isconstant(::Constant) = true

islinear(::MatrixType) = false
islinear(::Empty) = true
islinear(::Zeros) = true
islinear(::Identity) = true
islinear(::Invariant) = true
islinear(::Constant) = true
islinear(::Linear) = true


text:
module Dynamics
    using OrbitalTrajectories.SpiceUtils

    using SPICE
    using Plots
    using Unitful
    using DifferentialEquations
    using DiffEqBase
    using DiffEqSensitivity
    using SciMLBase
    using LinearAlgebra
    using ModelingToolkit
    using NLsolve
    using LineSearches
    using ForwardDiff
    using FiniteDiff
    using Memoize
    using SimpleTraits
    using StaticArrays
    using ProgressMeter
    using RecipesBase
    using PhysicalConstants.CODATA2014: NewtonianConstantOfGravitation

    # Required to ensure that we can precompile ODEFunctions (which needs to be
    # done in our own cache).
    using RuntimeGeneratedFunctions
    RuntimeGeneratedFunctions.init(@__MODULE__)

    #----------------#
    # ABSTRACT TYPES #
    #----------------#

    # Dynamical models
    abstract type Abstract_DynamicalModel end
    abstract type Abstract_ModelODEFunctions <: ModelingToolkit.AbstractODESystem end
    abstract type Abstract_R3BPModel <: Abstract_DynamicalModel end
    abstract type Abstract_R4BPModel <: Abstract_DynamicalModel end

    # Reference frames
    abstract type Abstract_ReferenceFrame end

    # Differential correctors
    abstract type Abstract_DifferentialCorrector end

    #----------#
    # INCLUDES #
    #----------#
    include("orbital_trajectories.jl")
    include("reference_frames.jl")
    include("sensitivities.jl")
    include("dynamical_models.jl")
    include("differential_correction.jl")
    include("plotting.jl")
end


text:
using OBC, MPM
using HDF5, JLD
using PyPlot
using Distributions

kmax = 1
D = 2
iters = 5000

loc = "/home/bana/largeresearch/rnaseq.jld"
if !isfile(loc)
    loc = "/home/bana/largeresearch/seq-data/rnaseq.jld"
end
@load loc luad_norm luad_norm_sampids luad_norm_geneids

#pick two random genes
data = luad_norm[[4,5],:]'
#data = luad_norm[[2,3],:]'

#println(luad_norm_geneids[[2,3]])
#plot(data[:,1], data[:,2], "go", label="brca")
#assert(false)

d = 100.0
######################################################################
# Class 0
######################################################################
cov = eye(D,D)
#cov = [0.5 -0.2; -0.2 1.0]
mu = ones(D,kmax)
#mu = [1.6 -1.0; 1.6 -1.0]
#mu = [0.2 0.2]'
prior = MPM.MPMPrior(D=2, kmax=kmax, kappa=10., S=eye(2)*10)
start = MPM.MPMParams(mu, #mu :: Matrix{Float64}
    cov, #sigma :: Matrix{Float64}
    ones(kmax)/kmax, #w :: Vector{Float64}
    clamp(log(data'/d),-3.0,Inf), #lam :: Matrix{Float64}
    1) #k :: Int

pmoves = MPM.MPMPropMoves()
obj_a = MPM.MPMSampler(prior, data, deepcopy(start), pmoves, d)
obj_a.usepriors = true

opt = Opt(:LN_COBYLA, 2)
lower_bounds!(opt, [-Inf, 0.])
xtol_rel!(opt,1e-4)

min_objective!(opt, myfunc)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)

(minf,minx,ret) = optimize(opt, [1.234, 5.678])
println("got $minf at $minx after $count iterations (returned $ret)")


text:
import Revise
import Agtor

using Profile, BenchmarkTools, OwnTime, Logging, Infiltrator

using CSV
using Dates
using DataStructures
using DataFrames
using Agtor

import Flatten

# Start with the environment variable set for multi-threading
# $ JULIA_NUM_THREADS=4 ./julia

# Windows:
# $ set JULIA_NUM_THREADS=4 && julia dev.jl

function setup_zone(data_dir::String="test/data/")
    zone_spec_dir::String = "$(data_dir)zones/"
    zone_specs::Dict{String, Dict} = load_yaml(zone_spec_dir)

    zone_params = generate_agparams("", zone_specs["Zone_1"])

    collated_specs::Array = []
    agparams = collect_agparams!(zone_params, collated_specs; ignore=["crop_spec"])

    climate_fn::String = "$(data_dir)climate/farm_climate_data.csv"
    climate::Climate = load_climate(climate_fn)

    return create(zone_params, climate), agparams
end

data_dir = "test/data/"

z1, collated = setup_zone(data_dir)

samples = DataFrame!(CSV.File("test/data/scenarios/sampled_params.csv"))

@info names(samples)

@info "before"
@info z1.fields[1].irrigation.capital_cost

update_params!(z1, samples[1, :])

@info "after"
@info z1.fields[1].irrigation.capital_cost

update_params!(z1, samples[2, :])

@info z1.fields[1].irrigation.capital_cost

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
@testset "convert" begin

bf = BigFloat(pi)
df = Double64(bf)
sqrt2 = sqrt(2.0)

@test BigFloat(df) == convert(BigFloat, df)
@test Double64(bf) == convert(Double64, bf)

@test Double32(sqrt2) == convert(Double32, sqrt2)
@test Double64(sqrt2) == convert(Double64, sqrt2)

end # convert

@testset "maxintfloat $T" for T in (Double32, Double64)
    @test isinteger(maxintfloat(T))

    # Previous integer is representable, next integer is not
    @test maxintfloat(T) == (maxintfloat(T) - one(T)) + one(T)
    @test maxintfloat(T) != (maxintfloat(T) + one(T)) - one(T)
end

@testset "predicates $T" for T in (Double32, Double64)

    @test iszero(zero(T))
    @test isone(one(T))
    @test ispositive(one(T))
    @test isnegative(-one(T))
    @test isnonpositive(-one(T))
    @test isnonpositive(zero(T))
    @test isnonnegative(one(T))
    @test isnonnegative(zero(T))

    @test isinf(T(Inf))
    @test isposinf(T(Inf))
    @test isneginf(T(-Inf))
    @test isnan(T(NaN))

    @test isnormal(one(T))
    @test isinteger(one(T))
    @test isfractional(one(T)/2)

end # predicates

@testset "double $T" for T in (Double32, Double64)

    @test HI(one(T)) == one(T)
    @test LO(one(T)) == zero(T)
    @test HILO(one(T)) == (one(T), zero(T))

    sqrt2  = sqrt(T(2))
    hi, lo = sqrt2.hi, sqrt2.lo
    hilo   = HILO(sqrt2)

    @test HI(sqrt2) == hi
    @test LO(sqrt2) == lo
    @test HILO(sqrt2) == (hi, lo)

    @test HI(hilo) == hi
    @test LO(hilo) == lo
    @test HILO(hilo) == hilo

end # Double

@testset "compare $T" for T in (Double32, Double64)

    two = 2
    fp2 = T(2)
    sqrt2  = sqrt(T(2))
    sqrt5  = sqrt(T(5))

    @test   sqrt2 == sqrt2
    @test   sqrt2 >= sqrt2
    @test   sqrt2 <= sqrt2
    @test !(sqrt2 != sqrt2)
    @test !(sqrt2 >  sqrt2)
    @test !(sqrt2 <  sqrt2)

    @test !(sqrt2 == sqrt5)
    @test !(sqrt2 >= sqrt5)
    @test   sqrt2 <= sqrt5
    @test   sqrt2 != sqrt5
    @test !(sqrt2 >  sqrt5)
    @test   sqrt2 <  sqrt5

    @test !(fp2 == sqrt5)
    @test !(fp2 >= sqrt5)
    @test   fp2 <= sqrt5
    @test   fp2 != sqrt5
    @test !(fp2 >  sqrt5)
    @test   fp2 <  sqrt5

    @test !(sqrt5 == fp2)
    @test !(sqrt5 <= fp2)
    @test   sqrt5 >= fp2
    @test   sqrt5 != fp2
    @test !(sqrt5 <  fp2)
    @test   sqrt5 >  fp2

    @test !(two == sqrt5)
    @test !(two >= sqrt5)
    @test   two <= sqrt5
    @test   two != sqrt5
    @test !(two >  sqrt5)
    @test   two <  sqrt5

    @test !(sqrt5 == two)
    @test !(sqrt5 <= two)
    @test   sqrt5 >= two
    @test   sqrt5 != two
    @test !(sqrt5 <  two)
    @test   sqrt5 >  two

end # compare


text:
function expr_arr_to_block(exprs)
  block = :(begin end)
  foreach(expr -> push!(block.args, expr), exprs)
  block
end


text:
function s_t_connectivity(nodes::Vector{Int}, source::Vector{Int}, target::Vector{Int})

    φ = function (system::Array{Float64,2}, x::Vector)
        A = copy(system)
        s = copy(source)
        t = copy(target)

        failed = setdiff(nodes, x)

        A[failed, :] .= 0
        A[:, failed] .= 0

        while true
            if length(intersect(s, t)) > 0
                return true
            end

            A[:, s] .= 0
            s = findall(x -> x > 0, vec(sum(A[s, :], dims = 1)))

            if length(s) == 0
                return false
            end
        end

    end

    return φ
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module DefineSave

export TrajLog
export trajSave, trajLoad, getSaveFileRoot
export readdirSavs, isSave

using RLESUtils, FileUtils, Loggers

import Base: getindex, copy!, haskey

immutable TrajLog
    log::TaggedDFLogger
end

function trajSave(fileroot::AbstractString, d::TrajLog)
    file = LogFile("$fileroot")
    save_log(file, d.log)
    name(file) #return
end

getindex(d::TrajLog, x) = getindex(d.log, x)

function trajLoad(logfile::AbstractString)
    file = LogFile(logfile)
    TrajLog(load_log(file)) #return
end

readdirSavs(dir::AbstractString=".") = readdir_ext(".sav", dir)
isSave(f::AbstractString) = endswith(f, ".sav")

function getSaveFileRoot(file::AbstractString)
    splitext(file)[1]
end

function copy!(dst::TrajLog, src::TrajLog)
    for k in keys(src.log)
        dst.log[k] = src.log[k]
    end
    dst
end

haskey(d::TrajLog, k::String) = haskey(d.log, k)

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
module TestVerbosity

using ParallelKMeans
using Random
using Test
using Suppressor

@testset "LLoyd: Testing verbosity of implementation" begin
    Random.seed!(2020)
    X = rand(4, 150)
    Random.seed!(2020)
    # Capture output and compare
    r = @capture_out kmeans(Lloyd(), X, 3; n_threads=1, max_iters=1, verbose=true)
    @test r == "Iteration 1: Jclust = 46.534795844478815\n"
end

@testset "Hamerly: Testing verbosity of implementation" begin
    Random.seed!(2020)
    X = rand(4, 150)
    Random.seed!(2020)
    # Capture output and compare
    r = @capture_out kmeans(Hamerly(), X, 3; n_threads=1, max_iters=1, verbose=true)
    @test r == "Iteration 1: Jclust = 46.534795844478815\n"
end

@testset "Elkan: Testing verbosity of implementation" begin
    Random.seed!(2020)
    X = rand(4, 150)
    Random.seed!(2020)
    # Capture output and compare
    r = @capture_out kmeans(Elkan(), X, 3; n_threads=1, max_iters=1, verbose=true)
    @test r == "Iteration 1: Jclust = 46.534795844478815\n"
end

end # module


text:
module test_massComputation

using  Modia3D
import Modia3D.ModiaMath
using Modia3D.StaticArrays
using Modia3D.Unitful

vmat1 = Modia3D.Material(color="Red", transparency=0.5)
vmat2 = Modia3D.Material(color="Blue", transparency=0.5)

rCM1 = zeros(3)
rCM1[2] = 9.0
rCM2 = zeros(3)
rCM2[1] = 5.0
rCM2[3] = 2.0
I1 = fill(2.0,(3,3))
I2 = fill(3.0,(3,3))
massProp1 = Modia3D.MassProperties(15.0, rCM1, I1)
massProp2 = Modia3D.MassProperties(5.0, rCM2, I2)

rCM3 = zeros(3)
rCM3[2] = 4.0
rCM4 = zeros(3)
rCM4[1] = 3.0
rCM4[3] = 9.0
I3 = fill(45.0,(3,3))
I4 = fill(100.0,(3,3))
massProp3 = Modia3D.MassProperties(30.0, rCM3, I3)
massProp4 = Modia3D.MassProperties(60.0, rCM4, I4)

@signal Sine(;A=1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Sine, sim::ModiaMath.SimulationState)
    signal.y.value  = signal.A*sin(2*pi*signal.freqHz*sim.time)
end

@signal LinearSig(;A=1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::LinearSig, sim::ModiaMath.SimulationState)
    signal.y.value  = signal.A*sim.time - signal.freqHz
end

@assembly DoublePendulum(;length=1.0, width=0.2) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))

   part1 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), massProp1, vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   part2 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), massProp2, vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])

   part3 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), massProp3, vmat2),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   part4 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), massProp4, vmat2),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])

   rev1 = Modia3D.Revolute(world          , part1.frames[1])
   Modia3D.connect(part1.frames[2], part2.frames[1]; R=ModiaMath.rot2(45u"°"))

   rev2 = Modia3D.Revolute(part2.frames[2], part3.frame0)
   Modia3D.connect(part3.frames[2], part4.frames[1] ; R=ModiaMath.rot2(-45u"°"))

   sig1 = Sine(A=1.5,freqHz=0.5)
   sig2 = LinearSig(A=2.5,freqHz=1.0)
   sine1     = Modia3D.SignalToFlangeAngle(sig1.y)
   sine2     = Modia3D.SignalToFlangeAngle(sig2.y)
   Modia3D.connect(sine1,   rev1)
   Modia3D.connect(sine2,   rev2)

end

doublePendulum = DoublePendulum(sceneOptions = Modia3D.SceneOptions(visualizeFrames=true,
                                                                    defaultFrameLength=0.3,
                                                                    enableContactDetection=false))
model = Modia3D.SimulationModel( doublePendulum, analysis=ModiaMath.KinematicAnalysis )
result = ModiaMath.simulate!(model, stopTime=3.0)

ModiaMath.plot(result, ("rev1.phi", "rev2.phi") )
# ModiaMath.plot(result, [("rev1.phi", "rev2.phi"), ("rev1.w"  , "rev2.w"), ("rev1.a"  , "rev2.a")])

println("... success of test_massComputation.jl!")
end


text:
module CURAND

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, libraryPropertyType, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK
using ..CUDA: libcurand, @retry_reclaim

using CEnum

using Memoize

using DataStructures

# core library
include("libcurand_common.jl")
include("error.jl")
include("libcurand.jl")

# low-level wrappers
include("wrappers.jl")

# high-level integrations
include("random.jl")

# cache for created, but unused handles
const rng_cache_lock = ReentrantLock()
const active_curand_rngs = Set{RNG}()
const active_gpuarray_rngs = Set{GPUArrays.RNG}()
const idle_curand_rngs = DefaultDict{CuContext,Vector{RNG}}(()->RNG[])
const idle_gpuarray_rngs = DefaultDict{CuContext,Vector{GPUArrays.RNG}}(()->GPUArrays.RNG[])

function default_rng()
    state = CUDA.active_state()
    rng, stream = get!(task_local_storage(), (:CURAND, state.context)) do
        new_rng = @lock rng_cache_lock begin
            new_rng = if isempty(idle_curand_rngs[state.context])
                RNG()
            else
                pop!(idle_curand_rngs[state.context])
            end

            # protect handles from the GC when the owning task is collected. we only
            # need to do this for CURAND, as handles typically don't have finalizers.
            push!(active_curand_rngs, new_rng)

            new_rng
        end

        finalizer(current_task()) do task
            @spinlock rng_cache_lock begin
                push!(idle_curand_rngs[state.context], new_rng)
                delete!(active_curand_rngs, new_rng)
            end
        end
        # TODO: curandDestroyGenerator to preserve memory, or at exit?

        curandSetStream(new_rng, state.stream)

        Random.seed!(new_rng)
        new_rng, state.stream
    end::Tuple{RNG,CuStream}

    if stream != state.stream
        curandSetStream(rng, state.stream)
        task_local_storage((:CURAND, state.context), (rng, state.stream))
    end

    return rng
end

function GPUArrays.default_rng(::Type{<:CuArray})
    ctx = context()
    get!(task_local_storage(), (:GPUArraysRNG, ctx)) do
        rng = @lock rng_cache_lock begin
            rng = if isempty(idle_gpuarray_rngs[ctx])
                dev = device()
                N = attribute(dev, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
                state = CuArray{NTuple{4, UInt32}}(undef, N)
                GPUArrays.RNG(state)
            else
                pop!(idle_gpuarray_rngs[ctx])
            end

            push!(active_gpuarray_rngs, rng)

            rng
        end

        finalizer(current_task()) do task
            @spinlock rng_cache_lock begin
                push!(idle_gpuarray_rngs[ctx], rng)
                delete!(active_gpuarray_rngs, rng)
            end
        end
        # TODO: destroy to preserve memory, or at exit?

        Random.seed!(rng)
        rng
    end::GPUArrays.RNG
end

@deprecate seed!() CUDA.seed!()
@deprecate seed!(seed) CUDA.seed!(seed)

end


text:
using FastScanning
using Test
using Unitful, AxisArrays, Images, ImagineInterface, ImagineProcedures
import Unitful:s, Hz, μm

ocpi2 = rigtemplate("ocpi-2")
pos_mon = getname(ocpi2, "axial piezo monitor")
pstart = 15.0μm
pstop = 115.0μm
stack_rate = 10.0Hz
z_spacing = 5.0μm
z_pad = 5.0μm
coms = FastScanning.pos_commands("ocpi-2", "axial piezo", pstart, pstop, stack_rate; sample_rate = 100000Hz)
pos = getpositioners(coms)[1]
fake_samps = get_samples(pos)
append!(pos_mon, "fake_samps", fake_samps) #perfect with no lag
pr = slicetiming_experiment(coms, [pos_mon;], "488nm laser", "camera1", z_spacing, z_pad; subpixel=false)
st_coms = outputs(pr)

lags = [0.0000s:0.0001s:0.002s...] #note: must match what's in experiments.jl
ntrials = 1

nconditions = length(lags)
slice_zs = ImagineInterface.slice_positions(pstart, pstop, z_spacing, z_pad)
nslices = length(slice_zs)
best_conds_fwd = fill(1, nslices)
best_conds_back = fill(2, nslices)
imgs = fake_slicetiming_run(nconditions, ntrials, nslices, best_conds_fwd, best_conds_back)

fwd_lags, back_lags, slice_zs2 = process(pr, imgs)
#@assert slice_zs2.==slice_zs #not equal due to finite sample rate
@test lags[best_conds_fwd] == fwd_lags
@test lags[best_conds_back] == back_lags

#### again with subpixel == true
pr = slicetiming_experiment(coms, [pos_mon;], "488nm laser", "camera1", z_spacing, z_pad; subpixel=true)
st_coms = outputs(pr)
imgs = fake_slicetiming_run(nconditions, ntrials, nslices, best_conds_fwd, best_conds_back)

fwd_lags, back_lags, slice_zs2 = process(pr, imgs)
@test lags[best_conds_fwd] == fwd_lags
@test lags[best_conds_back] == back_lags


text:
using HMCExamples

HMCExamples.main_FVGQ_2_joint(ARGS)

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
struct Currency{C} end

struct Cash{C<:Currency, T<:Real} <: Number
    value::T
end