text:
module PkgDashboards

using Pkg, PrettyTables, Weave, InteractiveUtils, Markdown

export dashboard, pkgdashboard

struct DashBoard
    markdowntable
    header
end

function markdownstring(db)
    io = IOBuffer()
    pretty_table(io, db.markdowntable, db.header, backend=:text, tf=markdown)
    md = String(take!(io))
end

function Base.show(io::IO, ::MIME"text/plain", db::DashBoard)
    pretty_table(io, db.markdowntable, db.header, backend=:text)
end

function Base.show(io::IO, m::MIME"text/markdown", db::DashBoard)
    md = markdownstring(db)
    MD = Markdown.parse(md)
    show(io,m,MD)
end

function Base.show(io::IO, m::MIME"text/html", db::DashBoard)
    md = markdownstring(db)
    MD = Markdown.parse(md)
    write(io, html(MD))
end

function getuser(ctx, uuid::Pkg.Types.UUID)
    urls = String[]
    for path = Pkg.Types.registered_paths(ctx.env, uuid)
        info = Pkg.Types.parse_toml(path, "Package.toml")
        repo = info["repo"]
        repo in urls || push!(urls, repo)
    end
    users = zip(getuser.(urls), urls)
end

function getuser(url::String)
    m = match(r"github\.com/(.+?)/", url)
    m === nothing ? nothing : string(m.captures[1])
end

"""
dashboard(target_users; output=:markdown, autoopen=false, stargazers=false, githubci=false, stars=true, activity=false)

Create a dashboard with all your badges!

### Arguments:
- `target_users`: a string or a vector of strings with github usernames or org names
- `output`: `:markdown` or `:html`
- `autoopen`: indicate whether or not the result is opened in editor or browser
- `githubci`: Also add a Badge for github CI
- `activity`: display number of commits per month
- `stars`: Add github star badge
- `stargazers`: include a plot with github stars over time

See also [`pkgdashboard`](@ref)
"""
function dashboard(target_users; kwargs...)
    target_users isa String && (target_users = [target_users])

    reg = Pkg.Types.collect_registries()[1]
    data = Pkg.Types.read_registry(joinpath(reg.path, "Registry.toml"))
    ctx = Pkg.Types.Context()

    markdownpage = []
    for (uuid, pkginfo) in data["packages"]
        name = pkginfo["name"]
        spec = PackageSpec(uuid=uuid)
        Uuid = Pkg.Types.UUID(uuid)
        users = getuser(ctx, Uuid)
        for (user,url) in users
            if user ∈ target_users
                entry = create_entry(user, name, url; kwargs...)
                push!(markdownpage, entry)
            end
        end
    end
    tab = table(markdownpage; kwargs...)
    write_output(tab; kwargs...)
    tab
end

"""
    pkgdashboard(packages; kwargs...)

Same as [`dashboard`](@ref) but accepts a list of package names instead.
"""
function pkgdashboard(packages; kwargs...)
    packages isa String && (packages = [packages])

    reg = Pkg.Types.collect_registries()[1]
    data = Pkg.Types.read_registry(joinpath(reg.path, "Registry.toml"))
    ctx = Pkg.Types.Context()

    markdownpage = []
    for (uuid, pkginfo) in data["packages"]
        name = pkginfo["name"]
        name ∈ packages || continue
        spec = PackageSpec(uuid=uuid)
        Uuid = Pkg.Types.UUID(uuid)
        users = getuser(ctx, Uuid)
        for (user,url) in users
            entry = create_entry(user, name, url; kwargs...)
            push!(markdownpage, entry)
        end
    end
    tab = table(markdownpage)
    write_output(tab; kwargs...)
    tab
end

function create_entry(user, name, url; output = :markdown, autoopen=false, stargazers=false, githubci=false, stars=true, activity=false)
    entry = ["[$(user)/$(name)]($(url))",
    "[![Build Status](https://travis-ci.org/$(user)/$(name).jl.svg?branch=master)](https://travis-ci.org/$(user)/$(name).jl)",
    "[![PkgEval](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/$(first(name))/$(name).svg)](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/report.html)",
    "[![codecov](https://codecov.io/gh/$(user)/$(name).jl/branch/master/graph/badge.svg)](https://codecov.io/gh/$(user)/$(name).jl)"]
    if stars
        push!(entry, "[![Stars](https://img.shields.io/github/stars/$(user)/$(name).jl.svg)](https://github.com/$(user)/$(name).jl/stargazers)")
    end
    if githubci
        push!(entry, "[![Build Status](https://github.com/$(user)/$(name).jl/workflows/CI/badge.svg)](https://github.com/$(user)/$(name).jl/actions)")
    end
    if activity
        push!(entry, "[![activity](https://img.shields.io/github/commit-activity/m/$(user)/$(name).jl)](https://github.com/$(user)/$(name).jl/pulse)")
    end
    if stargazers
        push!(entry, "[![Stargazers over time](https://starchart.cc/$(user)/$(name).jl.svg)](https://starchart.cc/$(user)/$(name).jl)")
    end

    entry
end

function table(markdownpage; kwargs...)
    if length(markdownpage) == 1
        markdowntable = copy(reshape(markdownpage[1], 1, :))
    else
        markdowntable = copy(permutedims(reduce(hcat, markdownpage), (2,1)))
    end
    I = sortperm(markdowntable[:,1])
    markdowntable = markdowntable[I,:]
    header = ["URL", "Build status", "PkgEval", "CodeCov"]
    get(kwargs, :stars, true) && push!(header, "Github stars")
    get(kwargs, :githubci, false) && push!(header, "Github CI")
    get(kwargs, :activity, false) && push!(header, "Commit activity")
    get(kwargs, :stargazers, false) && push!(header, "stargazers")
    DashBoard(markdowntable, header)
end

function write_output(db::DashBoard; kwargs...)

    path = mktempdir()
    markdownpath = joinpath(path, "dashboard.md")
    open(markdownpath, "w") do io
        pretty_table(io, db.markdowntable, db.header, backend=:text, tf=markdown)
    end

    if get(kwargs, :output, :markdown) == :html
        htmlpath = weave(markdownpath, doctype = "md2html")
        if get(kwargs, :autoopen, false)
            @async edit(htmlpath)
            @async run(`sensible-browser $htmlpath`)
        end
    else
        @info "Wrote markdown to $markdownpath"
        get(kwargs, :autoopen, false) && edit(markdownpath)
    end
end

end # module


text:
if haskey(ENV,"GROUP")
   group = ENV["GROUP"]
else
   group = "All"
end

is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )

if group == "All" || group == "Interface"
  tests = ["discontinuities.jl",
           "history_function.jl",
           "parameters.jl",
           "retcode.jl",
           "composite_solution.jl",
           "dependent_delays.jl",
           "reinit.jl",
           "saveat.jl",
           "save_idxs.jl",
           "events.jl",
           "units.jl",
           "unique_times.jl"
           ]
   for test in tests
       @time include(test)
   end
elseif group == "All" || group == "Integrators"
  tests =   ["constrained.jl",
             "unconstrained.jl",
             "residual_control.jl",
             "lazy_interpolants.jl",
             "sdirk_integrators.jl",
             "rosenbrock_integrators.jl"
             ]
   for test in tests
       @time include(test)
   end
end

ext: jl


text:
module Pathfinder

using AbstractDifferentiation: AD
using Distributions: Distributions
using Folds: Folds
# ensure that ForwardDiff is conditionally loaded by GalacticOptim
using ForwardDiff: ForwardDiff
using IrrationalConstants: log2π
using LinearAlgebra
using GalacticOptim: GalacticOptim
using GalacticOptimJL: GalacticOptimJL
using Optim: Optim, LineSearches
using PDMats: PDMats
using ProgressLogging: ProgressLogging
using PSIS: PSIS
using Random
using Requires: Requires
using SciMLBase: SciMLBase
using Statistics: Statistics
using StatsBase: StatsBase
using Transducers: Transducers
using UnPack: @unpack

export PathfinderResult, MultiPathfinderResult
export pathfinder, multipathfinder

# Note: we override the default history length to be shorter and the default line search
# to be More-Thuente, which keeps the approximate inverse Hessian positive-definite
const DEFAULT_HISTORY_LENGTH = 6
const DEFAULT_LINE_SEARCH = LineSearches.MoreThuente()
const DEFAULT_NDRAWS_ELBO = 5

function default_optimizer(history_length)
    return Optim.LBFGS(; m=history_length, linesearch=DEFAULT_LINE_SEARCH)
end

include("transducers.jl")
include("woodbury.jl")
include("optimize.jl")
include("inverse_hessian.jl")
include("mvnormal.jl")
include("elbo.jl")
include("resample.jl")
include("singlepath.jl")
include("multipath.jl")

function __init__()
    Requires.@require DynamicHMC = "bbc10e6e-7c05-544b-b16e-64fede858acb" begin
        include("integration/dynamichmc.jl")
    end
    Requires.@require AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d" begin
        include("integration/advancedhmc.jl")
    end
    Requires.@require Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" begin
        include("integration/turing.jl")
    end
end

end


text:
module MinimalPerfectHashes

export MinimalPerfectHash, hash

include("minperfhash.jl")

end


text:
@everywhere function testmodel(ydot, y, p, t)
    """
    Eq. 44-46 in http://onlinelibrary.wiley.com/doi/10.1002/aic.15577/full
    The data sets used for the parameter estimation, involving the time evolution
    of the dimensionless concentrations for species A and Q, are generated
    within the time interval  0 ≤ t ≤ 1 using the dynamic model.
    Subject to the initial conditions

    ya(0)=1
    yq(0)=0

    For the parameter values k₁=12, k₂=8, k₃=2 this is the test cases mentionrd
    in the original paper (ref 80).
    The “experimental” data used and the fit achieved with the method's
    obtained parameters are shown in Figure 8.
    Therefore, this demonstration involves a system of two ODEs, J = 2,
    zero algebraic equations, K = 0, with two sets of accompanying data, L = 2,
    and a total of three parameters to fit, N = 3.
    """
    k₁, k₂, k₃, k₄ = p

    ydot[1] = -k₁*y[1] + k₂*y[2]
    ydot[2] = k₁*y[1] - (k₂+k₃)*y[2] + k₄*y[3]
    ydot[3] = k₃*y[2] - k₄*y[3]
end

@everywhere function obj_func(dataset, model, params::Array{Float64, 1}, lbs, ubs)
    """
    :L: number of dataset
    :timepoints: timepoints of experimental data.
    :params: model parameters
    :eval_idx: indecies of components used for evaluating objective function value.
    """
    params = (lbs.*(ubs./lbs).^(params./10)) # rescaling parameters
    y0 = [1.0, 0.0, 0.0]
    alg = CVODE_BDF(max_hnil_warns = 0, max_error_test_failures = 0)

    L = 1
    rss_per_data = zeros(1)
    timepoints = linspace(0,1,51)
    tspan = (timepoints[1], timepoints[end])
    eval_idx = [1,2,3]
    res = zeros(length(timepoints), length(eval_idx))
    for i in 1:L   # doseの数だけ回す。
        prob = ODEProblem(model, y0, tspan, params)
        sol = solve(prob, alg, saveat=timepoints, maxiters=1e7,
                    abstol = 1e-8, reltol = 1e-8, verbose=false)
        if sol.retcode == :Success
            res[:] = vcat(sol.u'...)
            rss_per_data[i] = mean((dataset[i][:, eval_idx] - res[:, eval_idx]).^2)
        else
            rss_per_data[i] = -1
        end
    end

    if -1 in rss_per_data
        return Inf
    else
        rss = mean(rss_per_data)
        return sqrt(rss)
    end
end

@everywhere function obj_func2(dataset, model, params::Array{Float64, 1}, lbs, ubs)
    """
    :L: number of dataset
    :timepoints: timepoints of experimental data.
    :params: model parameters
    :eval_idx: indecies of components used for evaluating objective function value.
    """
    y0 = [1.0, 0.0, 0.0]
    alg = CVODE_BDF(max_hnil_warns = 0, max_error_test_failures = 0)

    L = 1
    rss_per_data = zeros(1)
    timepoints = linspace(0,1,51)
    tspan = (timepoints[1], timepoints[end])
    eval_idx = [1,2,3]
    res = zeros(length(timepoints), length(eval_idx))
    for i in 1:L   # doseの数だけ回す。
        prob = ODEProblem(model, y0, tspan, params)
        sol = solve(prob, alg, saveat=timepoints, maxiters=1e7,
                    abstol = 1e-8, reltol = 1e-8, verbose=false)
        if sol.retcode == :Success
            res[:] = vcat(sol.u'...)
            rss_per_data[i] = mean((dataset[i][:, eval_idx] - res[:, eval_idx]).^2)
        else
            rss_per_data[i] = -1
        end
    end

    if -1 in rss_per_data
        return Inf
    else
        rss = mean(rss_per_data)
        return sqrt(rss)
    end
end


text:
module Restacker

export restack

"""
    restack(x) -> x

An identity function in the sense `restack(x) === x`.  However, it
(recursively) re-construct `x` to beg the compiler to move everything
in the stack.
"""
restack
restack(x) = _restack(x, Val(false))

"""
    unsafe_restack(x) -> x

Like `restack(x)` but works for mutable objects, too. Thus, this is
the identity function only w.r.t (at most) `isequal`.
"""
unsafe_restack
unsafe_restack(x) = _restack(x, Val(true))

_getlength(::Type{<:NTuple{N,Any}}) where {N} = N
_getnames(::Type{<:NamedTuple{names}}) where {names} = names

# Not using dispatch, to check `issingletontype` first.
@generated function _restack(x, ::Val{hard}) where {hard}
    if Base.issingletontype(x)
        ex = x.instance::x
    elseif x <: Tuple
        N = _getlength(x)
        ex = :(($(map(i -> :(_restack(x[$i], Val{$hard}())), 1:N)...),))
    elseif x <: NamedTuple
        names = _getnames(x)
        ex = :((; $(map(n -> Expr(:kw, n, :(_restack(x.$n, Val{$hard}()))), names)...)))
    elseif parentmodule(x) === Core
        # Workaround the issue with `isstructtype(String)` and
        #  `isstructtype(Symbol)`.
        # https://github.com/JuliaLang/julia/issues/30210
        ex = :x
    elseif isstructtype(x)
        new = Expr(
            :new,
            x,
            map(
                n -> :(_restack(getfield(x, $(QuoteNode(n))), Val{$hard}())),
                fieldnames(x),
            )...,
        )
        if hard === true
            ex = new
        else
            ex = quote
                isimmutable(x) ? $new : x
            end
        end
    else
        ex = :x
    end
    return Expr(:block, Expr(:meta, :inline), ex)
end

end # module


text:
"""
    KernelMachine(kernel, augmenter, data, cs, dims)

Return a `KernelMachine`. `kernel` is required to have a method `kernel(u, v, k=nothing)`,
where `u` and `v` are matrices (datapoints are columns), and `k` is the kernel evaluated
in the previous hidden space. See [`additivegaussiankernel`](@ref) for an example.
`augmenter` is a matrix mapping the input to the machine space.
`cs` is the matrix of coefficient corresponding to `data`. Finally, `dims` is a tuple
containing the dimensionality of the spaces that compose the machine space. In particular,
`size(augmenter, 1) == sum(dims)`.

Given a kernel machine `km`, `km(v)` returns the value of the stable state of the machine
on the last component of machine space, as well as the square norm of the machine plus the
square Frobenius norm of the augmenter. This second quantity should be used to regularize
the machine.

    KernelMachine(kernel, data::M; dims, init=glorot_uniform)

A simpler constructor. `kernel` is option and defaults to [`additivegaussiankernel`](@ref).
`data` is the training data, or, more generally, the set of anchor points (datapoints are columns).
`dims` is a tuple containing the dimensionality of the spaces that compose the machine space.
`init` is used to initialize the `cs` matrix.
"""
struct KernelMachine{K, M<:AbstractMatrix, N}
    kernel::K
    augmenter::M
    data::M
    cs::M
    dims::NTuple{N, Int}
end

function KernelMachine(kernel, data::M; dims, init=glorot_uniform) where M
    augmenter_size = sum(dims)
    cs_size = augmenter_size - first(dims)
    augmenter::M = init(augmenter_size, size(data, 1))
    cs::M = init(cs_size, size(data, 2))
    return KernelMachine(kernel, augmenter, data, cs, dims)
end

function KernelMachine(data; kwargs...)
    KernelMachine(additivegaussiankernel, data; kwargs...)
end

Base.show(io::IO, d::KernelMachine) = print(io, "KernelMachine {...}")

# TODO: what if `kernel` is operator-valued?
# If it is a relevant use case, it would require a different kernel API.
function consume(kernel, xss, css, indices=axes(first(xss), 2))
    res = first(xss)
    k, sn = nothing, zero(eltype(res))
    # Iteratively update k (kernel), res (result so far), and sn (square norm)
    for (xs, cs) in zip(tail(xss), css)
        k = kernel(res[:, indices], res, k)
        val = cs * k
        res = xs + val
        sn += dot(val[:, indices], cs)
    end
    return res, sn
end

maybe_hcat(x, ::Nothing) = x
maybe_hcat(x, y) = hcat(x, y)

function (dm::KernelMachine)(input)
    kernel, augmenter, data, cs, dims =
        dm.kernel, dm.augmenter, dm.data, dm.cs, dm.dims
    full = maybe_hcat(data, input)
    xs = augmenter * full
    cost = sum(abs2, augmenter)
    xss, css = split_matrix(xs, dims), split_matrix(cs, tail(dims))
    res, sn = consume(kernel, xss, css, axes(data, 2))
    return res, cost + sn
end


text:
@testset "buffer" begin

@testset "SumTree" begin
    t = SumTree(3)
    @test capacity(t) == 3
    @test length(t) == 0

    push!(t, 1)
    @test capacity(t) == 3
    @test length(t) == 1
    @test size(t) == (1,)
    @test t[1] == 1

    t[1] = 2
    @test capacity(t) == 3
    @test length(t) == 1
    @test t[1] == 2

    push!(t, 3)
    @test capacity(t) == 3
    @test length(t) == 2
    @test t[1:end] == [2, 3]
    @test t.tree == [5., 5., 0., 2, 3, 0]
    @test Ju.get(t, 0) == (1, 2.)
    @test Ju.get(t, 1) == (1, 2.)
    @test Ju.get(t, 2) == (1, 2.)
    @test Ju.get(t, 3) == (2, 3.)
    @test Ju.get(t, 4) == (2, 3.)
    @test Ju.get(t, 5) == (2, 3.)

    t[1] = 1
    t[2] = 2
    @test capacity(t) == 3
    @test length(t) == 2
    @test t[1:end] == [1, 2]
    @test t.tree == [3., 3., 0., 1, 2, 0]
    @test Ju.get(t, 0) == (1, 1.)
    @test Ju.get(t, 1) == (1, 1.)
    @test Ju.get(t, 2) == (2, 2.)
    @test Ju.get(t, 3) == (2, 2.)

    push!(t, 3)
    @test capacity(t) == 3
    @test length(t) == 3
    @test t[1:end] == [1, 2, 3]
    @test t.tree == [6., 3., 3., 1., 2., 3.]
    @test Ju.get(t, 0) == (1, 1.)
    @test Ju.get(t, 1) == (1, 1.)
    @test Ju.get(t, 2) == (2, 2.)
    @test Ju.get(t, 3) == (2, 2.)
    @test Ju.get(t, 4) == (3, 3.)
    @test Ju.get(t, 5) == (3, 3.)
    @test Ju.get(t, 6) == (3, 3.)

    push!(t, 4)
    @test capacity(t) == 3
    @test length(t) == 3
    @test t[1:end] == [2, 3, 4]
    @test t.tree == [9., 6., 3., 4, 2, 3]
    @test Ju.get(t, 0) == (3, 4.)
    @test Ju.get(t, 1) == (3, 4.)
    @test Ju.get(t, 2) == (3, 4.)
    @test Ju.get(t, 3) == (3, 4.)
    @test Ju.get(t, 4) == (3, 4.)
    @test Ju.get(t, 5) == (1, 2.)
    @test Ju.get(t, 6) == (1, 2.)
    @test Ju.get(t, 9) == (2, 3.)
    @test Ju.get(t, 9) == (2, 3.)
    @test Ju.get(t, 9) == (2, 3.)

    empty!(t)
    @test capacity(t) == 3
    @test length(t) == 0

    push!(t, 1)
    push!(t, 2)
    push!(t, 3)
    @test capacity(t) == 3
    @test length(t) == 3
    @test t[1:end] == [1, 2, 3]
    @test t.first == 2
    @test t.tree == [6., 4., 2., 3., 1., 2.]
    @test Ju.get(t, 0) == (3, 3.)
    @test Ju.get(t, 1) == (3, 3.)
    @test Ju.get(t, 2) == (3, 3.)
    @test Ju.get(t, 3) == (3, 3.)
    @test Ju.get(t, 4) == (1, 1.)
    @test Ju.get(t, 5) == (2, 2.)
    @test Ju.get(t, 6) == (2, 2.)
end

@testset "CircularArrayBuffer" begin
    A = ones(2, 2)
    @testset "1D Int" begin
        b = CircularArrayBuffer{Int}(3)

        @test eltype(b) == Int
        @test capacity(b) == 3
        @test isfull(b) == false
        @test isempty(b) == true
        @test length(b) == 0
        @test size(b) == (0,)
        # element must has the exact same length with the element of buffer
        @test_throws DimensionMismatch push!(b, [1, 2])

        for x in 1:3 push!(b, x) end

        @test capacity(b) == 3
        @test isfull(b) == true
        @test length(b) == 3
        @test size(b) == (3,)
        @test b[1] == 1
        @test b[end] == 3
        @test b[1:end] == [1, 2, 3]

        for x in 4:5 push!(b, x) end

        @test capacity(b) == 3
        @test length(b) == 3
        @test size(b) == (3,)
        @test b[1] == 3
        @test b[end] == 5
        @test b[1:end] == [3, 4, 5]

        empty!(b)
        @test isfull(b) == false
        @test isempty(b) == true
        @test length(b) == 0
        @test size(b) == (0,)

        push!(b, 6)
        @test isfull(b) == false
        @test isempty(b) == false
        @test length(b) == 1
        @test size(b) == (1,)
        @test b[1] == 6

        push!(b, 7)
        push!(b, 8)
        @test isfull(b) == true
        @test isempty(b) == false
        @test length(b) == 3
        @test size(b) == (3,)
        @test b[[1,2,3]] == [6, 7, 8]

        push!(b, 9)
        @test isfull(b) == true
        @test isempty(b) == false
        @test length(b) == 3
        @test size(b) == (3,)
        @test b[[1,2,3]] == [7, 8, 9]
    end

    @testset "2D Float64" begin
        b = CircularArrayBuffer{Array{Float64,2}}(3, (2, 2))

        @test eltype(b) == Array{Float64, 2}
        @test capacity(b) == 3
        @test isfull(b) == false
        @test length(b) == 0
        @test size(b) == (2, 2, 0)
        # element must has the exact same length with the element of buffer
        @test_throws DimensionMismatch push!(b, [1., 2.])

        for x in 1:3 push!(b, x * A) end

        @test capacity(b) == 3
        @test isfull(b) == true
        @test length(b) == 3
        @test size(b) == (2, 2, 3)
        for i in 1:3 @test b[i] == i * A end
        @test b[end] == 3 * A

        for x in 4:5 push!(b, x * ones(4)) end  # collection is also OK

        @test capacity(b) == 3
        @test length(b) == 3
        @test size(b) == (2, 2, 3)
        @test b[1] == 3 * A
        @test b[end] == 5 * A

        @test b == reshape([c for x in 3:5 for c in x*A], 2, 2, 3)
    end
end

@testset "CircularTurnBuffer" begin
    @testset "CircularSARDBuffer" begin
        b = CircularSARDBuffer(3)
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true
        @test capacity(b) == 3

        push!(b, 1, 1)
        push!(b, 1., false, 2, 2)
        @test length(b) == 1
        @test size(b) == (1,)
        @test isfull(b) == false
        @test isempty(b) == false

        push!(b, 2., false, 3, 3)
        push!(b, 3., false, 4, 4)
        @test length(b) == 3
        @test isfull(b) == true
        @test isempty(b) == false
        @test b[1] == (state=1, action=1, reward=1., isdone=false, nextstate=2, nextaction=2)
        @test b[end] == (state=3, action=3, reward=3., isdone=false, nextstate=4, nextaction=4)

        push!(b, 4., false, 5, 5)
        @test length(b) == 3
        @test isfull(b) == true
        @test b[1] == (state=2, action=2, reward=2., isdone=false, nextstate=3, nextaction=3)
        @test b[end] == (state=4, action=4, reward=4., isdone=false, nextstate=5, nextaction=5)
    end
    @testset "CircularSARDBuffer of 2D element type" begin
        b = CircularSARDBuffer(3; state_type=Array{Int,2}, state_size=(2,2))
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true
        @test capacity(b) == 3

        push!(b, [1 1; 1 1], 1)
        push!(b, 1., false, [2 2; 2 2], 2)
        @test length(b) == 1
        @test isfull(b) == false
        @test isempty(b) == false

        push!(b, 2., false, [3 3; 3 3], 3)
        push!(b, 3., false, [4 4; 4 4], 4)
        @test length(b) == 3
        @test isfull(b) == true
        @test isempty(b) == false
        @test b[1] == (state=[1 1; 1 1], action=1, reward=1., isdone=false, nextstate=[2 2; 2 2], nextaction=2)
        @test b[end] == (state=[3 3; 3 3], action=3, reward=3., isdone=false, nextstate=[4 4; 4 4], nextaction=4)

        push!(b, 4., true, [5 5; 5 5], 5)
        @test length(b) == 3
        @test isfull(b) == true
        @test b[1] == (state=[2 2; 2 2], action=2, reward=2., isdone=false, nextstate=[3 3; 3 3], nextaction=3)
        @test b[end] == (state=[4 4; 4 4], action=4, reward=4., isdone=true, nextstate=[5 5; 5 5], nextaction=5)

        # now push a new state-action
        push!(b, [6 6; 6 6], 6)
        @test length(b) == 3
        @test isfull(b) == true
        @test b[1] == (state=[2 2; 2 2], action=2, reward=2., isdone=false, nextstate=[3 3; 3 3], nextaction=3)
        @test b[end] == (state=[4 4; 4 4], action=4, reward=4., isdone=true, nextstate=[6 6; 6 6], nextaction=6)

        batch_inds = [3, 2, 1]
        batch_state = reshape([4 4 4 4 3 3 3 3 2 2 2 2], 2, 2, 3)
        batch_action = [4, 3, 2]
        batch_reward = [4., 3., 2.]
        batch_isdone = [true, false, false]
        batch_nextstate = reshape([6 6 6 6 4 4 4 4 3 3 3 3], 2, 2, 3)
        @test Ju.batch_view(b, batch_inds) == (batch_state, batch_action, batch_reward, batch_isdone, batch_nextstate)
    end
    @testset "CircularSARDSBuffer" begin
        b = CircularSARDSBuffer(3)
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true
        @test capacity(b) == 3

        push!(b, 1, 1, 1.0, false, 2)
        push!(b, 3, 3, 3.0, false, 4)
        push!(b, 5, 5, 5.0, false, 6)
        @test length(b) == 3
        @test isfull(b) == true
        @test isempty(b) == false
        @test b[1] == (state=1, action=1, reward=1.0, isdone=false, nextstate=2)
        @test b[end] == (state=5, action=5, reward=5.0, isdone=false, nextstate=6)

        push!(b, 7, 7, 7.0, false, 8)
        @test b[1] == (state=3, action=3, reward=3.0, isdone=false, nextstate=4)
        @test b[end] == (state=7, action=7, reward=7.0, isdone=false, nextstate=8)
    end
    @testset "CircularSARDSABuffer" begin
        b = CircularSARDSABuffer(3)
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true
        @test capacity(b) == 3

        push!(b, 1, 1, 1.0, false, 2, 2)
        push!(b, 3, 3, 3.0, false, 4, 4)
        push!(b, 5, 5, 5.0, false, 6, 6)
        @test length(b) == 3
        @test isfull(b) == true
        @test isempty(b) == false
        @test b[1] == (state=1, action=1, reward=1.0, isdone=false, nextstate=2, nextaction=2)
        @test b[end] == (state=5, action=5, reward=5.0, isdone=false, nextstate=6, nextaction=6)

        push!(b, 7, 7, 7.0, false, 8, 8)
        @test b[1] == (state=3, action=3, reward=3.0, isdone=false, nextstate=4, nextaction=4)
        @test b[end] == (state=7, action=7, reward=7.0, isdone=false, nextstate=8, nextaction=8)
    end
end

@testset "EpisodeTurnBuffer" begin
    @testset "EpisodeSARDBuffer" begin
        b = EpisodeSARDBuffer()
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true
        @test capacity(b) == typemax(Int)

        push!(b, 1, 1)
        push!(b, 1., false, 2, 2)
        @test length(b) == 1
        @test isfull(b) == false
        @test isempty(b) == false

        # we can also push a turn infor together,
        # notice that state and action will be iginored
        push!(b, 2., true, 3, 3)
        @test length(b) == 2
        @test isfull(b) == true
        @test isempty(b) == false
        @test b[end] == (state=2, action=2, reward=2., isdone=true)

        push!(b, 3, 3)
        push!(b, 3., false, 4, 4)
        @test length(b) == 1
        @test isfull(b) == false
        @test isempty(b) == false
        @test b[end] == (state=3, action=3, reward=3., isdone=false)
    end
    @testset "EpisodeSARDSBuffer" begin
        b = EpisodeSARDSBuffer()
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true

        push!(b, 1, 1, 1.0, false, 2)
        push!(b, 3, 3, 3.0, false, 4)
        push!(b, 5, 5, 5.0, false, 6)
        @test length(b) == 3
        @test isfull(b) == false
        @test isempty(b) == false
        @test b[1] == (state=1, action=1, reward=1.0, isdone=false, nextstate=2)
        @test b[end] == (state=5, action=5, reward=5.0, isdone=false, nextstate=6)

        push!(b, 7, 7, 7.0, true, 8)
        @test length(b) == 4
        @test isfull(b) == true
        @test b[1] == (state=1, action=1, reward=1.0, isdone=false, nextstate=2)
        @test b[end] == (state=7, action=7, reward=7.0, isdone=true, nextstate=8)

        push!(b, 9, 9, 9.0, false, 10)
        @test length(b) == 1
        @test isfull(b) == false
        @test b[1] == b[end] == (state=9, action=9, reward=9.0, isdone=false, nextstate=10)
    end
    @testset "EpisodeSARDSABuffer" begin
        b = EpisodeSARDSABuffer()
        @test length(b) == 0
        @test isfull(b) == false
        @test isempty(b) == true

        push!(b, 1, 1, 1.0, false, 2, 2)
        push!(b, 3, 3, 3.0, false, 4, 4)
        push!(b, 5, 5, 5.0, false, 6, 6)
        @test length(b) == 3
        @test isfull(b) == false
        @test isempty(b) == false
        @test b[1] == (state=1, action=1, reward=1.0, isdone=false, nextstate=2, nextaction=2)
        @test b[end] == (state=5, action=5, reward=5.0, isdone=false, nextstate=6, nextaction=6)

        push!(b, 7, 7, 7.0, true, 8, 8)
        @test length(b) == 4
        @test isfull(b) == true
        @test b[1] == (state=1, action=1, reward=1.0, isdone=false, nextstate=2, nextaction=2)
        @test b[end] == (state=7, action=7, reward=7.0, isdone=true, nextstate=8, nextaction=8)

        push!(b, 9, 9, 9.0, false, 10, 10)
        @test length(b) == 1
        @test isfull(b) == false
        @test b[1] == b[end] == (state=9, action=9, reward=9.0, isdone=false, nextstate=10, nextaction=10)
    end
end

end


text:
@rule GammaMixture((:a, k), Marginalisation) (q_out::Any, q_switch::Any, q_b::GammaDistributionsFamily) = begin
    p = probvec(q_switch)[k]
    β = logmean(q_out) + logmean(q_b)
    γ = p * β
    return GammaShapeLikelihood(p, γ)
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.logic;
#
# import typedb.client.api.TypeDBTransaction;
# import graql.lang.pattern.Pattern;
#
# import javax.annotation.CheckReturnValue;
#
# public interface Rule {
#
#     @CheckReturnValue
#     String getLabel();
#
#     @CheckReturnValue
#     Pattern getWhen();
#
#     @CheckReturnValue
#     Pattern getThen();
#
#     @CheckReturnValue
#     Rule.Remote asRemote(TypeDBTransaction transaction);
#
#     @CheckReturnValue
#     boolean isRemote();
#
#     interface Remote extends Rule {
#
#         void setLabel(String label);
#
#         void delete();
#
#         @CheckReturnValue
#         boolean isDeleted();
#     }
# }


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(scene, alignmode = Outside(30))

    las = Array{LAxis, 2}(undef, 4, 4)

    for i in 1:4, j in 1:4
        las[i, j] = maingl[i, j] = LAxis(scene)
    end
end

las[4, 1].attributes.aspect = AxisAspect(1)
las[4, 2].attributes.aspect = AxisAspect(2)
las[4, 3].attributes.aspect = AxisAspect(0.5)
las[4, 4].attributes.aspect = nothing
las[1, 1].attributes.maxsize = (Inf, Inf)
las[1, 2].attributes.aspect = nothing
las[1, 3].attributes.aspect = nothing

begin
    subgl = gridnest!(maingl, 1, 1)
    cb1 = subgl[:, 2] = LColorbar(scene, width=30, height=Relative(0.66))
    sleep(0.5)

    subgl2 = gridnest!(maingl, 1:2, 1:2)
    cb2 = subgl2[:, 3] = LColorbar(scene, width=30, height=Relative(0.66))
    sleep(0.5)

    subgl3 = gridnest!(maingl, 1:3, 1:3)
    cb3 = subgl3[:, 4] = LColorbar(scene, width=30, height=Relative(0.66))
    sleep(0.5)

    subgl4 = gridnest!(maingl, 1:4, 1:4)
    cb4 = subgl4[:, 5] = LColorbar(scene, width=30, height=Relative(0.66))
    sleep(0.5)

end


text:
using SimilaritySearch, Random

function main(k, n=30_000_000)
    Random.seed!(0)

    res = KnnResult(k)
    push!(res, 0 => rand())
    @time for i in 1:n
        push!(res, i => rand())
    end
end

main(10)


text:
export FourierField, FourierMetaData, Transform
export plot, forward, backward, norm
import LinearAlgebra: norm
struct FourierField{D,S}
    data::D
    metadata::S
end

struct FourierMetaData{𝒩, 𝒢, 𝒯}
    name::𝒩
    grid::𝒢
    transform::𝒯
end

function Base.show(io::IO, ϕ::FourierField{S,T}) where {S, T <: FourierMetaData}
    printstyled(io, ϕ.metadata.name, color = 128 )
end

function forward(ϕ::FourierField{S,T}) where {S, T <: FourierMetaData}
    ϕ̂ = ϕ.metadata.transform.forward * ϕ.data
    fmd = FourierMetaData(nothing, ϕ.metadata.grid, ϕ.metadata.transform)
    return FourierField(ϕ̂, fmd)
end

function backward(ϕ::FourierField{S,T}) where{S, T <: FourierMetaData}
    ϕ̂ = ϕ.metadata.transform.backward * ϕ.data
    fmd = FourierMetaData(nothing, ϕ.metadata.grid, ϕ.metadata.transform)
    return FourierField(ϕ̂, fmd)
end

# convenience function to initialize array
function (ϕ::FourierField)(a::AbstractArray)
    f1 = ϕ.metadata.transform.forward * (a .+ 0im*ϕ.data)
    ϕ.data .= f1
    return nothing
end

norm(ϕ::FourierField) = norm(ϕ.data)

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
import LinearAlgebra: I, Bidiagonal, tr
import ForneyLab: unsafeCov, unsafeMean, unsafePrecision, VariateType,
				  collectNaiveVariationalNodeInbounds, assembleClamp!, ultimatePartner
include("util.jl")

export ruleVariationalGFXOutNPPPPP,
       ruleVariationalGFXIn1PNPPPP,
       ruleVariationalGFXIn2PPNPPP,
       ruleVariationalGFXIn3PPPNPP,
	   ruleVariationalGFXIn4PPPPNP,
	   ruleVariationalGFXIn5PPPPPN

function ruleVariationalGFXOutNPPPPP(Δt :: Float64,
									 marg_y :: Nothing,
									 marg_θ :: ProbabilityDistribution{Multivariate},
                                     marg_x :: ProbabilityDistribution{Multivariate},
                                     marg_η :: ProbabilityDistribution{Univariate},
                                     marg_u :: ProbabilityDistribution{Univariate},
                                     marg_τ :: ProbabilityDistribution{Univariate})

	# Extract moments of beliefs
	mθ = unsafeMean(marg_θ)
	mx = unsafeMean(marg_x)
	mη = unsafeMean(marg_η)
	mu = unsafeMean(marg_u)
	mτ = unsafeMean(marg_τ)

	# Set order of system
	order = dims(marg_θ)

	# Structure matrices
	S = Bidiagonal(ones(order,), Δt.*ones(order-1,), :U)
	s = zeros(order,); s[end] = Δt

	# Generate precision matrix
	EW = mτ*inv(noisecov(Δt, dims=order))

	# Compute transition matrices
	EA = S + s*mθ'
	EB = s*mη

	# Set outgoing message
	return Message(Multivariate, GaussianMeanPrecision, m=EA*mx + EB*mu, w=EW)
end

function ruleVariationalGFXIn1PNPPPP(Δt :: Float64,
									 marg_y :: ProbabilityDistribution{Multivariate},
                                     marg_θ :: Nothing,
									 marg_x :: ProbabilityDistribution{Multivariate},
                                     marg_η :: ProbabilityDistribution{Univariate},
                                     marg_u :: ProbabilityDistribution{Univariate},
                                     marg_τ :: ProbabilityDistribution{Univariate})

	# Extract moments of beliefs
	my = unsafeMean(marg_y)
	mx,Vx = unsafeMeanCov(marg_x)
	mη = unsafeMean(marg_η)
	mu = unsafeMean(marg_u)
	mτ = unsafeMean(marg_τ)

	# Set order of system
	order = dims(marg_x)

	# Structure matrices
	S = Bidiagonal(ones(order,), Δt.*ones(order-1,), :U)
	s = zeros(order,); s[end] = Δt

	# Generate precision matrix
	EW = mτ*inv(noisecov(Δt, dims=order))

	# Compute transition matrix
	EB = s*mη

	# Set parameters
	# ϕ = mx*s'*EW*(my - EB*mu) - EW*S'*(mx*mx' + Vx)*s
	# ϕ = mx*s'*EW*(my - EB*mu) - (mx*mx' + Vx)*S'*EW*s
	ϕ = mx*s'*EW*(my - EB*mu) - (mx*s'*EW*S*mx + Vx*S'*EW*s)

	# Φ = EW*(s'*(mx*mx' + Vx)*s)
	# Φ = (s'*s)*EW*(mx*mx' + Vx)
	Φ = mx*s'*EW*s*mx' + Vx*(s'*EW*s)

	# Set outgoing message
	return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalGFXIn2PPNPPP(Δt :: Float64,
									 marg_y :: ProbabilityDistribution{Multivariate},
									 marg_θ :: ProbabilityDistribution{Multivariate},
									 marg_x :: Nothing,
									 marg_η :: ProbabilityDistribution{Univariate},
								     marg_u :: ProbabilityDistribution{Univariate},
                                     marg_τ :: ProbabilityDistribution{Univariate})

   	# Extract moments of beliefs
	my = unsafeMean(marg_y)
	mθ,Vθ = unsafeMeanCov(marg_θ)
	mη = unsafeMean(marg_η)
	mu = unsafeMean(marg_u)
	mτ = unsafeMean(marg_τ)

	# Set order of system
	order = dims(marg_θ)

	# Structure matrices
	S = Bidiagonal(ones(order,), Δt.*ones(order-1,), :U)
	s = zeros(order,); s[end] = Δt

	# Generate precision matrix
	EW = mτ*inv(noisecov(Δt, dims=order))

	# Compute transition matrices
	EA = S + s*mθ'
	EB = s*mη

	# Set parameters
	ϕ = EA'*EW*(my - EB*mu)
	# Φ = S'*EW*S + S'*EW*s*mθ' + mθ*s'*EW*S + s'*EW*s*(Vθ + mθ*mθ')
	Φ = EA'*EW*EA + (s'*EW*s)*Vθ

	# Set outgoing message
	return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalGFXIn3PPPNPP(Δt :: Float64,
									 marg_y :: ProbabilityDistribution{Multivariate},
									 marg_θ :: ProbabilityDistribution{Multivariate},
									 marg_x :: ProbabilityDistribution{Multivariate},
							  	     marg_η :: Nothing,
								     marg_u :: ProbabilityDistribution{Univariate},
                                     marg_τ :: ProbabilityDistribution{Univariate})

	# Extract moments of beliefs
	my = unsafeMean(marg_y)
	mθ = unsafeMean(marg_θ)
	mx = unsafeMean(marg_x)
	mu = unsafeMean(marg_u)
	mτ = unsafeMean(marg_τ)

	# Set order
	order = dims(marg_θ)

	# Structure matrices
	S = Bidiagonal(ones(order,), Δt.*ones(order-1,), :U)
	s = zeros(order,); s[end] = Δt

	# Compute expected values
	EA = S + s*mθ'

	# Generate precision matrix
	EW = mτ*inv(noisecov(Δt, dims=order))

	# Set parameters
	ϕ = mu'*s'*EW*(my - EA*mx)
	Φ = mu'*(s'*EW*s)*mu

	# Set outgoing message
	return Message(Univariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalGFXIn4PPPPNP(Δt :: Float64,
									 marg_y :: ProbabilityDistribution{Multivariate},
									 marg_θ :: ProbabilityDistribution{Multivariate},
									 marg_x :: ProbabilityDistribution{Multivariate},
									 marg_η :: ProbabilityDistribution{Univariate},
								     marg_u :: Nothing,
                                     marg_τ :: ProbabilityDistribution{Univariate})

   error("The input u_k should always be clamped.")
end

function ruleVariationalGFXIn5PPPPPN(Δt :: Float64,
									 marg_y :: ProbabilityDistribution{Multivariate},
									 marg_θ :: ProbabilityDistribution{Multivariate},
									 marg_x :: ProbabilityDistribution{Multivariate},
									 marg_η :: ProbabilityDistribution{Univariate},
								     marg_u :: ProbabilityDistribution{Univariate},
                                     marg_τ :: Nothing)

	# Extract moments of beliefs
	my,Vy = unsafeMeanCov(marg_y)
	mθ,Vθ = unsafeMeanCov(marg_θ)
	mx,Vx = unsafeMeanCov(marg_x)
	mη,Vη = unsafeMeanCov(marg_η)
	mu,Vu = unsafeMeanCov(marg_u)

	# Set order
	order = dims(marg_θ)

	# Structure matrices
	S = Bidiagonal(ones(order,), Δt.*ones(order-1,), :U)
	s = zeros(order,); s[end] = Δt

	# Compute expected values
	EA = S + s*mθ'
	EB = s*mη

	t1 = my[order]^2 + Vy[order,order]
	t2 = -my[order]*(EA*mx)[order]
	t3 = -my[order]*(EB*mu)[order]
	t4 = -(EA*mx)[order]*my[order]
	t5 = mx'*(EA'*EA + Vθ)*mx + tr((EA'*EA + Vθ)*Vx)
	t6 = (EA*mx)[order]*(EB*mu)[order]
	t7 = -(EB*mu)[order]*my[order]
	t8 = (EB*mu)[order]*(EA*mx)[order]
	t9 = mη'*s'*(mu*mu' + Vu)*s*mη + s'*(mu*mu' + Vu)*s*Vη

	# Set parameters
	a = 3/2
	b = t1 + t2 + t3 + t4 + t5 + t6 + t7 + t8 + t9

	# Set outgoing message
	return Message(Univariate, Gamma, a=a, b=b)
end

function collectNaiveVariationalNodeInbounds(node::GeneralisedFilterX, entry::ScheduleEntry)
	inbounds = Any[]

	# Push function to calling signature (Δt needs to be defined in user scope)
	push!(inbounds, Dict{Symbol, Any}(:Δt => node.Δt, :keyword => false))

    target_to_marginal_entry = currentInferenceAlgorithm().target_to_marginal_entry

    for node_interface in entry.interface.node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if node_interface === entry.interface
            # Ignore marginal of outbound edge
            push!(inbounds, nothing)
        elseif (inbound_interface != nothing) && isa(inbound_interface.node, Clamp)
            # Hard-code marginal of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, ProbabilityDistribution))
        else
            # Collect entry from marginal schedule
            push!(inbounds, target_to_marginal_entry[node_interface.edge.variable])
        end
    end

    return inbounds
end

max_stars_repo_name: wmkouw/gfx-node

max_issues_repo_name: wmkouw/GFX

max_forks_repo_name: wmkouw/GFX


text:
# This file contains a function to read input data from a file

"""
    @skip()

Skip two lines and check for EOF each time.

# Usage

```julia
using Scats.Internal.Input: @skip

@skip() #= ===
begin
    readline(io)
    if eof(io)
        throw(ScatsInputEOF(file))
    end
    readline(io)
    if eof(io)
        throw(ScatsInputEOF(file))
    end
end =#
```
"""
macro skip()
    return esc(quote
        readline(io)
        if eof(io)
            throw(ScatsInputEOF(file))
        end
        readline(io)
        if eof(io)
            throw(ScatsInputEOF(file))
        end
    end)
end

"""
    @read(element::AbstractString, type::AbstractString)

Read a scalar element.

# Usage

```julia
using Scats.Internal.Input: @read

@read("N", "IT") #= ===
begin
    try
        Input.N = parse(IT, split(readline(io))[1])
    catch
        throw(ScatsInputWR_N(file))
    end
end =#
```
"""
macro read(element::AbstractString, type::AbstractString)
    return esc(
        Meta.parse("""
        try
            Input.$element = parse($type, split(readline(io))[1])
        catch
            throw(ScatsInputWR_$element(file))
        end
        """)
    )
end

"""
    @read_array(element::AbstractString, type::AbstractString)

Read an array element.

# Usage

```julia
using Scats.Internal.Input: @read_array

@read_array("t", "RT") #= ===
begin
    try
        Input.t = parse.(RT, split(readline(io))[1:Input.N])
    catch
        throw(ScatsInputWR_t(file))
    end
end =#
```
"""
macro read_array(element::AbstractString, type::AbstractString)
    return esc(
        Meta.parse("""
        try
            Input.$element = parse.($type, split(readline(io))[1:Input.N])
        catch
            throw(ScatsInputWR_$element(file))
        end
        """)
    )
end

"""
    read!(Input::InputStruct, file::AbstractString)

Read input data from a file to an instance of [`InputStruct`](@ref).

# Usage

```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Input.example(file)
s.Input.read!(file)

# output

```
"""
function read!(Input::InputStruct, file::AbstractString)

    # Strip the string
    file = strip(file)

    # Check if the file exists
    if !isfile(file)
        throw(ScatsInputNotAFile(file))
    end

    # Open the file for reading
    open(file, "r") do io

        # Check for EOF
        if eof(io)
            throw(ScatsInputEOF(file))
        end

        readline(io)

        # Check for EOF
        if eof(io)
            throw(ScatsInputEOF(file))
        end

        @read "N" "IT"
        @skip

        @read "Δt" "RT"
        @skip

        @read "q" "RT"
        @skip

        @read_array "t" "RT"
        @skip

        @read_array "x" "RT"

    end

    return nothing

end


text:
push!(LOAD_PATH, "../src")
using MT1DGeneticInversion
using MT1D
using Gadfly

"""
Forward models the apparent resistivity and phase of a 1D model and plots the
results using Gadfly. The plot is saved to `forward-modelled-response.html`.

## Arguments:
- `model` -- A matrix with depths in the first column and resistivities in the
  second. Ignore to use a default model.
"""
function testForwardModel(model=[])
    fs = logspace(-3, 3, 1024)
    if isempty(model)
        model = [0. 10;
                 1_000 100]
    end

    println("Forward modelling...")
    @time ρ, Φ = calc_response(model, fs)
    println("Done!")
    pRes = plot(x=fs.^-1,
                y=ρ,
                Geom.line,
                Guide.xlabel("Period (S)"),
                Guide.ylabel("ρa (Ωm)"),
                Guide.title("Apparent Resistivity"),
                Scale.x_log10,
                Scale.y_log10)
    pPhase = plot(x=fs.^-1,
                  y=Φ,
                  Geom.line,
                  Guide.xlabel("Period (S)"),
                  Guide.ylabel("Φ (°)"),
                  Guide.title("Phase"),
                  Scale.x_log10,
                  Scale.y_continuous(minvalue=0,maxvalue=90))
    pModel = plot(y=[model[:,2]; model[end,2]],
                  x=[model[:,1] / 1000; model[end,1] / 1000 + 0.5],
                  Stat.step(direction=:hv),
                  Geom.line,
                  Guide.ylabel("Resistivity (Ωm)"),
                  Guide.xlabel("Depth (km)"),
                  Guide.title("Model"),
                  Scale.y_log10)

    draw(SVG("forward-modelled-response.svg", 800px, 800px), hstack(vstack(pRes, pPhase), pModel))
    return
end

"""
Invert the data in the file ./SNO-96-106.data using the default inversion
parameters. This is data collected over the northern Canadian shield as part of
Lithoprobes's SNORCLE project. The data is the average of the TE and TM
measurements which were similar and so the structure is one-dimensional. The
inversion uses 4 layers and should produce a result with a drop in resistivity
at ~35 km, corresponding to the Moho discontinuity. See "The Electric Moho,
Jones & Ferguson (2001)" for a detailed analysis of this data set.

Arguments
=========

- `popsize::Integer`: Population size for the inversion
- `ngen::Integer`: Number of generations for inversion

Returns
=======

Nothing
"""
function testGASnorcle(popsize=150, ngen=3000)
    inDataFile = "./SNO-96-106.data"
    outDataFile = "./SNO-96-106.genetic.model"
    data = readdlm(inDataFile)

    nLayers = 4
    zBounds = [LayerBC(0, 0);
               fill(LayerBC(0, 100_000), (nLayers-1, 1))]
    rBounds = fill(LayerBC(0, 50_000), (nLayers, 1))

    I = Inversion(data, popsize, zBounds, rBounds)
    evolve!(I, ngen, verbose=true)
    print("Best model:\n$(I.pop[1].model)\n")
    p = plotResponseDataCompare(data, I.pop[1].model)
    draw(SVG("inversion-results.svg", 800px, 800px), p)
    return
end

function plotResponseDataCompare(data::Matrix, model::Matrix)
    # Calculate responses
    ρ, Φ = calc_response(model, data[:,1].^-1)

    resMin = data[:,2] - data[:,4]
    resMax = data[:,2] + data[:,4]
    dataResLayer = layer(x=data[:,1],
                         y=data[:,2],
                         ymin=resMin,
                         ymax=resMax,
                         Geom.point,
                         Geom.errorbar)
    phaseMin = data[:,3] - data[:,5]
    phaseMax = data[:,3] + data[:,5]
    dataPhaseLayer = layer(x=data[:,1],
                           y=data[:,3],
                           ymin=phaseMin,
                           ymax=phaseMax,
                           Geom.point,
                           Geom.errorbar)

    modelTheme = Theme(default_color=colorant"red")
    modelResLayer = layer(x=data[:,1],
                          y=ρ,
                          Geom.line,
                          modelTheme)
    modelPhaseLayer = layer(x=data[:,1],
                            y=Φ,
                            Geom.line,
                            modelTheme)

    resPlot = plot(dataResLayer,
                   modelResLayer,
                   Guide.xlabel("Period (S)"),
                   Guide.ylabel("ρa (Ωm)"),
                   Guide.title("Apparent Resistivity"),
                   Scale.x_log10,
                   Scale.y_log10)
    phasePlot = plot(dataPhaseLayer,
                     modelPhaseLayer,
                     Guide.xlabel("Period (S)"),
                     Guide.ylabel("Φ (°)"),
                     Guide.title("Phase"),
                     Scale.x_log10,
                     Scale.y_continuous(minvalue=0, maxvalue=90))

    modelPlot = plot(y=[model[:,2]; model[end,2]],
                     x=[model[:,1] / 1000; model[end,1] / 1000 + 1],
                     Stat.step(direction=:hv),
                     Geom.line,
                     Guide.ylabel("Resistivity (Ωm)"),
                     Guide.xlabel("Depth (km)"),
                     Guide.title("Model"),
                     Scale.y_log10,
                     modelTheme,
                     Guide.manual_color_key("Legend", ["Data", "Model"], ["deepskyblue", "red"]))

    return hstack(vstack(resPlot, phasePlot), modelPlot)
end


text:
struct DyadicRational <: AbstractFloat
    numerator::Int64
    R::Int64

    function DyadicRational(numerator, R)
        if R < 0
            throw(DomainError(R, "Scale must be non-negative"))
        end

        largest_power_of_two_dividing_numerator = trailing_zeros(numerator)
        common_power_of_two = min(R, largest_power_of_two_dividing_numerator)

        new(numerator >> common_power_of_two, R - common_power_of_two)
    end
end

numerator(x::DyadicRational) = x.numerator
resolution(x::DyadicRational) = x.R

function Base.show(io::IO, x::DyadicRational)
    if resolution(x) == 0
        print(io, numerator(x))
    elseif resolution(x) == 1
        print(io, numerator(x), "/2")
    else
        print(io, numerator(x), "/2^", resolution(x))
    end
end

Base.convert(::Type{DyadicRational}, x::Integer) = DyadicRational(x, 0)
Base.promote_rule(::Type{DyadicRational}, ::Type{T}) where {T <: Integer} = DyadicRational

Base.:*(x::DyadicRational, y::DyadicRational) = DyadicRational(numerator(x) * numerator(y), resolution(x) + resolution(y))
Base.:*(x::DyadicRational, k::Integer) = DyadicRational(k * numerator(x), resolution(x))
Base.:*(k::Integer, x::DyadicRational) = x * k

for op in [:+, :-]
    @eval begin
        Base.$op(k::Integer, x::DyadicRational) = $op(x, k)

        function Base.$op(x::DyadicRational, k::Integer)
            new_numerator = $op(numerator(x), k << resolution(x))
            DyadicRational(new_numerator, resolution(x))
        end

        function Base.$op(a::DyadicRational, b::DyadicRational)
            common_numerator = $op(numerator(a) << resolution(b), numerator(b) << resolution(a))
            common_resolution = resolution(a) + resolution(b)

            DyadicRational(common_numerator, common_resolution)
        end
    end
end

for op in [:<, :<=]
    @eval begin
        Base.$op(x::DyadicRational, y::DyadicRational) = $op(numerator(x) << resolution(y), numerator(y) << resolution(x))
    end
end

function Base.Integer(x::DyadicRational)
    if resolution(x) != 0
        throw(InexactError(:Integer, Int64, x))
    end

    numerator(x)
end

function all_dyadic_rationals(left, right, R)
    all_dyadic_rationals(convert(DyadicRational, left), convert(DyadicRational, right), R)
end

function all_dyadic_rationals(left::DyadicRational, right::DyadicRational, R::Integer)
    if R < resolution(left) || R < resolution(right)
        throw(DomainError(R, "Resolution must be at least that of the endpoints"))
    end

    common_resolution = max(resolution(left), resolution(right), R)
    left_numerator = numerator(left) << (common_resolution - resolution(left))
    right_numerator = numerator(right) << (common_resolution - resolution(right))

    DyadicRational.(left_numerator:right_numerator, R)
end

function AbstractFloat(x::DyadicRational)
    numerator(x) / 2.0^resolution(x)
end


text:
# LibFTD2XX.jl - Utility Module
#
# By Reuben Hill 2019, Gowerlabs Ltd, reuben@gowerlabs.co.uk
#
# Copyright (c) Gowerlabs Ltd.

module Util

export ntuple2string, versionnumber

"""
    function ntuple2string(input::NTuple{N, Cchar}) where N

Convert an NTuple of Cchars (optionally null terminated) to a julia string.

# Example

```jldoctest
julia> ntuple2string(Cchar.(('h','e','l','l','o')))
"hello"

julia> ntuple2string(Cchar.(('h','e','l','l','o', '\0', 'x'))) # null terminated
"hello"
```

"""
function ntuple2string(input::NTuple{N, Cchar}) where N
  if any(input .== 0)
    endidx = findall(input .== 0)[1]-1
  elseif all(input .> 0)
    endidx = N
  else
    throw(ArgumentError("No terminator or negative values!"))
  end
  String(UInt8.([char for char in input[1:endidx]]))
end

function versionnumber(hex)
  hex <= 0x999999 || throw(DomainError("Input must be less than 0x999999"))
  patchhex = UInt8( (hex & 0x0000FF) >> 0 )
  patchhex <= 0x99 || throw(DomainError("Patch field must be less than or equal to 0x99"))
  minorhex = UInt8( (hex & 0x00FF00) >> 8 )
  minorhex <= 0x99 || throw(DomainError("Minor field must be less than or equal to 0x99"))
  majorhex = UInt8( (hex & 0xFF0000) >> 16 )
  majorhex <= 0x99 || throw(DomainError("Minor field must be less than or equal to 0x99"))

  patchdec = 10(patchhex >> 4) + (patchhex & 0x0F)
  minordec = 10(minorhex >> 4) + (minorhex & 0x0F)
  majordec = 10(majorhex >> 4) + (majorhex & 0x0F)

  VersionNumber(majordec,minordec,patchdec)
end

end # module Util


text:
using Tune
using Test

@testset "Tune.jl" begin
    # Write your tests here.
end


text:
#!/usr/bin/env julia
# Geof Hannigan

#################
# Load Packages #
#################

# Start by setting up the required packages
using Pkg
Pkg.add("DataFrames")
Pkg.add("Queryverse")
Pkg.add("Gadfly")
using DataFrames
using Queryverse
using Statistics
using Gadfly

#####################
# Load In The Files #
#####################

# First set the file names
metafile = "data/references/genome_id_taxonomy.tsv"
easv_file = "data/processed/rrnDB.easv.count_tibble.tsv"

# Then bring those files in
metadata = DataFrame(load(metafile))
easv = DataFrame(load(easv_file))

# Take a look at the data dataframes (top 10 lines)
first(metadata, 10)
first(easv, 10)

# Merge together
metadata_easv = innerjoin(metadata, easv, on = [:genome_id => :genome])
first(metadata_easv, 10)

outdf = groupby(metadata_easv, [:region, :threshold, :genome_id])

summarydf = combine(outdf, :count => sum, nrow)

# Determine 95%-tile
tileoutdf = groupby(summarydf, [:region, :threshold, :count_sum])
combotile = combine(tileoutdf, :nrow => t -> quantile(t, .95))

# Plot the 95th percentile as a function of region,  threshold, #rrns
plot(combotile, ygroup=:region, x=:count_sum, y=:nrow_function, color=:threshold, Geom.subplot_grid(Geom.line))


text:
# See learn_Plotly.jl for more examples
# http://juliaplots.org/PlotlyJS.jl/stable/

using PlotlyJS, HTTP, CSV
import JSON
import DataFrames
#import RDatasets
import WebIO
#import UUIDs
import DefaultApplication

"""
    dataset(dir, name)

Replace the dataset function from RDatasets.  It stopped working.
"""
function dataset(dir, name)
    url = "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/$(dir)/$(name).csv"
    df = urldownload(url)

    return df
end

function urldownload(url)
    body = HTTP.get(url).body
    csv = CSV.File(body)
    df = DataFrames.DataFrame(csv)

    return df
end

#
# Building Blocks
#

"""

http://juliaplots.org/PlotlyJS.jl/stable/building_traces_layouts/
"""
function sec_Traces()
    # http://juliaplots.org/PlotlyJS.jl/stable/building_traces_layouts/
    t = PlotlyJS.scatter(
        x=[1,2,3,4,5],
        y=[1,6,3,6,1],
        mode="markers+text",
        name="Team A",
        text=["A1", "A2", "A3", "A4", "A5"],
        textposition="top center",
        textfont_family="Raleway, sans-serif", # _ for nested attributes
        marker_size=12
    )
    #print(JSON.json(t,2))

    # Access attributes
    @show getindex(t, :text)
    @show t[:text]
    @show t["textfont.family"]
    @show t["textfont_family"]

    # Setting attributes
    t["marker_color"] ="red"
    t["line_width"] = 5

    return t
end

function sec_Layouts()
    l = Layout(
        title="Penguins",
        xaxis_range=[0, 42.0],
        xaxis_title="fish",
        xaxis_showgrid=true,
        yaxis_title="Weight",
        yaxis_showgrid=true,
        legend_x=0.7,
        legend_y=1.15
    )

    # Using attr()
    l = Layout(
        title="Penguins",
        xaxis=attr(range=[0, 42.0], title="fish", showgrid=true),
        yaxis=attr(title="Weight", showgrid=true),
        legend_x=0.7,
        legend_y=1.15
    )
    return l
end

function sec_DataFrames()
    iris = RDatasets.dataset("datasets", "iris")
    trace = scatter(
        iris,
        x=:SepalLength,
        y=:SepalWidth,
        marker_color=:red
    )

    # TODO: group
    return iris
end

#
# Putting it Together
#

"""

http://juliaplots.org/PlotlyJS.jl/stable/syncplots/#Plot
"""
function sec_Plot()
    iris = dataset("datasets", "iris")
    layout = Layout(title="Iris", width=800, height=600)
    plot = Plot(iris, x=Symbol("Sepal.Length"), y=Symbol("Sepal.Width"), layout,
        mode="markers", marker_size=8,
        group=:Species
    )
    @info typeof(plot)

    return plot
end
#display(sec_Plot())

function sec_SyncPlots1()
    p = sec_Plot()
    sp = plot(p)
    display(sp)

    # Call back
    # only works when launch from console
    obs = WebIO.on(sp["hover"]) do data
        println("\nYou hovered over $data")
    end
end

function sec_SyncPlots2()
    p = plot(rand(10,4))
    display(p)
    WebIO.on(p["hover"]) do data
        println("\nYou hovered over: $data")
    end
end

function sec_SyncPlots3()
    colors = [fill("red", 10), fill("blue", 10)]
    symbols = [fill("circle", 10), fill("circle", 10)]

    ys = [rand(10), rand(10)]
    # list comprehension
    scatters = [
        scatter(y=y, marker=attr(color=c, symbol=s, size=15), line_color=c[1])
        for (y,c,s) in zip(ys, colors, symbols)
    ]
    p = plot(scatters)
    display(p)

    # Callback
    WebIO.on(p["click"]) do data
        colors = (fill("red", 10), fill("blue", 10))
        symbols = (fill("circle", 10), fill("circle", 10))
        for point in data["points"]
            # zero-offset -> unit-offset
            colors[point["curveNumber"] + 1][point["pointIndex"] + 1] = "gold"
            symbols[point["curveNumber"] + 1][point["pointIndex"] + 1] = "star"
        end
        restyle!(p, marker_color=colors, marker_symbol=symbols)
    end

    nothing
end

#
# Styles
#

function sec_Cyclers()
    style = Style(global_trace=attr(
        marker_color=Cycler(["green","red"]),
        line_color="blue"
        )
    )
    #p = plot(rand(10,3), style=style)
    p = plot(
    [
        scatter(y=rand(4)),
        scatter(y=rand(4), marker_color="black"),
        scatter(y=rand(4)),
        scatter(y=rand(4)),
    ],
    style=style
    )
    display(p)
end

#
# 3D
#
using PlotlyJS, DataFrames, Colors, Distributions, LinearAlgebra

function random_line()
    n = 400
    rw() = cumsum(randn(n))
    trace1 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#1f77b4", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="#1f77b4", width=1))
    trace2 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#9467bd", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="rgb(44, 160, 44)", width=1))
    trace3 = scatter3d(;x=rw(),y=rw(), z=rw(), mode="lines",
                        marker=attr(color="#bcbd22", size=12, symbol="circle",
                                    line=attr(color="rgb(0,0,0)", width=0)),
                        line=attr(color="#bcbd22", width=1))
    layout = Layout(autosize=false, width=500, height=500,
                    margin=attr(l=0, r=0, b=0, t=65))
    plot([trace1, trace2, trace3], layout)
end
#display(random_line())

function topo_surface()
    z = Vector[[27.80985, 49.61936, 83.08067, 116.6632, 130.414, 150.7206, 220.1871,
                156.1536, 148.6416, 203.7845, 206.0386, 107.1618, 68.36975, 45.3359,
                49.96142, 21.89279, 17.02552, 11.74317,   14.75226, 13.6671, 5.677561,
                3.31234, 1.156517, -0.147662],
               [27.71966, 48.55022, 65.21374, 95.27666, 116.9964, 133.9056, 152.3412,
                151.934, 160.1139, 179.5327, 147.6184, 170.3943, 121.8194, 52.58537,
                33.08871, 38.40972, 44.24843, 69.5786, 4.019351, 3.050024, 3.039719,
                2.996142, 2.967954, 1.999594],
               [30.4267, 33.47752, 44.80953, 62.47495, 77.43523, 104.2153, 102.7393, 137.0004,
                186.0706, 219.3173, 181.7615, 120.9154, 143.1835, 82.40501, 48.47132,
                74.71461, 60.0909, 7.073525, 6.089851, 6.53745, 6.666096, 7.306965, 5.73684,
                3.625628],
               [16.66549, 30.1086, 39.96952, 44.12225, 59.57512, 77.56929, 106.8925,
                166.5539, 175.2381, 185.2815, 154.5056, 83.0433, 62.61732, 62.33167,
                60.55916, 55.92124, 15.17284, 8.248324, 36.68087, 61.93413, 20.26867,
                68.58819, 46.49812, 0.2360095],
               [8.815617, 18.3516, 8.658275, 27.5859, 48.62691, 60.18013, 91.3286,
                145.7109, 116.0653, 106.2662, 68.69447, 53.10596, 37.92797, 47.95942,
                47.42691, 69.20731, 44.95468, 29.17197, 17.91674, 16.25515, 14.65559,
                17.26048, 31.22245, 46.71704],
               [6.628881, 10.41339, 24.81939, 26.08952, 30.1605, 52.30802, 64.71007,
                76.30823, 84.63686, 99.4324, 62.52132, 46.81647, 55.76606, 82.4099,
                140.2647, 81.26501, 56.45756, 30.42164, 17.28782, 8.302431, 2.981626,
                2.698536, 5.886086, 5.268358],
               [21.83975, 6.63927, 18.97085, 32.89204, 43.15014, 62.86014, 104.6657,
                130.2294, 114.8494, 106.9873, 61.89647, 55.55682, 86.80986, 89.27802,
                122.4221, 123.9698, 109.0952, 98.41956, 77.61374, 32.49031, 14.67344,
                7.370775, 0.03711011, 0.6423392],
               [53.34303, 26.79797, 6.63927, 10.88787, 17.2044, 56.18116, 79.70141,
                90.8453, 98.27675, 80.87243, 74.7931, 75.54661, 73.4373, 74.11694, 68.1749,
                46.24076, 39.93857, 31.21653, 36.88335, 40.02525, 117.4297, 12.70328,
                1.729771, 0],
               [25.66785, 63.05717, 22.1414, 17.074, 41.74483, 60.27227, 81.42432, 114.444,
                102.3234, 101.7878, 111.031, 119.2309, 114.0777, 110.5296, 59.19355,
                42.47175, 14.63598, 6.944074, 6.944075, 27.74936, 0, 0, 0.09449376, 0.07732264],
               [12.827, 69.20554, 46.76293, 13.96517, 33.88744, 61.82613, 84.74799,
                121.122, 145.2741, 153.1797, 204.786, 227.9242, 236.3038, 228.3655,
                79.34425, 25.93483, 6.944074, 6.944074, 6.944075, 7.553681, 0, 0, 0, 0],
               [0, 68.66396, 59.0435, 33.35762, 47.45282, 57.8355, 78.91689, 107.8275,
                168.0053, 130.9597, 212.5541, 165.8122, 210.2429, 181.1713, 189.7617,
                137.3378, 84.65395, 8.677168, 6.956576, 8.468093, 0, 0, 0, 0],
               [0, 95.17499, 80.03818, 59.89862, 39.58476, 50.28058, 63.81641, 80.61302,
                66.37824, 198.7651, 244.3467, 294.2474, 264.3517, 176.4082, 60.21857,
                77.41475, 53.16981, 56.16393, 6.949235, 7.531059, 3.780177, 0, 0, 0],
               [0, 134.9879, 130.3696, 96.86325, 75.70494, 58.86466, 57.20374, 55.18837,
                78.128, 108.5582, 154.3774, 319.1686, 372.8826, 275.4655, 130.2632, 54.93822,
                25.49719, 8.047439, 8.084393, 5.115252, 5.678269, 0, 0, 0],
               [0, 48.08919, 142.5558, 140.3777, 154.7261, 87.9361, 58.11092, 52.83869,
                67.14822, 83.66798, 118.9242, 150.0681, 272.9709, 341.1366, 238.664, 190.2,
                116.8943, 91.48672, 14.0157, 42.29277, 5.115252, 0, 0, 0],
               [0, 54.1941, 146.3839, 99.48143, 96.19411, 102.9473, 76.14089, 57.7844,
                47.0402, 64.36799, 84.23767, 162.7181, 121.3275, 213.1646, 328.482,
                285.4489, 283.8319, 212.815, 164.549, 92.29631, 7.244015, 1.167, 0, 0],
               [0, 6.919659, 195.1709, 132.5253, 135.2341, 89.85069, 89.45549, 60.29967,
                50.33806, 39.17583, 59.06854, 74.52159, 84.93402, 187.1219, 123.9673,
                103.7027, 128.986, 165.1283, 249.7054, 95.39966, 10.00284, 2.39255, 0, 0],
               [0, 21.73871, 123.1339, 176.7414, 158.2698, 137.235, 105.3089, 86.63255, 53.11591,
                29.03865, 30.40539, 39.04902, 49.23405, 63.27853, 111.4215, 101.1956,
                40.00962, 59.84565, 74.51253, 17.06316, 2.435141, 2.287471, -0.0003636982, 0],
               [0, 0, 62.04672, 136.3122, 201.7952, 168.1343, 95.2046, 58.90624, 46.94091,
                49.27053, 37.10416, 17.97011, 30.93697, 33.39257, 44.03077, 55.64542,
                78.22423, 14.42782, 9.954997, 7.768213, 13.0254, 21.73166, 2.156372,
                0.5317867],
               [0, 0, 79.62993, 139.6978, 173.167, 192.8718, 196.3499, 144.6611, 106.5424,
                57.16653, 41.16107, 32.12764, 13.8566, 10.91772, 12.07177, 22.38254,
                24.72105, 6.803666, 4.200841, 16.46857, 15.70744, 33.96221, 7.575688,
                -0.04880907],
               [0, 0, 33.2664, 57.53643, 167.2241, 196.4833, 194.7966, 182.1884, 119.6961,
                73.02113, 48.36549, 33.74652, 26.2379, 16.3578, 6.811293, 6.63927, 6.639271,
                8.468093, 6.194273, 3.591233, 3.81486, 8.600739, 5.21889, 0],
               [0, 0, 29.77937, 54.97282, 144.7995, 207.4904, 165.3432, 171.4047, 174.9216,
                100.2733, 61.46441, 50.19171, 26.08209, 17.18218, 8.468093, 6.63927,
                6.334467, 6.334467, 5.666687, 4.272203, 0, 0, 0, 0],
               [0, 0, 31.409, 132.7418, 185.5796, 121.8299, 185.3841, 160.6566, 116.1478,
                118.1078, 141.7946, 65.56351, 48.84066, 23.13864, 18.12932, 10.28531,
                6.029663, 6.044627, 5.694764, 3.739085, 3.896037, 0, 0, 0],
               [0, 0, 19.58994, 42.30355, 96.26777, 187.1207, 179.6626, 221.3898, 154.2617,
                142.1604, 148.5737, 67.17937, 40.69044, 39.74512, 26.10166, 14.48469,
                8.65873, 3.896037, 3.571392, 3.896037, 3.896037, 3.896037, 1.077756, 0],
               [0.001229679, 3.008948, 5.909858, 33.50574, 104.3341, 152.2165, 198.1988,
                191.841, 228.7349, 168.1041, 144.2759, 110.7436, 57.65214, 42.63504,
                27.91891, 15.41052, 8.056102, 3.90283, 3.879774, 3.936718, 3.968634,
                0.1236256, 3.985531, -0.1835741],
               [0, 5.626141, 7.676256, 63.16226, 45.99762, 79.56688, 227.311, 203.9287,
                172.5618, 177.1462, 140.4554, 123.9905, 110.346, 65.12319, 34.31887,
                24.5278, 9.561069, 3.334991, 5.590495, 5.487353, 5.909499, 5.868994,
                5.833817, 3.568177]]
    trace = surface(z=z)
    layout = Layout(title="Mt. Bruno Elevation", autosize=false, width=500,
                    height=500, margin=attr(l=65, r=50, b=65, t=90))
    plot(trace, layout)
end
#display(topo_surface())
function multiple_surface()
    z1 = Vector[[8.83, 8.89, 8.81, 8.87, 8.9, 8.87],
                [8.89, 8.94, 8.85, 8.94, 8.96, 8.92],
                [8.84, 8.9, 8.82, 8.92, 8.93, 8.91],
                [8.79, 8.85, 8.79, 8.9, 8.94, 8.92],
                [8.79, 8.88, 8.81, 8.9, 8.95, 8.92],
                [8.8, 8.82, 8.78, 8.91, 8.94, 8.92],
                [8.75, 8.78, 8.77, 8.91, 8.95, 8.92],
                [8.8, 8.8, 8.77, 8.91, 8.95, 8.94],
                [8.74, 8.81, 8.76, 8.93, 8.98, 8.99],
                [8.89, 8.99, 8.92, 9.1, 9.13, 9.11],
                [8.97, 8.97, 8.91, 9.09, 9.11, 9.11],
                [9.04, 9.08, 9.05, 9.25, 9.28, 9.27],
                [9, 9.01, 9, 9.2, 9.23, 9.2],
                [8.99, 8.99, 8.98, 9.18, 9.2, 9.19],
                [8.93, 8.97, 8.97, 9.18, 9.2, 9.18]]
    z2 = map(x->x.+1, z1)
    z3 = map(x->x.-1, z1)
    trace1 = surface(z=z1, colorscale="Viridis")
    trace2 = surface(z=z2, showscale=false, opacity=0.9, colorscale="Viridis")
    trace3 = surface(z=z3, showscale=false, opacity=0.9, colorscale="Viridis")
    plot([trace1, trace2, trace3])
end
#display(multiple_surface())

function clustering_alpha_shapes()
    # load data
    iris = dataset("datasets", "iris")
    nms = unique(iris[!,:Species])
    colors = [RGB(0.89, 0.1, 0.1), RGB(0.21, 0.50, 0.72), RGB(0.28, 0.68, 0.3)]

    data = GenericTrace[]

    for (i, nm) in enumerate(nms)
        df = iris[iris[!,:Species] .== nm, :]
        x=df[!,Symbol("Sepal.Length")]
        y=df[!,Symbol("Sepal.Width")]
        z=df[!,Symbol("Petal.Length")]
        trace = scatter3d(;name=nm, mode="markers",
                           marker_size=3, marker_color=colors[i], marker_line_width=0,
                           x=x, y=y, z=z)
        push!(data, trace)

        cluster = mesh3d(;color=colors[i], opacity=0.3, x=x, y=y, z=z)
        push!(data, cluster)
    end

    # notice the nested attrs to create complex JSON objects
    layout = Layout(width=800, height=550, autosize=false, title="Iris dataset",
                    scene=attr(xaxis=attr(gridcolor="rgb(255, 255, 255)",
                                          zerolinecolor="rgb(255, 255, 255)",
                                          showbackground=true,
                                          backgroundcolor="rgb(230, 230,230)"),
                               yaxis=attr(gridcolor="rgb(255, 255, 255)",
                                           zerolinecolor="rgb(255, 255, 255)",
                                           showbackground=true,
                                           backgroundcolor="rgb(230, 230,230)"),
                               zaxis=attr(gridcolor="rgb(255, 255, 255)",
                                           zerolinecolor="rgb(255, 255, 255)",
                                           showbackground=true,
                                           backgroundcolor="rgb(230, 230,230)"),
                               aspectratio=attr(x=1, y=1, z=0.7),
                               aspectmode = "manual"))

    # Need to use Plot to work with Dash
    # return plot(data, layout)
    return Plot(data, layout)
end
#display(clustering_alpha_shapes())

function scatter_3d()
    Σ = fill(0.5, 3, 3) + Diagonal([0.5, 0.5, 0.5])
    obs1 = rand(MvNormal(zeros(3), Σ), 200)'
    obs2 = rand(MvNormal(zeros(3), 0.5Σ), 100)'

    trace1 = scatter3d(;x=obs1[:, 1], y=obs1[:, 2], z=obs1[:, 3],
                        mode="markers", opacity=0.8,
                        marker_size=12, marker_line_width=0.5,
                        marker_line_color="rgba(217, 217, 217, 0.14)")

    trace2 = scatter3d(;x=obs2[:, 1], y=obs2[:, 2], z=obs2[:, 3],
                        mode="markers", opacity=0.9,
                        marker=attr(color="rgb(127, 127, 127)",
                                    symbol="circle", line_width=1.0,
                                    line_color="rgb(204, 204, 204)"))

    layout = Layout(margin=attr(l=0, r=0, t=0, b=0))

    return Plot([trace1, trace2], layout)
end
#display(scatter_3d())

function trisurf()
    facecolor = repeat([
    	"rgb(50, 200, 200)",
    	"rgb(100, 200, 255)",
    	"rgb(150, 200, 115)",
    	"rgb(200, 200, 50)",
    	"rgb(230, 200, 10)",
    	"rgb(255, 140, 0)"
    ], inner=[2])

    t = mesh3d(
        x=[0, 0, 1, 1, 0, 0, 1, 1],
        y=[0, 1, 1, 0, 0, 1, 1, 0],
        z=[0, 0, 0, 0, 1, 1, 1, 1],
        i=[7, 0, 0, 0, 4, 4, 2, 6, 4, 0, 3, 7],
        j=[3, 4, 1, 2, 5, 6, 5, 5, 0, 1, 2, 2],
        k=[0, 7, 2, 3, 6, 7, 1, 2, 5, 5, 7, 6],
        facecolor=facecolor)

    plot(t)
end
#display(trisurf())

function meshcube()
    t = mesh3d(
        x=[0, 0, 1, 1, 0, 0, 1, 1],
        y=[0, 1, 1, 0, 0, 1, 1, 0],
        z=[0, 0, 0, 0, 1, 1, 1, 1],
        i=[7, 0, 0, 0, 4, 4, 6, 6, 4, 0, 3, 2],
        j=[3, 4, 1, 2, 5, 6, 5, 2, 0, 1, 6, 3],
        k=[0, 7, 2, 3, 6, 7, 1, 1, 5, 5, 7, 6],
        intensity=range(0, stop=1, length=8),
        colorscale=[
            [0, "rgb(255, 0, 255)"],
            [0.5, "rgb(0, 255, 0)"],
            [1, "rgb(0, 0, 255)"]
        ]
    )
    plot(t)
end
#display(meshcube())

#
# Contour
#
function contour1()

end

#
# Line and Scatter
#
using PlotlyJS, DataFrames, CSV, Dates
function linescatter1()
    l = Layout(
        #title="Penguins",
        #xaxis=attr(range=[0, 42.0], title="fish", showgrid=true),
        #yaxis=attr(title="Weight", showgrid=true),
        legend=attr(x=0.1, y=0.1)
    )
    trace1 = scatter(;x=1:4, y=[10, 15, 13, 17], mode="markers")
    trace2 = scatter(;x=2:5, y=[16, 5, 11, 9], mode="lines")
    trace3 = scatter(;x=1:4, y=[12, 9, 15, 12], mode="lines+markers")
    #p = Plot([trace1, trace2, trace3], l)
    plot([trace1, trace2, trace3], l)
end
#display(linescatter1())

function linescatter2()
    trace1 = scatter(;x=1:5, y=[1, 6, 3, 6, 1],
                      mode="markers", name="Team A",
                      text=["A-1", "A-2", "A-3", "A-4", "A-5"],
                      marker_size=12)

    trace2 = scatter(;x=1:5+0.5, y=[4, 1, 7, 1, 4],
                      mode="markers", name= "Team B",
                      text=["B-a", "B-b", "B-c", "B-d", "B-e"])
    # setting marker.size this way is _equivalent_ to what we did for trace1
    trace2["marker"] = Dict(:size => 12)

    data = [trace1, trace2]
    layout = Layout(;title="Data Labels Hover", xaxis_range=[0.75, 5.25],
                     yaxis_range=[0, 8])
    plot(data, layout)
end
#display(linescatter2())

function linescatter3()
    trace1 = scatter(;x=1:5, y=[1, 6, 3, 6, 1],
                      mode="markers+text", name="Team A",
                      textposition="top center",
                      text=["A-1", "A-2", "A-3", "A-4", "A-5"],
                      marker_size=12, textfont_family="Raleway, sans-serif")

    trace2 = scatter(;x=1:5+0.5, y=[4, 1, 7, 1, 4],
                      mode="markers+text", name= "Team B",
                      textposition="bottom center",
                      text= ["B-a", "B-b", "B-c", "B-d", "B-e"],
                      marker_size=12, textfont_family="Times New Roman")

    data = [trace1, trace2]

    layout = Layout(;title="Data Labels on the Plot", xaxis_range=[0.75, 5.25],
                     yaxis_range=[0, 8], legend_y=0.5, legend_yref="paper",
                     legend=attr(family="Arial, sans-serif", size=20,
                                 color="grey"))
    plot(data, layout)
end
#display(linescatter3())

function linescatter4()
    trace1 = scatter(;y=fill(5, 40), mode="markers", marker_size=40,
                      marker_color=0:39)
    layout = Layout(title="Scatter Plot with a Color Dimension")
    plot(trace1, layout)
end
#display(linescatter4())

function linescatter5()

    country = ["Switzerland (2011)", "Chile (2013)", "Japan (2014)",
               "United States (2012)", "Slovenia (2014)", "Canada (2011)",
               "Poland (2010)", "Estonia (2015)", "Luxembourg (2013)",
               "Portugal (2011)"]

    votingPop = [40, 45.7, 52, 53.6, 54.1, 54.2, 54.5, 54.7, 55.1, 56.6]
    regVoters = [49.1, 42, 52.7, 84.3, 51.7, 61.1, 55.3, 64.2, 91.1, 58.9]

    # notice use of `attr` function to make nested attributes
    trace1 = scatter(;x=votingPop, y=country, mode="markers",
                      name="Percent of estimated voting age population",
                      marker=attr(color="rgba(156, 165, 196, 0.95)",
                                  line_color="rgba(156, 165, 196, 1.0)",
                                  line_width=1, size=16, symbol="circle"))

    trace2 = scatter(;x=regVoters, y=country, mode="markers",
                      name="Percent of estimated registered voters")
    # also could have set the marker props above by using a dict
    trace2["marker"] = Dict(:color => "rgba(204, 204, 204, 0.95)",
                           :line => Dict(:color=> "rgba(217, 217, 217, 1.0)",
                                         :width=> 1),
                           :symbol => "circle",
                           :size => 16)

    data = [trace1, trace2]
    layout = Layout(Dict{Symbol,Any}(:paper_bgcolor => "rgb(254, 247, 234)",
                                     :plot_bgcolor => "rgb(254, 247, 234)");
                    title="Votes cast for ten lowest voting age population in OECD countries",
                    width=600, height=600, hovermode="closest",
                    margin=Dict(:l => 140, :r => 40, :b => 50, :t => 80),
                    xaxis=attr(showgrid=false,
                               showline=true,
                               linecolor="rgb(102, 102, 102)",
                               titlefont_color="rgb(204, 204, 204)",
                               tickfont_color="rgb(102, 102, 102)",
                               autotick=false, dtick=10, ticks="outside",
                               tickcolor="rgb(102, 102, 102)"),
                    legend=attr(font_size=10, yanchor="middle",
                                xanchor="right"),
    )

    return Plot(data, layout)
end
#display(linescatter5())

function linescatter6()
    trace1 = scatter(;x=[52698, 43117], y=[53, 31],
                      mode="markers",
                      name="North America",
                      text=["United States", "Canada"],
                      marker=attr(color="rgb(164, 194, 244)", size=12,
                                  line=attr(color="white", width=0.5))
                      )

    trace2 = scatter(;x=[39317, 37236, 35650, 30066, 29570, 27159, 23557, 21046, 18007],
                      y=[33, 20, 13, 19, 27, 19, 49, 44, 38],
                      mode="markers", name="Europe",
                      marker_size=12, marker_color="rgb(255, 217, 102)",
                      text=["Germany", "Britain", "France", "Spain", "Italy",
                            "Czech Rep.", "Greece", "Poland", "Portugal"])

    trace3 = scatter(;x=[42952, 37037, 33106, 17478, 9813, 5253, 4692, 3899],
                      y=[23, 42, 54, 89, 14, 99, 93, 70],
                      mode="markers",
                      name="Asia/Pacific",
                      marker_size=12, marker_color="rgb(234, 153, 153)",
                      text=["Australia", "Japan", "South Korea", "Malaysia",
                            "China", "Indonesia", "Philippines", "India"])

    trace4 = scatter(;x=[19097, 18601, 15595, 13546, 12026, 7434, 5419],
                      y=[43, 47, 56, 80, 86, 93, 80],
                      mode="markers", name="Latin America",
                      marker_size=12, marker_color="rgb(142, 124, 195)",
                      text=["Chile", "Argentina", "Mexico", "Venezuela",
                            "Venezuela", "El Salvador", "Bolivia"])

    data = [trace1, trace2, trace3, trace4]

    layout = Layout(;title="Quarter 1 Growth",
                     xaxis=attr(title="GDP per Capital", showgrid=false, zeroline=false),
                     yaxis=attr(title="Percent", zeroline=false))

    plot(data, layout)
end
#display(linescatter6())

function dumbell()
    # reference: https://plot.ly/r/dumbbell-plots/
    # read Data into dataframe
    nm = tempname()
    url = "https://raw.githubusercontent.com/plotly/datasets/master/school_earnings.csv"
    download(url, nm)
    df = CSV.read(nm, DataFrame)
    rm(nm)

    # sort dataframe by male earnings
    df = sort(df, :Men, rev=false)

    men = scatter(;y=df[!,:School], x=df[!,:Men], mode="markers", name="Men",
                   marker=attr(color="blue", size=12))
    women = scatter(;y=df[!,:School], x=df[!,:Women], mode="markers", name="Women",
                     marker=attr(color="pink", size=12))

    lines = map(eachrow(df)) do r
        scatter(y=fill(r[:School], 2), x=[r[:Women], r[:Men]], mode="lines",
                name=r[:School], showlegend=false, line_color="gray")
    end

    data = Base.typed_vcat(GenericTrace, men, women, lines)
    layout = Layout(width=650, height=650, margin_l=100, yaxis_title="School",
                    xaxis_title="Annual Salary (thousands)",
                    title="Gender earnings disparity")

    plot(data, layout)
end
#display(dumbell())

function errorbars1()
    trace1 = scatter(;x=vcat(1:10, 10:-1:1),
                     y=vcat(2:11, 9:-1:0),
                     fill="tozerox",
                     fillcolor="rgba(0, 100, 80, 0.2)",
                     line_color="transparent",
                     name="Fair",
                     showlegend=false)

    trace2 = scatter(;x=vcat(1:10, 10:-1:1),
                     y=[5.5, 3.0, 5.5, 8.0, 6.0, 3.0, 8.0, 5.0, 6.0, 5.5, 4.75,
                        5.0, 4.0, 7.0, 2.0, 4.0, 7.0, 4.4, 2.0, 4.5],
                     fill="tozerox",
                     fillcolor="rgba(0, 176, 246, 0.2)",
                     line_color="transparent",
                     name="Premium",
                     showlegend=false)

    trace3 = scatter(;x=vcat(1:10, 10:-1:1),
                     y=[11.0, 9.0, 7.0, 5.0, 3.0, 1.0, 3.0, 5.0, 3.0, 1.0,
                        -1.0, 1.0, 3.0, 1.0, -0.5, 1.0, 3.0, 5.0, 7.0, 9.],
                     fill="tozerox",
                     fillcolor="rgba(231, 107, 243, 0.2)",
                     line_color="transparent",
                     name="Fair",
                     showlegend=false)

    trace4 = scatter(;x=1:10, y=1:10,
                     line_color="rgb(00, 100, 80)",
                     mode="lines",
                     name="Fair")

    trace5 = scatter(;x=1:10,
                     y=[5.0, 2.5, 5.0, 7.5, 5.0, 2.5, 7.5, 4.5, 5.5, 5.],
                     line_color="rgb(0, 176, 246)",
                     mode="lines",
                     name="Premium")

    trace6 = scatter(;x=1:10, y=vcat(10:-2:0, [2, 4,2, 0]),
                     line_color="rgb(231, 107, 243)",
                     mode="lines",
                     name="Ideal")
    data = [trace1, trace2, trace3, trace4, trace5, trace6]
    layout = Layout(;paper_bgcolor="rgb(255, 255, 255)",
                    plot_bgcolor="rgb(229, 229, 229)",

                    xaxis=attr(gridcolor="rgb(255, 255, 255)",
                               range=[1, 10],
                               showgrid=true,
                               showline=false,
                               showticklabels=true,
                               tickcolor="rgb(127, 127, 127)",
                               ticks="outside",
                               zeroline=false),

                    yaxis=attr(gridcolor="rgb(255, 255, 255)",
                               showgrid=true,
                               showline=false,
                               showticklabels=true,
                               tickcolor="rgb(127, 127, 127)",
                               ticks="outside",
                               zeroline=false))

    plot(data, layout)
end
#display(errorbars1())

#t = sec_Traces()
#l = sec_Layouts()
#iris = sec_DataFrames()
#p = sec_Plot()
#display(plot(p))

# Open in Pane when run from VS Code
# Open in Electron when run from console
#sec_SyncPlots()

# Open in default browser
#savefig(p, "plotlyjs.html")
#savefig(p, "plotlyjs.png")
#DefaultApplication.open("plotlyjs.html")

#sec_Cyclers()

nothing


text:
TikzDocument(elements) = PGFPlotsX.TikzDocument("\\sffamily\\sansmath",elements;use_default_preamble = false, preamble = CUSTOM_PREAMBLE);

function plot(coefplot::Coefplot, filename::String; verbose::Bool=false)
    td = TikzDocument(TikzPicture(to_plotable(coefplot)))
    verbose ? print_tex(td) : nothing
    pgfsave(filename,td)
end

function plot(coefplot::Coefplot; verbose::Bool=false)
    td = TikzDocument(TikzPicture(to_plotable(coefplot)))
    if verbose
        @show td
        print_tex(td)
    else
        display_on_different_IDEs(td)
    end
end

function plot(mcoefplot::MultiCoefplot; offset::Union{Real,Missing}=missing, verbose::Bool=false)
    td = TikzDocument(TikzPicture(to_plotable(mcoefplot, offset)))
    if verbose
        @show td
        print_tex(td)
    else
        display_on_different_IDEs(td)
    end
end

function plot(mcoefplot::MultiCoefplot, filename::String; offset::Union{Real,Missing}=missing, verbose::Bool=false)
    td = TikzDocument(TikzPicture(to_plotable(mcoefplot, offset)))
    verbose ? print_tex(td) : nothing
    pgfsave(filename,td)
end

function plot(regmodel::SupportedEstimation, filename::String; verbose::Bool=false)
    parsed_model = parse(regmodel)
    plot(parsed_model,filename; verbose)
end

function plot(regmodel::SupportedEstimation; verbose::Bool=false)
    parsed_model = parse(regmodel)
    plot(parsed_model; verbose)
end

function display_on_different_IDEs(td::PGFPlotsX.TikzDocument)
    if PGFPlotsX._is_ijulia()
        display("image/png",td)
    else
        display(td)
    end
end


text:
#=

Plots reservation wage against the discount factor

=#

using Plots, LaTeXStrings
pyplot()

include("mccall_bellman_iteration.jl")
include("compute_reservation_wage.jl")

grid_size = 25
beta_vals = linspace(0.8, 0.99, grid_size)
w_bar_vals = similar(beta_vals)

mcm = McCallModel()

for (i, beta) in enumerate(beta_vals)
    mcm.beta = beta
    w_bar = compute_reservation_wage(mcm)
    w_bar_vals[i] = w_bar
end

plot(beta_vals,
    w_bar_vals,
    lw=2,
    alpha=0.7,
    xlabel="discount rate",
    ylabel="reservation wage",
    label=L"$\bar w$ as a function of $\beta$")

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
module ArcadeLearningEnvironment
using ArcadeLearningEnvironment_jll
using Pkg.Artifacts

include("aleinterface.jl")

ROM_PATH = artifact"atari_roms"

export ALEInterface,
    ALEPtr,
    # Functions
    ALE_new,
    ALE_del,
    getInt,
    getBool,
    getFloat,
    setString,
    setInt,
    setBool,
    setFloat,
    loadROM,
    act,
    game_over,
    reset_game,
    getLegalActionSet,
    getLegalActionSize,
    getMinimalActionSet,
    getMinimalActionSize,
    getFrameNumber,
    lives,
    getEpisodeFrameNumber,
    getScreen,
    getRAM,
    getRAMSize,
    getScreenWidth,
    getScreenHeight,
    getScreenRGB,
    getScreenGrayscale,
    saveState,
    loadState,
    cloneState,
    restoreState,
    cloneSystemState,
    restoreSystemState,
    deleteState,
    saveScreenPNG,
    encodeState,
    encodeStateLen,
    decodeState,
    setLoggerMode!,
    getROMList

end


text:
using SSDPClient

# const ARDUINOMATCH = "SERVER: Arduino.+(http://[^/]+)/description"
# println(ssdpquery(ARDUINOMATCH)[1])

# println(ssdpquery("Roku"))

println(ssdpquery())

true


text:
using Documenter, Wright

makedocs(
    modules = [Wright],
    sitename = "Wright.jl",
    format = Documenter.HTML(
        assets = ["assets/favicon.ico"],
    ),

    pages = [
        "index.md",
    ]
)

deploydocs(
    repo = "github.com/SciFracX/Wright.jl.git";
)


text:
export set_lagrange_term!, set_mayer_term!

function set_lagrange_term!( model::Model, lagrange::Union{Nothing,Function}; dx::CallbackArg=:autodiff, du::CallbackArg=:autodiff, dp::CallbackArg=:autodiff )
    # Allow the lagrange term to be disabled
    if isa( lagrange, Nothing )
        model.integralCost = false

        model.callbacks.lfct_cb = nothing
        model.callbacks.dldx_cb = nothing
        model.callbacks.dldu_cb = nothing
        model.callbacks.dldp_cb = nothing

        return
    else
        model.integralCost = true
        model.callbacks.lfct_cb = lagrange
    end

    if isa( dx, Function )
        model.callbacks.dldx_cb = dx
    elseif dx == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( du, Function )
        model.callbacks.dldu_cb = du
    elseif du == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(du)) for du" ) )
    end

    if isa( dp, Function )
        model.callbacks.dldp_cb = dp
    elseif dp == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end
end

function set_mayer_term!( model::Model, mayer::Union{Nothing,Function}; dx::CallbackArg=:autodiff, dp::CallbackArg=:autodiff, dT::CallbackArg=:autodiff )
    # Allow the mayer term to be disabled
    if isa( mayer, Nothing )
        model.terminalCost = false

        model.callbacks.Vfct_cb = nothing
        model.callbacks.dVdx_cb = nothing
        model.callbacks.dVdp_cb = nothing
        model.callbacks.dVdT_cb = nothing

        return
    else
        model.terminalCost = true
        model.callbacks.Vfct_cb = mayer
    end

    if isa( dx, Function )
        model.callbacks.dVdx_cb = dx
    elseif dx == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( dp, Function )
        model.callbacks.dVdp_cb = dp
    elseif dp == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end

    if isa( dT, Function )
        model.callbacks.dVdT_cb = dT
    elseif dT == :autodiff
        # TODO
    else
        throw( ArgumentError( "Unknown option $(String(dT)) for dT" ) )
    end
end


text:
# TODO special methods for array types

import Distributions: params, ncategories, probs
import Distributions: rand, pdf, logpdf, cdf
import Distributions: insupport, mode, modes
import Distributions: mgf, cf
import Base: mean, var, minimum, maximum, std
import StatsBase: entropy, moment, skewness, kurtosis

# for now all these methods just call from underlying categorical
params(d::SampleDistribution) = params(d.dist)
ncategories(d::SampleDistribution) = ncategories(d.dist)
probs(d::SampleDistribution) = probs(d.dist)

rand(d::SampleDistribution) = getsample(d, rand(d.dist))
pdf(d::SampleDistribution{T,K}, x::K) where {T,K} = pdf(d.dist, d.map[x])
logpdf(d::SampleDistribution{T,K}, x::K) where {T,K} = log(pdf(d, x))

insupport(d::SampleDistribution{T,K}, x::K) where {T,K} = (x ∈ supportkeys(d))

# requires ≥
function cdf(d::SampleDistribution{T,K}, x::K) where {T,K}
    sum(pdf(d, k) for k ∈ supportkeys(d) if k ≥ x)
end

# it would be up to users to define quantile

# requires * and +
mean(d::SampleDistribution) = expectation(d, x -> x)
# requires ≤
minimum(d::SampleDistribution) = minimum(supportkeys(d))
# requlres ≥
maximum(d::SampleDistribution) = maximum(supportkeys(d))
mode(d::SampleDistribution) = getsample(d, mode(d.dist))
modes(d::SampleDistribution) = getsample.(d, modes(d.dist))

# written this way so users only need to define - and *
var(d::SampleDistribution{T,K}, μ=mean(d)) where {T,K} = expectation(d, x -> (x - μ)*(x - μ))

std(d::SampleDistribution{T,K}, μ=mean(d)) where {T,K} = sqrt(var(d, μ))

# TODO these probably really inefficient, why mismatch with StatsBase?
moment(d::SampleDistribution, k::Integer, μ=mean(d)) = expectation(d, x -> (x - μ)^k)
skewness(d::SampleDistribution, μ=mean(d)) = moment(d, 3, μ)/std(d, μ)^3
kurtosis(d::SampleDistribution, μ=mean(d)) = moment(d, 4, μ)/var(d, μ)^2

entropy(d::SampleDistribution) = entropy(d.dist)

# requires exp(t*k)
mgf(d::SampleDistribution, t::Real) = expectation(d, x -> exp(t*x))
# requires exp(im*t*k)
cf(d::SampleDistribution, t::Real) = expectation(d, x -> exp(im*t*x))


text:
using BinaryBuilder

# Collection of sources required to build Nettle
name = "dmg2img"
version = v"1.6.7"
sources = [
    GitSource("https://github.com/Lekensteyn/dmg2img.git", "f16f247d30f868e84f31e24792b4464488f1c009"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/dmg2img*/

make -j${nproc} CFLAGS="-O2 -Wall -I${prefix}/include" LDFLAGS="-L${prefix}/lib -lssl"
make install DESTDIR=${prefix}
mv ${prefix}/usr/bin/* ${prefix}/bin/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(!Sys.isfreebsd, platforms)
filter!(!Sys.iswindows, platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
    Dependency("Bzip2_jll"),
    Dependency("OpenSSL_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# ASuMa Sept, 2018
# This script takes an AdaGram model and evalautes 4 similarity measures
# using the SCWS dataset (Huang et al. 2012)

# usage:
#   evaluate_SCWS.sh model(file or dir)
#                   output-dir [--joiner=@] [--min-prob=0.3] [--window=4]
# see arg table for meaning of parameters

#push!(LOAD_PATH, "./src/")

function disambiguate_word(vm, dict, word, valid_context, window)
"""
    Disambiguates a valid word using valid context
"""
        return disambiguate(vm, dict, word, valid_context)
end

function similarity(vm, dict, w1, s1, w2, s2)
"""
    Calculates cosine distance between two sense vectors
"""
    v1 = dict.word2id[w1]
    v2 = dict.word2id[w2]
    sense_vec1 = vec(vm, v1, s1)
    sense_vec2 = vec(vm, v2, s2)
    return dot(sense_vec1, sense_vec2)
end

function evaluate_entry(vm, dict, word1, word2, context1, context2; min_prob=1e-2)
"""
    Evaluates 4 similarity measures for a pair or words
"""
    priors1 = expected_pi(vm, dict.word2id[word1])
    priors2 = expected_pi(vm, dict.word2id[word2])
    posteriors1 = disambiguate(vm, dict, word1, context1)
    posteriors2 = disambiguate(vm, dict, word2, context2)

    # AvgSim & MaxSim (no context involved)
    # AvgSimC & MaxSimC (context taken into account)
    similarities = []
    similaritiesC = []
    for s1 in 1:T(vm)
        if priors1[s1] > min_prob
            for s2 in 1:T(vm)
                if priors2[s2] > min_prob
                    sim = similarity(vm, dict, word1, s1, word2, s2)
                    push!(similarities, sim)
                    push!(similaritiesC, sim * posteriors1[s1] * posteriors2[s2])
                end
            end
        end
    end
    AvgSim = mean(similarities)
    MaxSim = maximum(similarities)
    AvgSimC = sum(similaritiesC)

    best_sense1 = findmax(posteriors1)[2]
    best_sense2 = findmax(posteriors2)[2]
    MaxSimC = similarity(vm, dict, word1, best_sense1, word2, best_sense2)

    return AvgSim, MaxSim, AvgSimC, MaxSimC
end

function words_in_model(dict, word1, word2)
"""
    Checks if both words in pair are included in vector model
"""
    if (word1 in dict.id2word) && (word2 in dict.id2word)
        return true
    else
        return false
    end
end

function load_SCWS(path)
"""
    Loads SCWS database in path into a matrix where each row is of the form:
    [<id> <word1> <POS of word1> <word2> <POS of word2> <word1 in context> <word2 in context> <average human rating> <10 individual human ratings>]
"""
    database = []
    open(path, "r") do fi
        database = readdlm(fi, '\t')
    end
    return database
end

function build_context(dict, tokens, window)
"""
    Builds a valid context array
"""
    context = tokens
    filter!(x -> x in dict.id2word, context)
    if length(context) > window
        context = context[1:window]
    end
    return context
end

function filter_context(dict, text, window, word)
"""
    Prepares context of size 2 * window with only valid words
"""
    split_text = split(text, "<b> " * word * " </b>") # remove current word
    left_text = reverse(split(split_text[1]))
    right_text = split(split_text[2])
    left_context = build_context(dict, left_text, window)
    right_context = build_context(dict, right_text, window)
    return vcat(left_context, right_context)
end

function evaluate_model(local_model, window, SCWS_PATH)
"""
    Evaluates the given model using a context of size window
"""

    data = load_SCWS(SCWS_PATH)
    vm, dict = load_model(local_model);

    Avg_list = Float32[]
    Max_list = Float32[]
    AvgC_list = Float32[]
    MaxC_list = Float32[]
    similarities_GS = Float32[]

    for i = 1:size(data)[1]
        entry = data[i,:]
        word1 = lowercase(string(entry[2]))
        word2 = lowercase(string(entry[4]))
        text1 = lowercase(entry[6])
        text2 = lowercase(entry[7])
        if !words_in_model(dict, word1, word2)
            #println("At least one word not in model for entry: ", i)
            continue
        end

        context1 = filter_context(dict, text1, window, word1)
        context2 = filter_context(dict, text2, window, word2)
        if length(context1) == 0 || length(context2) == 0
            println("Context invalid for entry id: ", i)
            continue
        end

        scores = evaluate_entry(vm, dict, word1, word2, context1, context2)
        push!(Avg_list, scores[1])
        push!(Max_list, scores[2])
        push!(AvgC_list, scores[3])
        push!(MaxC_list, scores[4])
        push!(similarities_GS, entry[8])
    end

    println("--------------------------------------------")
    println("Processed pairs: ", length(Avg_list))
    println("--------------------------------------------")
    println("Spearman's coeff for AvgSim: ", corspearman(Avg_list, similarities_GS))
    println("Spearman's coeff for MaxSim: ", corspearman(Max_list, similarities_GS))
    println("Spearman's coeff for AvgSimC: ", corspearman(AvgC_list, similarities_GS))
    println("Spearman's coeff for MaxSimC: ", corspearman(MaxC_list, similarities_GS))
    println("--------------------------------------------")
end

using ArgParse
using AdaGram
using StatsBase

s = ArgParseSettings()

@add_arg_table s begin
  "model"
    help = "File or directory with AdaGram model(s)"
    arg_type = AbstractString
    required = true
  "--window"
    help = "Size of context window"
    arg_type = Int64
    default = 4
end

args = parse_args(ARGS, s)

#SCWS_PATH = "/home/andres/SCWS_dataset/test"
SCWS_PATH = "/home/andres/SCWS_dataset/ratings_unquoted.txt"
win = args["window"]
model = args["model"]
temp_folder = "temp_folder"

# if model is single file, copy it in a temp folder
if !isdir(model)
    if isdir(temp_folder)
        rm(temp_folder, recursive = true)
    end
    mkdir(temp_folder)
    cp(model, temp_folder * "/" * model)
    model = temp_folder
end

for curr_model in readdir(model)
    println("Evaluating model: " * curr_model)
    evaluate_model(model * "/" * curr_model, win, SCWS_PATH)
end

# if temp folder was created, delete it
if isdir(temp_folder)
    rm(temp_folder, recursive=true)
end

max_stars_repo_name: akolonin/language-learning

max_issues_repo_name: akolonin/language-learning

max_forks_repo_name: akolonin/language-learning


text:
using PyCall: pyimport_conda
using Conda

function installpypackage()
	try
		pyimport_conda("sklearn", "scikit-surprise")
	catch
		try
			Conda.add("scikit-surprise", channel="conda-forge")
		catch
			println("scikit-learn failed to install")
		end
	end
end

installpypackage()