text:
reload("ODBC")

@show ODBC.listdrivers()
@show ODBC.listdsns()

using HTTP

db = HTTP.get("http://api.postgression.com")
dsn = ODBC.DSN("Driver={PostgreSQL};Server=")

# Check some basic queries

# setup a test database
# ODBC.execute!(dsn, "drop database if exists testdb")
# ODBC.execute!(dsn, "create database testdb")

# @test data.data[1][1] === Nullable(Int64(1))
# @test data.data[2][1] === Nullable(Int8(1))
# @test data.data[3][1] === Nullable(ODBC.DecFP.Dec64(1))
# @test data.data[4][1] === Nullable(Int32(1))
# @test data.data[5][1] === Nullable(ODBC.DecFP.Dec64(1))
# @test data.data[6][1] === Nullable(Int16(1))
# @test data.data[7][1] === Nullable(Int32(1))
# @test data.data[8][1] === Nullable(Int8(1))
# @test data.data[9][1] === Nullable(Float32(1.2))
# @test data.data[10][1] === Nullable(Float64(1.2))
# @test data.data[11][1] === Nullable(ODBC.API.SQLDate(2016,1,1))
# @test data.data[12][1] === Nullable(ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0))
# @test data.data[13][1] === Nullable(ODBC.API.SQLTimestamp(2016,1,1,1,1,1,0))
# @test data.data[14][1] === Nullable(ODBC.API.SQLTime(1,1,1))
# @test data.data[15][1] === Nullable(Int16(2016))
# @test string(get(data.data[16][1])) == "A"
# @test string(get(data.data[17][1])) == "hey there sailor"
# @test get(data.data[18][1]) == UInt8[0x31,0x32]
# @test isnull(data.data[19][1])
# @test get(data.data[20][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x61,0x62,0x72,0x61,0x68,0x61,0x6d]
# @test get(data.data[21][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x62,0x69,0x6c,0x6c]
# @test get(data.data[22][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x63,0x68,0x61,0x72,0x6c,0x69,0x65]
# @test get(data.data[23][1]) == UInt8[0x68,0x65,0x79,0x20,0x74,0x68,0x65,0x72,0x65,0x20,0x64,0x61,0x6e]
# @test string(get(data.data[24][1])) == "hey there ephraim"
# @test string(get(data.data[25][1])) == "hey there frank"
# @test string(get(data.data[26][1])) == "hey there george"
# @test string(get(data.data[27][1])) == "hey there hank"

ODBC.Source(dsn, "insert test1 VALUES
                    (1, -- bigint
                     1, -- bit
                     1.0, -- decimal
                     1, -- int
                     1.0, -- numeric
                     1, -- smallint
                     1, -- mediumint
                     1, -- tinyint
                     1.2, -- float
                     1.2, -- double
                     '2016-01-01', -- date
                     '2016-01-01 01:01:01', -- datetime
                     '2016-01-01 01:01:01-05:00', -- timestamp
                     '01:01:01', -- time
                     2016, -- year
                     'A', -- char(1)
                     'hey there sailor', -- varchar
                     cast(123456 as binary(2)), -- binary
                     NULL, -- varbinary
                     'hey there abraham', -- tinyblob
                     'hey there bill', -- blob
                     'hey there charlie', -- mediumblob
                     'hey there dan', -- longblob
                     'hey there ephraim', -- tinytext
                     'hey there frank', -- text
                     'hey there george', -- mediumtext
                     'hey there hank' -- longtext
                    )")
data = ODBC.query(dsn, "select * from test1")
@test size(data) == (2,27)

# ODBC.execute!(dsn, """
# CREATE TABLE test2
# (
#     ID INT NOT NULL PRIMARY KEY,
#     first_name VARCHAR(25),
#     last_name VARCHAR(25),
#     Salary DECIMAL,
#     `hourly rate` real,
#     hireDate DATE,
#     `last clockin` DATETIME
# );""")
# ODBC.execute!(dsn, "load data local infile '/Users/jacobquinn/Downloads/randoms.csv' into table test2
#                     fields terminated by ',' lines terminated by '\n'
#                     (id,first_name,last_name,salary,`hourly rate`,hiredate,`last clockin`)")

data = ODBC.query(dsn, "select count(*) from test2")
@test size(data) == (1,1)
@test data.data[1][1] === Nullable(70000)

@time data = ODBC.query(dsn, "select * from test2");
@test size(data) == (70000,7)
@test data.data[1].values == [1:70000...]
@test data.data[end][1] === Nullable(ODBC.API.SQLTimestamp(2002,1,17,21,32,0,0))

# test exporting test1 to CSV
source = ODBC.Source(dsn, "select * from test1")
csv = CSV.Sink("test1.csv")
Data.stream!(source, csv)
open("test1.csv") do f
    @test readline(f) == "\"test_bigint\",\"test_bit\",\"test_decimal\",\"test_int\",\"test_numeric\",\"test_smallint\",\"test_mediumint\",\"test_tiny_int\",\"test_float\",\"test_real\",\"test_date\",\"test_datetime\",\"test_timestamp\",\"test_time\",\"test_year\",\"test_char\",\"test_varchar\",\"test_binary\",\"test_varbinary\",\"test_tinyblob\",\"test_blob\",\"test_mediumblob\",\"test_longblob\",\"test_tinytext\",\"test_text\",\"test_mediumtext\",\"test_longtext\"\n"
    @test readline(f) == "1,1,+1E+0,1,+1E+0,1,1,1,1.2,1.2,2016-01-01,2016-01-01T01:01:01,2016-01-01T01:01:01,01:01:01,2016,\"A\",\"hey there sailor\",UInt8[49,50],\"\",UInt8[104,101,121,32,116,104,101,114,101,32,97,98,114,97,104,97,109],UInt8[104,101,121,32,116,104,101,114,101,32,98,105,108,108],UInt8[104,101,121,32,116,104,101,114,101,32,99,104,97,114,108,105,101],UInt8[104,101,121,32,116,104,101,114,101,32,100,97,110],\"hey there ephraim\",\"hey there frank\",\"hey there george\",\"hey there hank\"\n"
end
rm("test1.csv")

# # test exporting test2 to CSV
source = ODBC.Source(dsn, "select * from test2")
csv = CSV.Sink("test2.csv")
Data.stream!(source, csv)
open("test2.csv") do f
    @test readline(f) == "\"ID\",\"first_name\",\"last_name\",\"Salary\",\"hourly rate\",\"hireDate\",\"last clockin\"\n"
    @test readline(f) == "1,\"Lawrence\",\"Powell\",+87217E+0,26.47,2002-04-09,2002-01-17T21:32:00\n"
end

# test exporting test1 to SQLite
db = SQLite.DB()
source = ODBC.Source(dsn, "select * from test1")
sqlite = SQLite.Sink(source, db, "test1")
Data.stream!(source, sqlite)

data = SQLite.query(db, "select * from test1")
@test size(data) == (2,27)
@test data.data[1][1] === Nullable(1)
@test data.data[3][1] === Nullable(1.0)
@test data.data[11][1] === Nullable(ODBC.API.SQLDate(2016,1,1))

# test exporting test2 to SQLite
source = ODBC.Source(dsn, "select * from test2")
sqlite = SQLite.Sink(source, db, "test2")
Data.stream!(source, sqlite)

data = SQLite.query(db, "select * from test2")
@test size(data) == (70000,7)
@test data.data[1].values == [1:70000...]
@test data.data[end][1] === Nullable(ODBC.API.SQLTimestamp(2002,1,17,21,32,0,0))

ODBC.execute!(dsn, "drop table if exists test1")
# ODBC.Source(dsn, "drop table if exists test2")


text:
using Random, Distributions

numactions = 3
oppstrategy = [0.4, 0.3, 0.3];

function pos(x)
    x > 0.0 ? x : 0.0
end

function normalize(xs)
    total = sum(xs)
    n = length(xs)
    if total > 0
        return map(x -> x / total, xs)
    else
         return fill(1.0 / n, n)
    end
end

function regrets2strat(regrets)
    strategy = map(pos, regrets)
    normalize(strategy)
end

function action(strategy)
    d = Categorical(strategy)
    rand(d)
end

function utility(a, b)
    if a == b return 0 end
    if a == 1 && b == 3 return 1 end
    if b == 1 && a == 3 return -1 end
    a > b ? 1 : -1
end

function train(n, f1, f2)
    regrets = (zeros(numactions), zeros(numactions))
    strategys = (zeros(numactions), zeros(numactions))
    for _ in 1:n
        # Get regret-matched mixed-strategy actions
        strats = (f1(regrets[1]), f2(regrets[2]))
        strategys = strategys .+ strats
        actions = map(action, strats)
        # Accumulate action regrets
        for a in 1:numactions
            regrets[1][a] += utility(a, actions[2]) - utility(actions[1], actions[2])
            regrets[2][a] += -utility(actions[1], a) + utility(actions[2], actions[1])
        end
    end
    strats = map(regrets2strat, regrets)
    map(normalize, strategys .+ strats)
end

ext: jl

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


text:
module TestInterface
    using Test
    using CBSOData3

    include("mock_getodata.jl")

    tbls = CBSOData3.get_tables()
    @test length(tbls) > 0

end


text:
"""
    Battery
"""
struct Battery <: Resource
    index::Int           # Index of the resource
    num_timesteps::Int   # Number of time-steps in the battery's operation

    #==================================================
        Battery dynamics
    ==================================================#
    soc_min::Vector{Float64}     # Minimum state of charge for each time step
    soc_max::Vector{Float64}     # Maximum state of charge for each time step
    soc_init::Float64            # Initial state of charge
    selfdischargerate::Float64   # Self-discharging rate
    charge_pwr_min::Float64      # Minimum charging power rate
    charge_pwr_max::Float64      # Maximum charging power rate
    discharge_pwr_min::Float64   # Minimum discharging power rate
    discharge_pwr_max::Float64   # Maximum discharging power rate
    charge_eff::Float64          # Charging efficiency
    discharge_eff::Float64       # Discharging efficiency

    function Battery(;
        index::Integer=0,
        T::Integer=0,
        soc_min::Vector{Float64}=[0.0],
        soc_max::Vector{Float64}=[0.0],
        soc_init::Float64=0.0,
        selfdischargerate::Float64=1.0,
        charge_pwr_min::Float64=0.0,
        charge_pwr_max::Float64=0.0,
        discharge_pwr_min::Float64=0.0,
        discharge_pwr_max::Float64=0.0,
        charge_eff::Float64=0.0,
        discharge_eff::Float64=0.0
    )

        # Dimension checks
        T == size(soc_min, 1) || throw(DimensionMismatch("Invalid soc_min"))
        T == size(soc_max, 1) || throw(DimensionMismatch("Invalid soc_max"))
        0.0 <= charge_eff <= 1.0 || throw(DomainError("charge_eff must be in [0.0, 1.0]"))
        0.0 <= discharge_eff <= 1.0 || throw(DomainError("discharge_eff must be in [0.0, 1.0]"))
        0.0 <= selfdischargerate <= 1.0 || throw(DomainError("selfdischargerate must be in [0.0, 1.0]"))

        return new(
            index, T,
            copy(soc_min), copy(soc_max), soc_init,
            selfdischargerate,
            charge_pwr_min, charge_pwr_max,
            discharge_pwr_min, discharge_pwr_max,
            charge_eff, discharge_eff
        )
    end
end

function add_resource_to_model!(
    m::MOI.ModelLike,
    h::House,
    bat::Battery,
    var2idx, con2idx
)
    T = bat.num_timesteps

    # ==========================================
    #    I. Add local variables
    # ==========================================
    pnet = MOI.add_variables(m, T)  # battery net power
    pchg = MOI.add_variables(m, T)  # battery charging power
    pdis = MOI.add_variables(m, T)  # battery discharging power
    bsoc = MOI.add_variables(m, T)  # State of charge
    uchg = MOI.add_variables(m, T)  # Charge indicator
    udis = MOI.add_variables(m, T)  # Discharge indicator

    # Variable bounds
    for t in 1:T
        MOI.add_constraint(m, MOI.SingleVariable(pchg[t]), MOI.GreaterThan(0.0))
        MOI.add_constraint(m, MOI.SingleVariable(pdis[t]), MOI.GreaterThan(0.0))
        MOI.add_constraint(m, MOI.SingleVariable(bsoc[t]),
            MOI.Interval(bat.soc_min[t], bat.soc_max[t])
        )
        MOI.add_constraint(m, MOI.SingleVariable(uchg[t]), MOI.ZeroOne())
        MOI.add_constraint(m, MOI.SingleVariable(udis[t]), MOI.ZeroOne())
    end

    # Update variable indices
    for t in 1:T
        var2idx[(:bat, bat.index, :pnet, t)] = pnet[t]
        var2idx[(:bat, bat.index, :pchg, t)] = pchg[t]
        var2idx[(:bat, bat.index, :pdis, t)] = pdis[t]
        var2idx[(:bat, bat.index, :bsoc, t)] = bsoc[t]
        var2idx[(:bat, bat.index, :uchg, t)] = uchg[t]
        var2idx[(:bat, bat.index, :udis, t)] = udis[t]
    end

    # ==========================================
    #    II. Update linking constraints
    # ==========================================
    for t in 1:T
        cidx = con2idx[(:house, h.index, :link, t)]
        MOI.modify(m, cidx, MOI.ScalarCoefficientChange(pnet[t], -1.0))
    end

    # ==========================================
    #    III. Add local constraints
    # ==========================================
    for t in 1:T
        # Net power
        #   `p[t] - pchg[t] + pdis[t] = 0`
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm( 1.0, pnet[t]),
                MOI.ScalarAffineTerm(-1.0, pchg[t]),
                MOI.ScalarAffineTerm( 1.0, pdis[t])
            ], 0.0),
            MOI.EqualTo(0.0)
        )

        # Charging power
        #   `pchg[t] >= uchg[t] * P_chg_min`
        #   `pchg[t] <= uchg[t] * P_chg_max`
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm( 1.0, pchg[t]),
                MOI.ScalarAffineTerm(-bat.charge_pwr_min, uchg[t])
            ], 0.0),
            MOI.GreaterThan(0.0)
        )
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm( 1.0, pchg[t]),
                MOI.ScalarAffineTerm(-bat.charge_pwr_max, uchg[t])
            ], 0.0),
            MOI.LessThan(0.0)
        )

        # Discharging power
        #   `pdis[t] >= udis[t] * P_dis_min`
        #   `pdis[t] <= udis[t] * P_dis_max`
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm( 1.0, pdis[t]),
                MOI.ScalarAffineTerm(-bat.discharge_pwr_min, udis[t])
            ], 0.0),
            MOI.GreaterThan(0.0)
        )
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm( 1.0, pdis[t]),
                MOI.ScalarAffineTerm(-bat.discharge_pwr_max, udis[t])
            ], 0.0),
            MOI.LessThan(0.0)
        )

        # Charge-discharge
        #   `uchg[t] + udis[t] <= 1.0`
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(1.0, uchg[t]),
                MOI.ScalarAffineTerm(1.0, udis[t])
            ], 0.0),
            MOI.LessThan(1.0)
        )
    end

    # Energy conservation at time t=1
    #   soc[1] - soc_init = bat.charge_eff * pchg[1] - (1.0 / bat.discharge_eff) * pdis[1]
    MOI.add_constraint(m, MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, bsoc[1]),
            MOI.ScalarAffineTerm(-bat.charge_eff, pchg[1]),
            MOI.ScalarAffineTerm(1.0 / bat.discharge_eff, pdis[1]),
        ], 0.0),
        MOI.EqualTo(bat.soc_init * bat.selfdischargerate)
    )
    # Energy conservation at time t>1
    #   soc[t] - soc[t-1] == bat.charge_eff * pchg[t] - (1.0 / bat.discharge_eff) * pdis[t]
    for t in 2:T
        MOI.add_constraint(m, MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, bsoc[t]),
            MOI.ScalarAffineTerm(-bat.selfdischargerate, bsoc[t-1]),
            MOI.ScalarAffineTerm(-bat.charge_eff, pchg[t]),
            MOI.ScalarAffineTerm(1.0 / bat.discharge_eff, pdis[t]),
        ], 0.0),
        MOI.EqualTo(0.0)
    )
    end

    # ==========================================
    #    IV. Add sub-resources to current model
    # ==========================================
    # (None here)

    return nothing
end


text:
module Help

using Requests
using TextAnalysis
using Base.Markdown
using JLD

import Base: show, showall

const G_STEMMER = Stemmer("porter")
const HELP_DATA_PATH = joinpath(Pkg.dir("Help"), "helpdata.jld")

G_pkgnames = nothing
G_colidx = nothing
G_tfidf = nothing

if isfile(HELP_DATA_PATH)
    loaded = load(HELP_DATA_PATH)
    global G_pkgnames = loaded["G_pkgnames"]
    global G_colidx = loaded["G_colidx"]
    global G_tfidf = loaded["G_tfidf"]
end

stemlist(list) = map(word -> stem(G_STEMMER, word), list)

function process_strdoc(sd)
    todo = [strip_punctuation, strip_numbers, strip_case, strip_pronouns, strip_non_letters, strip_stopwords]
    for item in todo
n        prepare!(sd, item)
    end
    prepare!(sd, strip_patterns, skip_words = Set{AbstractString}(["julia", "package"]))
    return stemlist(tokens(sd))
end

function process_md_node!(node::Markdown.MD, vec)
    for item in node.content
        process_md_node!(item, vec)
    end
end

function process_md_node!(node::Markdown.Paragraph, vec)
    for item in node.content
        process_md_node!(item, vec)
    end
end

function process_md_node!(node::Markdown.List, vec)
    for bunch in node.items
        for item in bunch
            process_md_node!(item, vec)
        end
    end
end

process_md_node!(node::AbstractString, vec) = append!(vec, node |> StringDocument |> process_strdoc)

process_md_node!(node, vec) = nothing

function process_md(parsed)
    vec = AbstractString[]
    process_md_node!(parsed, vec)
    return vec
end

function process_pkgs()
    metadatadir = Pkg.dir("METADATA")
    files = split(readall(`find $metadatadir -name url`))
    pkgnames = AbstractString[]
    vecs = TokenDocument[]
    for file in files
        gitlink = open(file) do f; readall(f); end
        pkgname = match(r"/[^/]+/url", file).match[2:end-4]
        readme = ""
        try
            rawlink = "http://raw.githubusercontent.com" * match(Regex("/[^/]+/$pkgname"), gitlink).match * ".jl/master/README.md"
            readme = bytestring(get(rawlink).data)
        catch
            println("Errored: $pkgname")
            continue
        end
        push!(pkgnames, pkgname)
        push!(vecs, readme |> Markdown.parse |> process_md |> TokenDocument)
        println("Processed $pkgname $(length(readme))")
    end
    dtmat = getdtm(vecs)
    return pkgnames, dtmat.column_indices, tf_idf(dtmat)
end

function getdtm(vecs)
    cps = vecs |> Vector{GenericDocument} |> Corpus
    update_lexicon!(cps)
    return DocumentTermMatrix(cps)
end

function get_vector(idx, v)
    vec = zeros(Int, length(idx))
    for word in v
        try
            vec[idx[word]] += 1
        end
    end
    return vec
end

cosine(a, b::SparseMatrixCSC) = cosine(a, b |> full |> vec)
cosine(a, b) = dot(a, b) / (norm(a) * norm(b))

"""
    update()

Build the necessary data for searching packages.
"""
function update()
    global G_pkgnames, G_colidx, G_tfidf
    G_pkgnames, G_colidx, G_tfidf = process_pkgs()
    @save "$HELP_DATA_PATH" G_pkgnames G_colidx G_tfidf
    nothing
end

function get_exported_functions(mod)
    eval(parse("using $mod"))
    filter(x -> isa(eval(x), Function), names(mod))
end

process_funcsym(func) = func |> eval |> Base.doc |> process_md |> TokenDocument

function process_module(mod)
    funcnames = get_exported_functions(mod)
    tokens = map(process_funcsym, funcnames)
    dtmat = getdtm(tokens)
    return map(string, funcnames), dtmat.column_indices, tf_idf(dtmat)
end

type HelpResults
    results::Array{Tuple{AbstractString, Real}, 1}
end

Base.show(io::IO, h::HelpResults) = showall(h, num=5, io=io)

function Base.showall(h; num=typemax(Int), io=STDOUT)
    str = "Help Results\n---------\n"
    for i in 1:min(length(h.results), num)
        str = str * "$i. $(h.results[i][1])\n"
    end
    print(io, str)
end

process_user_query(query, colidx) = get_vector(colidx, query |> StringDocument |> process_strdoc)

function get_ratings(docnames, tfidf, qvec)
    ratings = [(docnames[i], cosine(qvec, tfidf[i, :])) for i in 1:size(tfidf)[1]]
    filter!(x -> x[2] != 0 && !isnan(x[2]), ratings)
    sort!(ratings, rev=true, lt = (x, y) -> x[2] < y[2])
    return ratings
end

function _help(query, funcnames, colidx, tfidf)
    qvec = process_user_query(query, colidx)
    ratings = get_ratings(funcnames, tfidf, qvec)
    if length(ratings) == 0
        println("Sorry, couldn't find any matching packages.")
        return nothing
    end
    return HelpResults(ratings[1 : threshold(ratings)])
end

"""
    help(query::AbstractString) -> HelpResults

Look for a package that best suits the query.
"""
function help(query)
    (G_pkgnames == nothing || G_colidx == nothing || G_tfidf == nothing) && error("Document vectors not built. Please run update()")
    return _help(query, G_pkgnames, G_colidx, G_tfidf)
end

# Return everything less than the average.
function threshold_avg(ratings)
    avg = sum(x -> x[2], ratings) / length(ratings)
    for i in length(ratings)
        ratings[i][2] < avg && return i - 1
    end
    return length(ratings)
end

function threshold(ratings)
    len = length(ratings)
    if len > 2
        return 2 + (map(x -> x[2], ratings) |> diff |> diff |> indmax)
    end
    return len
end

G_pkgdata = Dict()
PACKAGE_DATA_PATH = joinpath(Pkg.dir("Help"), "package_data")

get_pkg_data_file(pkg) = joinpath(PACKAGE_DATA_PATH, string(pkg) * ".jld")

"""
    update_pkg_data(pkg::Module)

Build search data for a package.
"""
function update_pkg_data(pkg)
    funcnames, colidx, tfidf = process_module(pkg)
    jldfile = get_pkg_data_file(pkg)
    @save "$jldfile" funcnames colidx tfidf
    G_pkgdata[pkg] = funcnames, colidx, tfidf
    nothing
end

function get_search_data(pkg)
    haskey(G_pkgdata, pkg) && return G_pkgdata[pkg]
    jldfile = get_pkg_data_file(pkg)
    if !isfile(jldfile)
        update_pkg_data(pkg)
    else
        loaded = load(jldfile)
        G_pkgdata[pkg] = loaded["funcnames"], loaded["colidx"], loaded["tfidf"]
    end
    return G_pkgdata[pkg]
end

"""
    help(pkg::Module, query::AbstractString) -> HelpResults

Look for a function in the package that best suits the query.

The package must be installed.
"""
function help(pkg, query)
    funcnames, colidx, tfidf = get_search_data(pkg)
    return _help(query, funcnames, colidx, tfidf)
end

export help, update, update_pkg_data

end # module


text:
using DataStructures

function main()
    H, W = parse.(Int, split(readline()))
    mp = [readline() for _ = 1:H]

    dx = [1, 0, -1, 0]
    dy = [0, 1, 0, -1]

    xₛ, yₛ, xₑ, yₑ = -1, -1, -1, -1
    for i = 1:H, j = 1:W
        if mp[i][j] == 's'
            xₛ, yₛ = i, j
        elseif mp[i][j] == 'g'
            xₑ, yₑ = i, j
        end
    end

    que = Deque{Tuple{Int, Int}}()
    pushfirst!(que, (xₛ, yₛ))
    dists = 10^9 * ones(Int, H, W)
    dists[xₛ, yₛ] = 0

    while !isempty(que)
        x, y = popfirst!(que)

        for d = 1:4
            x′ = x + dx[d]
            y′ = y + dy[d]

            if 1 ≤ x′ ≤ H && 1 ≤ y′ ≤ W
                if mp[x′][y′] == '#'
                    if dists[x′, y′] > dists[x, y] + 1
                        dists[x′, y′] = dists[x, y] + 1
                        push!(que, (x′, y′))
                    end
                else
                    if dists[x′, y′] > dists[x, y]
                        dists[x′, y′] = dists[x, y]
                        pushfirst!(que, (x′, y′))
                    end
                end
            end
        end
    end

    if dists[xₑ, yₑ] ≤ 2
        println("YES")
    else
        println("NO")
    end

end

main()


text:
# -*- coding: utf-8 -*-
# + {}
""" SemanticModels

provides the AbstractModel type and constructors for building hierarchical model representations.

"""
module SemanticModels

using Unitful
using Unitful: Units

export Parsers, Extraction, Graphs

abstract type AbstractModel end

abstract type Model end
abstract type EpiModel <: Model end
#abstract type Unit end
abstract type NumberClass end
abstract type Amount <: NumberClass end
abstract type Rate <: NumberClass end
abstract type BirthRate <: Rate end
abstract type DeathRate <: Rate end
abstract type TransitionRate <: Rate end
abstract type Equation end
abstract type Expression end
abstract type Variable end

struct NumParameter{U, C}
    name::Symbol
    unit::U
    class::C
end

struct NumVariable{U, C}
    name::Symbol
    unit::U
    class::C
end

struct ODE{P,V} <: Equation
    f::Expr
    parameters::Vector{P}
    variables::Vector{V}
end

struct SIR{E, Vp, Vt} <: EpiModel
    equations::Vector{E}
    population::Vp
    domain::Vt
end

@unit person "person" Person 1 true
Unitful.register(SemanticModels)
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(SemanticModels)
end

"""    BasicSIR

defines a representation of Susceptible Infected Recovered models using Unitful
numbers and expressions.

"""
function BasicSIR()
    β = NumParameter(:β, u"person/s", TransitionRate)
    γ = NumParameter(:γ, u"person/s", TransitionRate)
    @show β, γ

    S = NumVariable(:S, u"person", Amount)
    I = NumVariable(:I, u"person", Amount)
    R = NumVariable(:R, u"person", Amount)
    sir = SIR([ODE(:(dS/dt -> -βSI/N),
                [β], [S]),
                ODE(:(dI/dt -> βSI/N -γI),
                    [β, γ], [S, I]),
            ODE(:(dR/dt -> γI), [γ], [I])],
            NumVariable.([:N,:S,:I,:R], u"person", Amount),
            NumVariable(:t, u"s", Amount))
    return sir
end

include("cassette.jl")

"""    CombinedModel

represents a model that has a fixed set of dependencies.
It is the basic building block of a model DAG.

"""
mutable struct CombinedModel{F,S} <: Model
    deps::F
    target::S
end

"""    solve(m::AbstractModel)

executes the solving of a model for models with dependencies, those deps are
executed first and then the node is solved. This is the function that evaluates
the model DAG.

"""
function solve(m::CombinedModel)
    return solve(m.target(m, solve.(m.deps)))
end

include("parse.jl")
include("modeltools/ModelTools.jl")

end
# -


text:
"""
The `AbstractOptimizer` interface allows to specify optimization procedures.
It consists of three methods:
* `initialized(optimizer, bm)`: May be used for creating an optimizer that is
   specifically initialized for the Boltzmann machine `bm`.
   In particular it may be used to allocate reusable space for the gradient.
   The default implementation simply returns the unmodified `optimizer`.
* `computegradient!(optimizer, v, vmodel, h, hmodel, rbm)` or
  `computegradient!(optimizer, meanfieldparticles, gibbsparticles, dbm)`
   needs to be implemented for computing the gradient given the samples
   from the positive and negative phase.
* `updateparameters!(bm, optimizer)` needs to be specified for taking the
   gradient step. The default implementation for RBMs expects the fields
   `learningrate` and `gradient` and adds `learningrate * gradient` to the
   given RBM.
"""
abstract type AbstractOptimizer{R<:AbstractBM}
end

abstract type
      AbstractLoglikelihoodOptimizer{R<:AbstractRBM} <: AbstractOptimizer{R}
end

struct NoOptimizer <: AbstractOptimizer{AbstractBM}
end

"""
Implements the `AbstractOptimizer` interface for optimizing the loglikelihood
with stochastic gradient descent.
"""
mutable struct LoglikelihoodOptimizer{R<:AbstractRBM} <: AbstractLoglikelihoodOptimizer{R}
   gradient::R
   negupdate::Matrix{Float64}
   learningrate::Float64
   sdlearningrate::Float64
end

function converttodbmoptimizer(optimizer::NoOptimizer, dbm::MultimodalDBM)
   optimizer
end

function converttodbmoptimizer(optimizer::AbstractOptimizer{<:AbstractRBM},
      dbm::MultimodalDBM)
   StackedOptimizer(map(i -> deepcopy(optimizer), eachindex(dbm)))
end

function converttodbmoptimizer(optimizer::AbstractOptimizer{MultimodalDBM},
      dbm::MultimodalDBM)
   optimizer
end

function loglikelihoodoptimizer(;
      learningrate::Float64 = 0.0, sdlearningrate::Float64 = 0.0)

   LoglikelihoodOptimizer(NoRBM(), Matrix{Float64}(undef, 0,0),
         learningrate, sdlearningrate)
end

function loglikelihoodoptimizer(rbm::R;
      learningrate::Float64 = 0.0, sdlearningrate::Float64 = 0.0) where {R<:AbstractRBM}

   LoglikelihoodOptimizer{R}(deepcopy(rbm),
         Matrix{Float64}(undef, size(rbm.weights)),
         learningrate, sdlearningrate)
end

function loglikelihoodoptimizer(prbm::PartitionedRBM;
      learningrate::Float64 = 0.0, sdlearningrate::Float64 = 0.0)

   PartitionedOptimizer(map(
         rbm -> loglikelihoodoptimizer(rbm;
               learningrate = learningrate, sdlearningrate = sdlearningrate),
         prbm.rbms))
end

struct PartitionedOptimizer <: AbstractOptimizer{PartitionedRBM}
   optimizers::Vector{AbstractOptimizer}
end

mutable struct BeamAdversarialOptimizer{R<: AbstractRBM} <: AbstractOptimizer{R}
   gradient::R
   negupdate::Matrix{Float64}
   critic::Vector{Float64}
   learningrate::Float64
   sdlearningrate::Float64
   knearest::Int
end

"""
    StackedOptimizer(optimizers)
Can be used for optimizing a stack of RBMs / a DBM by using the given the vector
of `optimizers` (one for each RBM).
For more information about the concept of optimizers, see `AbstractOptimizer`.
"""
struct StackedOptimizer <: AbstractOptimizer{MultimodalDBM}
   optimizers::Vector{AbstractOptimizer{<:AbstractRBM}}
end

struct CombinedOptimizer{R<: AbstractRBM,
         G1 <: AbstractOptimizer{R},
         G2 <: AbstractOptimizer{R}} <: AbstractOptimizer{R}

   part1::G1
   part2::G2
   gradient::R
   weight1::Float64
   weight2::Float64
   learningrate::Float64
   sdlearningrate::Float64
end

# TODO document
function beamoptimizer(;learningrate::Float64 = 0.05,
      sdlearningrate::Float64 = 0.0,
      adversarialweight::Float64 = 0.1,
      knearest::Int = 5)

   llstep = loglikelihoodoptimizer(
         learningrate = learningrate, sdlearningrate = sdlearningrate)
   advstep = BeamAdversarialOptimizer(
         NoRBM(), Matrix{Float64}(0,0), Vector{Float64}(),
         learningrate, sdlearningrate, knearest)

   CombinedOptimizer(advstep, llstep, NoRBM(),
         adversarialweight, 1.0 - adversarialweight,
         learningrate, sdlearningrate)
end

"""
    initialized(optimizer, rbm)
Returns an `AbstractOptimizer` similar to the given `optimizer`
that can be used to optimize the `AbstractRBM` `rbm`.
"""
function initialized(optimizer::AbstractOptimizer, rbm::AbstractRBM)
   optimizer
end

function initialized(optimizer::LoglikelihoodOptimizer, rbm::R
      ) where {R <: AbstractRBM}

   LoglikelihoodOptimizer(deepcopy(rbm),
         Matrix{Float64}(undef, size(rbm.weights)),
         optimizer.learningrate, optimizer.sdlearningrate)
end

function initialized(optimizer::CombinedOptimizer, rbm::R
      ) where {R <: AbstractRBM}

   CombinedOptimizer(initialized(optimizer.part1, rbm),
         initialized(optimizer.part2, rbm),
         deepcopy(rbm),
         optimizer.weight1, optimizer.weight2,
         optimizer.learningrate, optimizer.sdlearningrate)
end

function initialized(optimizer::BeamAdversarialOptimizer{R1}, rbm::R2
      ) where {R1 <: AbstractRBM, R2 <: AbstractRBM}

   BeamAdversarialOptimizer{R2}(deepcopy(rbm),
         Matrix{Float64}(undef, size(rbm.weights)),
         Vector{Float64}(),
         optimizer.learningrate, optimizer.sdlearningrate,
         optimizer.knearest)
end

function initialized(optimizer::PartitionedOptimizer, prbm::PartitionedRBM)
   PartitionedOptimizer(map(
         i -> initialized(optimizer.optimizers[i], prbm.rbms[i]),
         eachindex(prbm)))
end

function initialized(optimizer::AbstractOptimizer{R}, dbm::MultimodalDBM
      ) where {R <: AbstractRBM}
   # transform optimizer for RBM into stacked optimizer for DBM
   StackedOptimizer(map(rbm -> initialized(optimizer, rbm), dbm))
end

function initialized(stackedoptimizer::StackedOptimizer, dbm::MultimodalDBM)
   StackedOptimizer(map(
         i -> initialized(stackedoptimizer.optimizers[i], dbm[i]),
         eachindex(dbm)))
end

"""
    computegradient!(optimizer, v, vmodel, h, hmodel, rbm)
Computes the gradient of the RBM `rbm` given the
the hidden activation `h` induced by the sample `v`
and the vectors `vmodel` and `hmodel` generated by sampling from the model.
The result is stored in the `optimizer` in such a way that it can be applied
by a call to `updateparameters!`. There is no return value.

For RBMs (excluding PartitionedRBMs), this means saving the gradient
in a RBM of the same type in the field `optimizer.gradient`.
"""
function computegradient!(
      optimizer::AbstractLoglikelihoodOptimizer{R},
      v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::R
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2},
            R <: AbstractRBM}

   computegradientsweightsandbiases!(optimizer, v, vmodel, h, hmodel, rbm)
end

function computegradient!(
      optimizer::LoglikelihoodOptimizer{GaussianBernoulliRBM},
      v::M1, vmodel::M1, h::M2, hmodel::M2, gbrbm::GaussianBernoulliRBM
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2}}

   # See bottom of page 15 in [Krizhevsky, 2009].

   if optimizer.sdlearningrate > 0.0
      optimizer.gradient.sd .=
            sdupdateterm(gbrbm, v, h) - sdupdateterm(gbrbm, vmodel, hmodel)
   end

   v = v ./ gbrbm.sd'
   vmodel = vmodel ./ gbrbm.sd'

   computegradientsweightsandbiases!(optimizer, v, vmodel, h, hmodel, gbrbm)
   optimizer.gradient.visbias ./= gbrbm.sd

   nothing
end

function computegradient!(
      optimizer::AbstractLoglikelihoodOptimizer{GaussianBernoulliRBM2},
      v::M1, vmodel::M1, h::M2, hmodel::M2, gbrbm::GaussianBernoulliRBM2
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2}}

   # See Cho,
   # "Improved learning of Gaussian-Bernoulli restricted Boltzmann machines"
   sdsq = gbrbm.sd .^ 2

   if optimizer.sdlearningrate > 0.0
      sdgrads = vmodel .* (hmodel * gbrbm.weights')
      sdgrads .-= v .* (h * gbrbm.weights')
      sdgrads .*= 2.0
      sdgrads .+= (v .- gbrbm.visbias') .^ 2
      sdgrads .-= (vmodel .- gbrbm.visbias') .^ 2
      optimizer.gradient.sd .= vec(mean(sdgrads, dims = 1))
      optimizer.gradient.sd ./= sdsq
      optimizer.gradient.sd ./= gbrbm.sd
   end

   v = v ./ sdsq'
   vmodel = vmodel ./ sdsq'

   computegradientsweightsandbiases!(optimizer, v, vmodel, h, hmodel, gbrbm)
end

function computegradientsweightsandbiases!(
      optimizer::AbstractLoglikelihoodOptimizer,
      v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::R
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2},
            R <: AbstractRBM}

   npossamples = size(v, 1)
   nnegsamples = size(vmodel, 1)

   mul!(optimizer.gradient.weights, transpose(v), h)
   if npossamples != 1
      optimizer.gradient.weights ./= npossamples
   end
   mul!(optimizer.negupdate, transpose(vmodel), hmodel)
   if nnegsamples != 1
      optimizer.negupdate ./= nnegsamples
   end
   optimizer.gradient.weights .-= optimizer.negupdate

   optimizer.gradient.hidbias .= vec(mean(h, dims = 1) - mean(hmodel, dims = 1))
   optimizer.gradient.visbias .= vec(mean(v, dims = 1) - mean(vmodel, dims = 1))

   nothing
end

function computegradientsweightsandbiases!(
   optimizer::AbstractLoglikelihoodOptimizer,
   v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::R
   ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2},
         R <: NegativeBinomialBernoulliRBM}

   npossamples = size(v, 1)
   nnegsamples = size(vmodel, 1)

   mul!(optimizer.gradient.weights, transpose(v), h)
   if npossamples != 1
      optimizer.gradient.weights ./= npossamples
   end
   mul!(optimizer.negupdate, transpose(vmodel), hmodel)
   if nnegsamples != 1
      optimizer.negupdate ./= nnegsamples
   end
   optimizer.gradient.weights .-= optimizer.negupdate

   optimizer.gradient.hidbias .= vec(mean(h, dims = 1) - mean(hmodel, dims = 1))
   optimizer.gradient.visbias .= vec(mean(v, dims = 1) - mean(vmodel, dims = 1))
   nothing
end

function computegradient!(
      optimizer::BeamAdversarialOptimizer{R},
      v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::AbstractRBM
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2},
            R <:AbstractRBM}

   optimizer.critic = nearestneighbourcritic(h, hmodel, optimizer.knearest)

   nvisible = nvisiblenodes(rbm)
   nhidden = nhiddennodes(rbm)

   for i = 1:nvisible
      for j = 1:nhidden
         optimizer.gradient.weights[i, j] =
               cov(optimizer.critic, vmodel[:, i] .* hmodel[:, j])
      end
   end

   for i = 1:nvisible
      optimizer.gradient.visbias[i] = cov(optimizer.critic, vmodel[:, i])
   end

   for j = 1:nhidden
      optimizer.gradient.hidbias[j] = cov(optimizer.critic, hmodel[:, j])
   end

   nothing
end

function computegradient!(
      optimizer::BeamAdversarialOptimizer{GaussianBernoulliRBM2},
      v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::GaussianBernoulliRBM2
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2}}

   invoke(computegradient!,
         Tuple{BeamAdversarialOptimizer{GaussianBernoulliRBM2}, M1, M1, M2, M2, AbstractRBM},
         optimizer, v, vmodel, h, hmodel, rbm)

   sdsq = rbm.sd .^ 2

   weights = optimizer.gradient.weights
   visbias = optimizer.gradient.visbias
   nvisible = size(weights, 1)

   weights ./= sdsq
   visbias ./= sdsq

   if optimizer.sdlearningrate > 0
      for i = 1:nvisible
         optimizer.gradient.sd[i] =
               cov(optimizer.critic, (vmodel[:, i] .- visbias[i]) .^ 2 -
                     2.0 .* vmodel[:, i] .* (hmodel * weights[i, :]))
      end
   end

   optimizer.gradient.sd ./= sdsq .* rbm.sd

   nothing
end

function computegradient!(optimizer::CombinedOptimizer{R},
      v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::AbstractRBM
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2},
            R <:AbstractRBM}

   computegradient!(optimizer.part1, v, vmodel, h, hmodel, rbm)
   computegradient!(optimizer.part2, v, vmodel, h, hmodel, rbm)

   grad1 = optimizer.part1.gradient
   grad2 = optimizer.part2.gradient

   optimizer.gradient.weights .=
         grad1.weights .* optimizer.weight1 .+
         grad2.weights .* optimizer.weight2

   optimizer.gradient.visbias .=
         grad1.visbias .* optimizer.weight1 .+
         grad2.visbias .* optimizer.weight2

   optimizer.gradient.hidbias .=
         grad1.hidbias .* optimizer.weight1 .+
         grad2.hidbias .* optimizer.weight2

   nothing
end

function computegradient!(optimizer::CombinedOptimizer{R},
      v::M1, vmodel::M1, h::M2, hmodel::M2, rbm::R
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2},
            R <:Union{GaussianBernoulliRBM, GaussianBernoulliRBM2}}

   invoke(computegradient!,
         Tuple{CombinedOptimizer{R}, M1, M1, M2, M2, AbstractRBM},
         optimizer, v, vmodel, h, hmodel, rbm)

   if optimizer.sdlearningrate > 0.0
      optimizer.gradient.sd .=
            optimizer.part1.gradient.sd .* optimizer.weight1 .+
            optimizer.part2.gradient.sd .* optimizer.weight2
   end

   nothing
end

function computegradient!(optimizer::PartitionedOptimizer,
      v::M1, vmodel::M1, h::M2, hmodel::M2, prbm::PartitionedRBM
      ) where {M1 <: AbstractArray{Float64, 2}, M2 <: AbstractArray{Float64, 2}}

   for i in eachindex(prbm.rbms)
      computegradient!(optimizer.optimizers[i],
            view(v, :, prbm.visranges[i]),
            view(vmodel, :, prbm.visranges[i]),
            view(h, :, prbm.hidranges[i]),
            view(hmodel, :, prbm.hidranges[i]),
            prbm.rbms[i])
   end
   nothing
end

function computegradient!(stackedoptimizer::StackedOptimizer,
      meanfieldparticles::Particles, gibbsparticles::Particles,
      dbm::MultimodalDBM)

   for i in eachindex(dbm)
      computegradient!(stackedoptimizer.optimizers[i],
            meanfieldparticles[i], gibbsparticles[i],
            meanfieldparticles[i + 1], gibbsparticles[i + 1],
            dbm[i])
   end
   nothing
end

"""
    updateparameters!(rbm, optimizer)
Updates the RBM `rbm` by walking a step in the direction of the gradient that
has been computed by calling `computegradient!` on `optimizer`.
"""
function updateparameters!(rbm::AbstractRBM,
      optimizer::AbstractOptimizer{R}) where {R <: AbstractRBM}

   updateweightsandbiases!(rbm, optimizer)
end

function updateparameters!(rbm::Binomial2BernoulliRBM,
      optimizer::AbstractOptimizer{Binomial2BernoulliRBM})

   # To train a Binomial2BernoulliRBM exactly like
   # training a BernoulliRBM where each two nodes share the weights,
   # use half the learning rate in the visible nodes.
   learningratehidden = optimizer.learningrate
   learningrate = optimizer.learningrate / 2.0

   rbm.weights .+= learningrate .* optimizer.gradient.weights
   rbm.visbias .+= learningrate .* optimizer.gradient.visbias
   rbm.hidbias .+= learningratehidden .* optimizer.gradient.hidbias
   rbm
end

"""
Update Parameters of Negative Binomial RBM
Restrict Parameters according to
"Exponential Family Restricted Boltzmann Machines and Annealed Importance
Sampling" by Li et al.
"""

function updateparameters!(rbm::NegativeBinomialBernoulliRBM,
      optimizer::AbstractOptimizer{NegativeBinomialBernoulliRBM})

   updateweightsandbiases!(rbm, optimizer)
   rbm.weights .= min.(rbm.weights, 0.0)
   rbm.visbias .= min.(rbm.visbias, -10.0^-10.0)
   rbm
end

function updateparameters!(rbm::R, optimizer::AbstractOptimizer{R}
      ) where {R <: Union{GaussianBernoulliRBM, GaussianBernoulliRBM2}}

   updateweightsandbiases!(rbm, optimizer)

   if optimizer.sdlearningrate > 0.0
      rbm.sd .+= optimizer.sdlearningrate * optimizer.gradient.sd
   end
   rbm
end

function updateparameters!(prbm::PartitionedRBM, optimizer::PartitionedOptimizer)
   for i in eachindex(prbm.rbms)
      updateparameters!(prbm.rbms[i], optimizer.optimizers[i])
   end
   prbm
end

function updateparameters!(dbm::MultimodalDBM, stackedoptimizer::StackedOptimizer)
   for i in eachindex(dbm)
      updateparameters!(dbm[i], stackedoptimizer.optimizers[i])
   end
   dbm
end

function updateweightsandbiases!(rbm::R,
      optimizer::AbstractOptimizer{R}) where {R <: AbstractRBM}

   rbm.weights .+= optimizer.learningrate .* optimizer.gradient.weights
   rbm.visbias .+= optimizer.learningrate .* optimizer.gradient.visbias
   rbm.hidbias .+= optimizer.learningrate .* optimizer.gradient.hidbias
   rbm
end


text:
# This is a general type for Adaptive selection based on Area-under-curve
# multi-armed bandit as described in the paper:
#   Francois-Michel De Rainville et al, "Sustainable Cooperative Coevolution
#   with a Multi-Armed Bandit", 2013
#

# An AdaptiveSelector selects between a set of integers that represent different
# operators, values etc. It does not have any info about what the things it selects
# between are, below we call them items. The set of items the selector selects
# from is typically static but might not be.
abstract type AdaptiveSelector end

# The default AdaptiveSelector is a random selector.
mutable struct RandomSelector <: AdaptiveSelector
  item_map::Dict{Int, Any}
  next_index::Int

  RandomSelector(items = []) = begin
    rs = new(Dict{Int, Any}(), 0)
    add_items(rs, items)
    rs
  end
end

# Returns an item map that maps from the indices of the current item set to
# each item.
function item_map(a::AdaptiveSelector)
  a.item_map # Override this if you do not save the item map in an instance var
end

# Get next free index. Returns unique index number each time it is called.
function get_next_free_index(a::AdaptiveSelector)
  a.next_index += 1
end

# Add an item. Returns its index.
function add_item(a::AdaptiveSelector, item)
  index = get_next_free_index(a)
  item_map(a)[index] = item
  index
end

function add_items(a::AdaptiveSelector, items)
  [add_item(a, item) for item in items]
end

# Remove an item given its index.
function remove_item_with_index(a::AdaptiveSelector, item_index::Int)
  delete!(a.item_map, item_index)
end

# Remove an item.
function remove_item(a::AdaptiveSelector, item)
  not_equals_item = (k,v) -> v != item
  filter!( not_equals_item, a.item_map )
end

# Select one item in the set of things that we select from. Returns both its
# index and itself / its value.
function select(a::AdaptiveSelector)
  ks = keys(a.item_map)
  randkey = ks[rand(1:length(ks))]
  a.item_map[randkey]
end

# Add reward for one item given its index. Default is to not do anything since
# the default selector is a random selector that need not care about rewards.
# A reward should be a positive value with larger values indicating more
# benefit when using the item.
function reward(a::AdaptiveSelector, index, reward)
  # Override for selectors that use rewards.
end

# We need a queue of fixed size to keep info about recent rewards.
mutable struct FixedSizeQueue
  queue::Array{Any,1}
  size::Int
  maxsize::Int
  index::Int

  FixedSizeQueue(maxsize) = begin
    q = Any[]
    [Base.push!(q, 0) for i in 1:maxsize]
    new(q, 0, maxsize, 0)
  end
end

size(fsq::FixedSizeQueue) = fsq.size
maxsize(fsq::FixedSizeQueue) = fsq.maxsize
is_full(fsq::FixedSizeQueue) = fsq.size == fsq.maxsize

function push!(fsq::FixedSizeQueue, item)
  fsq.queue[fsq.index + 1] = item
  fsq.index = mod(fsq.index + 1, fsq.maxsize)
  fsq.size = min(fsq.maxsize, fsq.size + 1)
end

function as_vector(fsq::FixedSizeQueue)
  v = Any[]
  for elem = fsq
    push!(v, elem)
  end
  v
end

##############
# Define FixedSizeQueue to be an iterator by defining start, done and next:
##############

Base.start(fsq::FixedSizeQueue) = (fsq.size, mod(fsq.index - fsq.size, fsq.maxsize))

Base.done(fsq::FixedSizeQueue, state) = state[1] == 0

function Base.next(fsq::FixedSizeQueue, state)
  elem = fsq.queue[state[2] + 1]
  return elem, (state[1]-1, mod(state[2]+1, fsq.maxsize))
end

mutable struct AucBanditSelector <: AdaptiveSelector
  window::FixedSizeQueue  # Window of latest (itemindex, reward) pairs
  ns::Int[]             # Number of times each item has been used within the current window.
  items::Any[]
  next_index::Int
  decay_factor::Float64
  exploration_factor::Float64

  AucBanditSelector(items = []; window_size = 50, decay_factor = 1.0,
    exploration_factor = 1.0) = begin

    aucb = new(FixedSizeQueue(window_size), Int[], Any[], 0,
      decay_factor, exploration_factor)

    add_items(aucb, items)

    aucb

  end
end

function add_item(a::AucBanditSelector, item)
  push!(a.ns, 0) # Add a zero to indicate how often it has been used => none!
  push!(a.items, item)
  length(a.items)
end

# Remove an item given its index.
function remove_item_with_index(a::AucBanditSelector, item_index::Int)
  splice!(a.items, item_index)
  splice!(a.ns, item_index)
end

# Remove an item.
function remove_item(a::AucBanditSelector, item)
  for i in 1:length(a.items)
    if a.items[i] == item
      remove_item_with_index(a, i)
      return i
    end
  end
end

function reward(a::AucBanditSelector, index, reward)
  # Update the ns counts by subtracting one from the one that is about the exit
  # the window.
  if is_full(a.window)
    (index, reward) = first(a.window)
    a.ns[index] -= 1
  end
  push!(a.window, (index, reward))
  a.ns[index] += 1
end

function select(a::AucBanditSelector)
  unused = a.ns == 0
  if any(unusued)
    # Randomly select one of the unused
    indices_to_unused = collect(1:length(a.ns))[unused]
    indices_to_unused[rand(1:length(indices_to_unused))]
  else
    # Select via the multi-armed bandit formula
    qs = calc_quality_values(a)
    vs = qs + a.exploration_factor * sqrt( 2 * log(sum(a.ns)) * (1 / a.ns) )
    items[argmin(vs)]
  end
end

function calc_quality_values(a:AucBanditSelector)
  qs = zeros(length(a.ns))
  for i in 1:length(a.ns)
    qs[i] = area_under_curve_credit_assignment(a, i)
  end
  qs
end

function sort_window(a::AucBanditSelector)
  sort(as_vector(a.window), by = (index, reward) -> reward, rev = true)
end

# We should extend this to calculate all credits in one go instead of sorting
# and looping for each index...
function area_under_curve_credit_assignment(a::AucBanditSelector, index)
  sorted_window = sort_window(a)
  wlength = length(sorted_window)
  y = q = 0.0
  rank = 1
  for (j, reward) in sorted_window
    # The AUC bandit skips the actual reward values and relies only on rank...
    reward = a.decay_factor^(rank-1) * (wlength - (rank - 1))
    if index == j
      y += reward
    else
      q += (y*reward)
    end
    rank += 1
  end
  q
end


text:
module DummyTapir
using Base.Threads: @spawn
const Output = Ref{Any}
end

module TapirFolds
using Base: Tapir
include("folds_impl.jl")
end

module ThreadsFolds
using ..TapirBenchmarks: DummyTapir
const Tapir = DummyTapir
include("folds_impl.jl")
end


text:
using Flux

# ResidualBlock is used with `SkipConnection`
# using Parallel, first extract features, then combine them together
# separate filters
struct ResidualBlock
    block
end
Flux.@functor ResidualBlock
(b::ResidualBlock)(x) = x |> b.block .|> leakyrelu

function basic_block()
    layer = Chain(Conv((3, 3, 3), 16 => 16, leakyrelu; stride = 1, pad = 1),
                  Conv((3, 3, 3), 16 => 16; stride = 1, pad = 1))
    return ResidualBlock(SkipConnection(layer, +))
end

function residual_block()
    layer = Chain(Conv((1, 1, 1), 8 => 16, leakyrelu; stride = 1),
                  basic_block(),
                  Conv((1, 1, 1), 16 => 8; stride = 1))
    return ResidualBlock(SkipConnection(layer, +))
end

function residual_net()

    init_conv = Conv((1, 1, 1), 1 => 8, leakyrelu; stride = 1)
    resblock1 = residual_block()
    resblock2 = residual_block()
    resblock3 = residual_block()
    final_conv = Conv((1, 1, 1), 8 => 1; stride = 1)

    return Chain(init_conv, resblock1, resblock2, resblock3, final_conv)
end

# final_conv = Conv((1, 1, 1), 16 => 2; stride = 1)
# basic_conv = Chain(Conv((3, 3, 3), 16 => 16, leakyrelu; stride = 1, pad = 1),
#               Conv((3, 3, 3), 16 => 16; stride = 1, pad = 1))
#
# function mymodel()
#         return Chain(Parallel((x, y) -> cat(x, y, dims = 4),
#                      residual_net(),
#                      residual_net(),
#                      ),
#                      basic_block(),
#                      basic_block(),
#                      basic_block(),
#                      final_conv
#                      )
# end

# xx = randn(Float32, 16, 16, 16, 2, 1)
# xtrue = ones(Float32, 16, 16, 16, 2, 1)
# model = residual_net()
# model(xx)
# ps = Flux.params(model)
# opt = ADAM(0.001)
# loss(x, y) = Flux.Losses.mse(model(x), y)
# evalcb() = @show(loss(xx, xtrue))
# data = [(xx, xtrue)]
# @epochs 500 Flux.train!(loss, ps, data, opt, cb = evalcb)


text:
#####
##### Simple CLI for AlphaZero.jl
##### This file can also be included directly in the REPL.
#####

# ENV["CUARRAYS_MEMORY_LIMIT"] = 7_500_000_000
ENV["JULIA_CUDA_MEMORY_POOL"] = "split" # "binned" / "split"

# Enables running the script on a distant machine without an X server
ENV["GKSwstype"]="nul"

using AlphaZero
include("games.jl")

#####
##### Parse arguments
#####

using ArgParse
argstab = ArgParseSettings()
@add_arg_table! argstab begin
  "--game"
    help = "Select a game ($(join(AVAILABLE_GAMES, "/")))"
  "train"
    action = :command
    help = "Resume the training session"
  "play"
    action = :command
    help = "Play against the current system"
  "explore"
    action = :command
    help = "Use the interactive exploration system"
  "replot"
    action = :command
    help = "Regenerate all the session plots from the JSON reports"
  "--save-intermediate"
    action = :store_true
    help = "Save all intermediate states during training"
end
args = parse_args(isempty(ARGS) ? ["train"] : ARGS, argstab)
!isnothing(args["game"]) && (ENV["GAME"] = args["game"])
cmd = args["%COMMAND%"]

if !haskey(ENV, "GAME")
  println(stderr, "You must specify a game.")
  exit()
end

#####
##### Main
#####

GAME = ENV["GAME"]
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

const SESSION_DIR = joinpath("sessions", GAME)

params = Training.params
netparams = Training.netparams
benchmark = Training.benchmark

include("lib/dummy_run.jl")

if get(ENV, "DUMMY_RUN", "false") == "true"
  @warn "Running dummy run"
  params, benchmark = dummy_run_params(Training.params, Training.benchmark)
end

session = Session(
  Game, Training.Network{Game}, params, netparams, benchmark=benchmark,
  dir=SESSION_DIR, autosave=true, save_intermediate=args["save-intermediate"])

if cmd == "train"
  resume!(session)
elseif cmd == "explore"
  start_explorer(session)
elseif cmd == "play"
  play_interactive_game(session)
elseif cmd == "replot"
  AlphaZero.UserInterface.regenerate_plots(session)
end


text:
cd(@__DIR__)

using Pkg
Pkg.activate(".")
Pkg.instantiate()
##

using DataFrames
using CSV
using Plots
using StatsPlots
using Plots.PlotMeasures
using MCPhylo
using LinearAlgebra
using RCall
using Formatting
using ProgressMeter
using Statistics
using Random
Random.seed!(818480685100856001)

##

fpairs = CSV.read("../data/fpairs.txt", DataFrame, header=false)[:,1]

##
function getEqulibirium(rates)
      q = zeros(4, 4)
      indices = [
            (2, 1),
            (3, 1),
            (4, 1),
            (1, 2),
            (3, 2),
            (4, 2),
            (1, 3),
            (2, 3),
            (4, 3),
            (1, 4),
            (2, 4),
            (3, 4),
      ]
      q[[CartesianIndex(x...) for x in indices]] = rates
      dia = sum(q, dims = 2)
      q[diagind(4, 4)] = -dia
      eq = svd(q).U[:,end]
      eq /= sum(eq)
end

##

function cor4(eq)
      cv = eq[1]*eq[4] - eq[2]*eq[3]
      vx = (eq[1]+eq[2])*(eq[3]+eq[4])
      vy = (eq[1]+eq[3])*(eq[2]+eq[4])
      cv/sqrt(vx*vy)
end

##
nrc = findall(occursin.("NRc",fpairs))

##
eqs_ = []

for charNum = 1:28
      sim = read("modelFitting/output/universal_$(lpad(charNum,2,"0")).jls", ModelChains)
      s_rates = vcat([sim[:, :s_rates, :].value[:, :, k] for k = 1:2]...)
      posteriorEq = mapslices(getEqulibirium, exp.(s_rates), dims = 2)
      if charNum ∈ nrc
            posteriorEq = posteriorEq[:, [2, 1, 4, 3]]
      end
      posteriorCor = mapslices(cor4, posteriorEq, dims = 2)[:, 1]
      push!(eqs_, mapslices(mean, posteriorEq, dims = 1)[1, :])
end

eqs = DataFrame(Tuple.(eqs_))
rename!(eqs, ["00", "01", "10", "11"])
insertcols!(eqs, 1, :fpair => fpairs)

##

cors_ = []

for charNum = 1:28
      sim = read("modelFitting/output/universal_$(lpad(charNum,2,"0")).jls", ModelChains)
      s_rates = vcat([sim[:, :s_rates, :].value[:, :, k] for k = 1:2]...)
      posteriorEq = mapslices(getEqulibirium, exp.(s_rates), dims = 2)
      if charNum ∈ nrc
            posteriorEq = posteriorEq[:, [2, 1, 4, 3]]
      end
      posteriorCor = mapslices(cor4, posteriorEq, dims = 2)[:, 1]
      push!(
            cors_,
            vcat(
                  hpd(posteriorCor, alpha=0.05),
                  hpd(posteriorCor, alpha = 0.5),
                  [median(posteriorCor)],
            ),
      )
end

cors = DataFrame(Tuple.(cors_))
rename!(cors, [:hpd95_l, :hpd95_u, :hpd50_l, :hpd50_u, :median])
insertcols!(cors, 1, :fpair => fpairs)

sort!(cors, :median)

##

theme(:solarized_light)
#
upscale = 2 #8x upscaling in resolution
fntsm = Plots.font("sans-serif", pointsize=round(15.0*upscale))
fntlg = Plots.font("sans-serif", pointsize=round(20.0*upscale))
default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
default(size=(800*upscale,1000*upscale)) #Plot canvas size
sl_palette = PlotThemes._themes[:solarized_light].defaults[:palette]

##

linreg_ = []

for (g,charNum) in enumerate(indexin(cors.fpair, fpairs))
      charName = fpairs[charNum]
      features = string.(split(charName, "-"))
      xNames = [features[1][1:1], features[1][2:end]]
      yNames = [features[2][1:1], features[2][2:end]]
      sim = read("modelFitting/output/universal_$(lpad(charNum,2,"0")).jls", ModelChains)
      s_rates = vcat([sim[:, :s_rates, :].value[:, :, k] for k = 1:2]...)
      posteriorEq = mapslices(getEqulibirium, exp.(s_rates), dims = 2)
      if charNum ∈ nrc
            posteriorEq = posteriorEq[:, [2, 1, 4, 3]]
      end
      ms = 150 * upscale
      p = plot(
            aspectratio=:equal,
            legend = false,
            xlim = [-0.5, 1.5],
            ylim = [-0.5, 1.5],
            xticks = ([0, 1], [xNames[2] * xNames[1], xNames[1] * xNames[2]]),
            yticks = ([0, 1], [yNames[2] * yNames[1], yNames[1] * yNames[2]]),
            left_margin = upscale * 40mm,
            right_margin = upscale * 20mm,
            size=(4*500, 500*7),
      )
#      title!(fpairs[charNum])
      for g in 1:100
            i = rand(1:size(posteriorEq, 1))
            eqs = posteriorEq[i, :]
            scatter!(
                  [0],
                  [0],
                  markersize = ms * sqrt(eqs[1]),
                  alpha = 0.05,
                  markerstrokewidth = 0,
                  color = sl_palette.colors[1],
            )
            scatter!(
                  [0],
                  [1],
                  markersize = ms * sqrt(eqs[2]),
                  alpha = 0.05,
                  markerstrokewidth = 0,
                  color = sl_palette.colors[2],
            )
            scatter!(
                  [1],
                  [0],
                  markersize = ms * sqrt(eqs[3]),
                  alpha = 0.05,
                  markerstrokewidth = 0,
                  color = sl_palette.colors[3],
            )
            scatter!(
                  [1],
                  [1],
                  markersize = ms * sqrt(eqs[4]),
                  alpha = 0.05,
                  markerstrokewidth = 0,
                  color = sl_palette.colors[4],
            )
            a = eqs[2]/(eqs[1]+eqs[2])
            b = eqs[4]/(eqs[3]+eqs[4]) - a
            Plots.abline!(b, a, color=:gray, linewidth=upscale=0.5)
            end
      push!(linreg_, p)
end
plot(linreg_...,layout=grid(7,4, heights=repeat([1/7],28)))

try
      mkdir("../img/")
catch e
end

savefig("../img/linearRegression.pdf")
##

sort!(cors, :median)
p = plot(legend = false, left_margin = upscale*10mm)
for charNum = 1:28
      plot!(
            p,
            Array(cors[charNum, 2:3]),
            [charNum, charNum],
            linewidth = 3 * upscale,
            color = :black,
      )
      plot!(
            p,
            Array(cors[charNum, 4:5]),
            [charNum, charNum],
            linewidth = 10 * upscale,
            color = :black,
      )
      scatter!(
            p,
            [cors[charNum, 6]],
            [charNum],
            markersize = 5 * upscale,
            color = :white,
      )
end

yticks!(p, (1:28, cors.fpair))
vline!(p,[0])
xlabel!("correlation (posterior distribution)")

display(p)

savefig("../img/correlations.pdf")

##

sort!(cors, :median, rev=true)

output = DataFrame(
      "feature pair" => cors.fpair,
      "median" => format.("{1:.2f}",cors.median),
      "HPD" => mapslices(
            x -> "(" * join(x, ", ") * ")",
            format.("{1:.2f}", Array(cors[:, 2:3])),
            dims = 2,
      )[:,1],
)

try
      mkdir("../data/tables/")
catch e
end

CSV.write("../data/tables/correlations.tex", output, newline="\\\\\n", delim="&")

##

prior = rand(Normal(0,1), (100000, 12))

priorEqs = mapslices(getEqulibirium, exp.(prior), dims=2)

priorCor = mapslices(cor4, priorEqs, dims=2)[:,1]

##
bf_ = []

for charNum = 1:28
      sim = read(
            "modelFitting/output/universal_$(lpad(charNum, 2, "0")).jls",
            ModelChains,
      )
      s_rates = vcat([sim[:, :s_rates, :].value[:, :, k] for k = 1:2]...)
      posteriorEq = mapslices(getEqulibirium, exp.(s_rates), dims = 2)
      if charNum ∈ nrc
            posteriorEq = posteriorEq[:, [2, 1, 4, 3]]
      end
      posteriorCor = mapslices(cor4, posteriorEq, dims = 2)[:, 1]
      bfc = convert(
            Float64,
            R"library('LRO.utilities')
            log(savage_dickey($posteriorCor, $priorCor, Q=0, plot=F)$BF01)",
      )
      push!(bf_, bfc)
end
##

bf = DataFrame(:fpair => fpairs, :bf => bf_)

sort!(bf, :bf, rev = true)

p0 = invlogit.(bf.bf)

bf[:, :cum_pv] = [minimum([1.,x]) for x in cumsum(1 .- p0)]

output = DataFrame(
      "feature pair" => bf.fpair,
      "(log) Bayes Factor" => format.("{1:.2f}",round.(bf.bf, digits = 2)),
      "cumulative posterior probability" =>
            rpad.(round.(bf.cum_pv, sigdigits = 3), 5, "0"),
)

CSV.write("../data/tables/bfCorr.tex", output, delim="&", newline="\\\\\n")

##
upscale = 2

sort!(bf, :bf)

scatter(bf.bf, 1:28,
    legend = false,
    yticks= (1:28,bf.fpair),
    xlab="(log) Bayes factor",
    left_margin = upscale*5mm,
    markersize=upscale*3,
    )

savefig("../img/bayesfactorCorr.pdf")

max_stars_repo_name: gerhardJaeger/phylogeneticTypology

max_issues_repo_name: gerhardJaeger/phylogeneticTypology

max_forks_repo_name: gerhardJaeger/phylogeneticTypology


text:
module ExBroadcast

export @mtb, @tab, @mtab, @lzb

const NTHREADS = Ref(Threads.nthreads())
@inline num_threads() = NTHREADS[]

"""
    ExBroadcast.set_num_threads(n)
Set the threads' num.
"""
function set_num_threads(n::Int)
    NTHREADS[] = max(1, min(n, Threads.nthreads()))
    nothing
end
include("thread_provider.jl")

include("util.jl")
include("macro_kernal\\mtb.jl")
include("macro_kernal\\tab.jl")
include("macro_kernal\\mtab.jl")
include("macro.jl")

# include("experiment\\modified_base.jl")

using Requires
@init @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("other_support\\gpusupport.jl")
@init @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" include("other_support\\staticarray.jl")
@init @require Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" include("other_support\\interp.jl")

end


text:
# Invert for the two glacier which have radar

println("Loading BITEModel.jl")
@time using BITEModel
const BM=BITEModel
import KissMCMC
using Base.Test
import JLD2
import Plots

i = 4
niter = [10^3, 3*10^4, 10^5, 3*10^5][i]
nthin = [1, 10, 50, 100][i]

with_good_radar = [409, 428]
region_svalbard = 7
gid = BM.RGIGlacier.(with_good_radar, region_svalbard)[2]

plotyes = false
verbose = false

gl,gb,pp,pm,pn,pl = BM.init_forward(gid, verbose)
# Inversion
###########
n_1d_theta = 3

run=[:h, :iv, :h_iv][3]
(theta0, logposterior, logposterior1d, logprior, logprior_debug,
 loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
 pmcmc_defaults, fit_target) =
     BM.init_inverse(gb, pp, pm, pn, n_1d_theta=n_1d_theta, run=run,
                     theta0_dict = Dict(:btilde=> zeros(BM.F,n_1d_theta), # NOTE in m/y (error)
                                        :fsl=> zeros(BM.F,n_1d_theta), # error
                                        :temp=> zeros(BM.F,n_1d_theta), # error
                                        :dist_exp=> pm.dist_exp,
                                        :iv_h_exp=> pm.iv_h_exp,
                                        :iv_dist_exp1=> pm.iv_dist_exp1,
                                        #:iv_dist_exp2=> pm.iv_dist_exp2,
                                        )
                     )
theta0_prior = deepcopy(theta0)

# some tests:
print("Time to run forward model:")
@inferred fwdm_fn(theta0.th0);
@time hs1d, ivs1d, hs2d, ivs2d, pm_, taus1d, taus1d_l, taus2d = fwdm_fn(theta0.th0);

print("Time to calculate posterior:")
@inferred logposterior(theta0.th0)
@time logposterior(theta0.th0)[1];

print("Prior value:")
@show logprior(theta0.th0)
print("Likelihood value:")
@show loglikelihood(theta0.th0)[1]
print("Posterior value:")
@show logposterior(theta0.th0)[1]

# import NLopt

# # using NLopt
# #
# # Much faster than Optim.jl. LN_SBPLX seems better than LN_NELDERMEAD

# ic = copy(theta0.th0)

# od_ll = (th, grad) -> (logprior(th)[1]>-Inf ? loglikelihood(th)[1] : -Inf)
# od_post = (th, grad) -> logposterior(th)[1]

# # https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#local-derivative-free-optimization
# alg = [:LN_COBYLA, :LN_NELDERMEAD, :LN_BOBYQA, :LN_PRAXIS, :LN_SBPLX][end]
# opt = NLopt.Opt(alg, length(ic))

# NLopt.max_objective!(opt, od_post)
# tol = 1e-6
# NLopt.ftol_rel!(opt, tol)
# NLopt.ftol_abs!(opt, tol)
# NLopt.xtol_abs!(opt, tol)
# NLopt.xtol_rel!(opt, tol)
# NLopt.maxeval!(opt, 5000)

# @time maxf, theta_start, ret = NLopt.optimize(opt, ic)
# good_returns = [:SUCCESS, :STOPVAL_REACHED, :FTOL_REACHED, :XTOL_REACHED]
# if !(ret in good_returns)
#     println("Optimization not successful: $ret\n   Best max found=$maxf)")
# else
#     println("Optimization successful.\n     Max found=$maxf")
# end
# theta0.th0[:] = theta_start

#############
# MCMC
############
pmcmc = BM.MCMCNum(niter = niter,
                   nthin = nthin,
                   sigma_ppdf = ones(length(theta0.th0))*3e-2,
                   alg=[:emcee, :metro][1]) # metro needs many more iterations:
                                            # - emcee good with 10^5
                                            # - metro needs 10^?

# sample posterior
thetas, theta_expect, theta_mode, accept_ratio, blobs, lpost, mc_hs1d, mc_ivs1d,
  mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d = mcmc(logposterior, theta0, pmcmc; verbose=true)

# sample the prior to get that distribution too:
pmcmc_prior = BM.MCMCNum(niter=5*10^5, nthin=1000, ball_radius=0.45, drop_low_accept_ratio=false)
thetas_prior, theta_expect_prior, theta_mode_prior, accept_ratio_prior, _, _, mc_hs1d_prior, mc_ivs1d_prior =
    mcmc(logprior, theta0_prior, pmcmc_prior; verbose=true)

varnames = BM.get_varnames(theta0)
KissMCMC.print_results(thetas, accept_ratio, names=varnames)
KissMCMC.print_results(mc_hs1d, 0.0, names=["h" for i=1:size(mc_hs1d,1)])

JLD2.@save("out-$(BM.getrgi_nr(gid)).jld2",
                  gid, pmcmc, pm, pn, run, n_1d_theta,
           #theta0, gb, don't work because of functions,
           # blobs, # does't work in JLD2
                  thetas, theta_expect, accept_ratio, mc_hs1d, mc_ivs1d,
                  mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d,
                  thetas_prior, theta_expect_prior, accept_ratio_prior)

JLD.@save("out-$(BM.getrgi_nr(gid)).jld",
                  # gid, pm, pn, run, n_1d_theta, #pmcmc does not work in JLD
           #theta0, gb, don't work because of functions,
           blobs,
                  # thetas, theta_expect, accept_ratio, mc_hs1d, mc_ivs1d,
                  # mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d,
          # thetas_prior, theta_expect_prior, accept_ratio_prior
          )

# JLD2.@load("out-409.jld2",
#            gb, pmcmc, pp, pm, pn, run,
#            theta0,
#            thetas, theta_expect, accept_ratio, blobs, mc_hs1d, mc_ivs1d,
#            mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d,
#            thetas_prior, theta_expect_prior, accept_ratio_prior)

pyplot()
fig = BM.plottheta(thetas, theta0, toplot=:all)
Plots.savefig("thetas-$(BM.getrgi_nr(gid)).png", fig)

if plotyes

    display(BM.plotinv1d(gb, blobs, reuse=false))
    display(BM.plotinv1d_iverr(gb, blobs, reuse=false))

    display(BM.plotinv2d(gl, blobs, reuse=false))
    display(BM.plotinv2d_h(gl, blobs, reuse=false))
    display(BM.plotinv2d_iverr(gl, blobs, reuse=false))

    display(BM.plottheta(thetas, theta0, toplot=:all))
    display(BM.plottheta(thetas_prior, theta0, toplot=:all))
    # display(BM.plottheta(thetas, theta0, toplot=:btilde, reuse=false))
    # display(BM.plottheta(thetas, theta0, toplot=:fsl, reuse=false))
    # display(BM.plottheta(thetas, theta0, toplot=:temp, reuse=false))

    display(BM.plottheta_violin((thetas, thetas_prior), theta0, :btilde, reuse=false, width=1))
    display(BM.plottheta_violin((thetas, thetas_prior), theta0, :fsl, reuse=false, width=1))
    display(BM.plottheta_violin((thetas, thetas_prior), theta0, :temp, reuse=false, width=1))

    display(BM.plot2d_h(gl, mh2d))

end

nothing


text:
using MultiFloats: Float64x3
using StaticArrays: SVector
using UUIDs: UUID

push!(LOAD_PATH, @__DIR__)
using PCREO
using PointGroups

@assert isdir(ENV["PCREO_DATABASE_DIRECTORY"])

lsdir(path...) = filter(!startswith('.'),
    readdir(joinpath(path...); sort=false))

function build_pcreo_cache(n_min::Int, n_max::Int)
    result = Dict{Int,Vector{Tuple{Float64x3,UUID}}}()
    println("Building PCREO cache...")
    flush(stdout)
    for num_dir in lsdir(ENV["PCREO_DATABASE_DIRECTORY"])
        num_points = parse(Int, num_dir)
        @assert !haskey(result, num_points)
        if n_min <= num_points <= n_max
            println("Scanning directory $num_dir...")
            flush(stdout)
            entry_list = Tuple{Float64x3,UUID}[]
            for entry_name in lsdir(ENV["PCREO_DATABASE_DIRECTORY"], num_dir)
                reference_name = "PCREO-03-$num_dir-$entry_name.csv"
                reference = PCREORecord(joinpath(
                    ENV["PCREO_DATABASE_DIRECTORY"], num_dir,
                    entry_name, reference_name))
                push!(entry_list, (reference.energy, UUID(entry_name)))
            end
            sort!(entry_list)
            result[num_points] = entry_list
        end
    end
    println("Finished building PCREO cache.")
    flush(stdout)
    return result
end

const PCREO_CACHE = build_pcreo_cache(parse(Int, ARGS[1]), parse(Int, ARGS[2]))

const UUID_REGEX = Regex("[0-9a-fA-F]{8}-[0-9a-fA-F]{4}-" *
                         "[0-9a-fA-F]{4}-[0-9a-fA-F]{4}-[0-9a-fA-F]{12}")

function to_point_vector(points::AbstractMatrix{T}) where {T}
    dimension, num_points = size(points)
    @assert dimension == 3
    return [SVector{3,T}(view(points, :, i)) for i = 1 : num_points]
end

function add_to_database!(filepath::String)
    filename = basename(filepath)
    uuid = UUID(match(UUID_REGEX, filename).match)
    record = PCREORecord(filepath)
    dim_str = lpad(record.dimension, 2, '0')
    num_str = lpad(record.num_points, 8, '0')
    num_dir = joinpath(ENV["PCREO_DATABASE_DIRECTORY"], num_str)
    if !isdir(num_dir)
        mkdir(num_dir)
    end
    if !haskey(PCREO_CACHE, record.num_points)
        PCREO_CACHE[record.num_points] = Tuple{Float64x3,UUID}[]
    end
    for (reference_energy, reference_uuid) in PCREO_CACHE[record.num_points]
        if isapprox(record.energy, reference_energy, rtol=Float64x3(2^-80))
            entry_dir = joinpath(num_dir, string(reference_uuid))
            reference = PCREORecord(joinpath(entry_dir,
                "PCREO-$dim_str-$num_str-$reference_uuid.csv"))
            @assert reference.energy == reference_energy
            if isometric(to_point_vector(record.points),
                         to_point_vector(reference.points), 2^-40)
                newpath = joinpath(entry_dir, filename)
                mv(filepath, newpath)
                return newpath
            end
        end
    end
    push!(PCREO_CACHE[record.num_points], (record.energy, uuid))
    sort!(PCREO_CACHE[record.num_points])
    new_dir = joinpath(num_dir, string(uuid))
    mkdir(new_dir)
    newpath = joinpath(new_dir, filename)
    mv(filepath, newpath)
    return newpath
end

function in_range(n_min, n_max)
    return filename -> n_min <= parse(Int, filename[10:17]) <= n_max
end

function main(n_min::Int, n_max::Int)
    filenames = filter(in_range(n_min, n_max),
                filter(startswith("PCREO-03-"),
                readdir("/Users/dzhang314/pcreo-output-rice")))
    for filename in filenames
        src = joinpath("/Users/dzhang314/pcreo-output-rice", filename)
        try
            dst = add_to_database!(src)
            println(src, " => ", dst)
            flush(stdout)
        catch e
            if e isa AssertionError
                println(filename, ": ", e)
                flush(stdout)
            else
                rethrow(e)
            end
        end
    end
end

main(parse(Int, ARGS[1]), parse(Int, ARGS[2]))


text:
# script to write data to mat file
#
# The script is called as follows
#
# genmodelsims.jl [model] [objective] [ratid] [extendedstats]
#
# to generate the simulated behavior in the mat files. [ratid] and
# [extendedstats] are optional. [ratid] defaults to 0 if not given.
# [extendedstats] defaults to "false" if not given.
#
# [extendedstats] is either "true" or "false" if given. If true, it adds
# additional statistics to the data file. This only works for a (small) subset
# of models, and not for the interleaved condition.
#
# If [ratid] is given, the simulated behavior is written to
# ../fitdata/modelsims/[model]_[objective]_[categ/ident]_rat[ratid].mat
# Otherwise, it is written to
# ../fitdata/modelsims/[model]_[objective]_[categ/ident].mat

using Random, Printf, MAT, CSV

include("common.jl")
include("simmodel.jl")
include("behaviorstats.jl")
include("models.jl")

const nrats = 1
const fitpath = "modelsims"
const datapath = "../fitdata/$(fitpath)"

# returns model for given name
function createmodel(modelname)
    modeldict = Dict("opt" => OptimalModel,
                     "optlapse" => OptimalLapseModel,
                     "optcollapse" => OptimalColLapseModel,
                     "optbiaslapse" => OptimalBiasLapseModel,
                     "optcolbiaslapse" => OptimalColBiasLapseModel,
                     "adflapse" => ADFLapseModel,
                     "adfcollapse" => ADFColLapseModel,
                     "adfbiaslapse" => ADFBiasLapseModel,
                     "adfcolbiaslapse" => ADFColBiasLapseModel,
                     "randomcolbiaslapse" => RandomColBiasLapseModel,
                     "randomcollapse" => RandomColLapseModel,
                     "gammalapse" => GammaLapseModel,
                     "gammacollapse" => GammaColLapseModel,
                     "gammabiaslapse" => GammaBiasLapseModel,
                     "gammacolbiaslapse" => GammaColBiasLapseModel,
                     "deltalapse" => DeltaLapseModel,
                     "deltacollapse" => DeltaColLapseModel,
                     "deltabiaslapse" => DeltaBiasLapseModel,
                     "deltacolbiaslapse" => DeltaColBiasLapseModel,
                     "deltaexpcolbiaslapse" => DeltaExpColBiasLapseModel,
                     "lcalapse" => LCALapseModel,
                     "lcacollapse" => LCAColLapseModel)
    if modelname in keys(modeldict)
        return modeldict[modelname]()
    else
        error("Uknown model identifier $modelname")
    end
end

# process arguments
if !(length(ARGS) ∈ 2:4)
    error("Script needs to be called with two to four argument")
end
const ratid = length(ARGS) ≥ 3 ? parse(Int64, ARGS[3]) : 0
const useratid = ratid > 0
const extendedstats = length(ARGS) == 4 && lowercase(ARGS[4]) == "true"
const modelname = lowercase(ARGS[1])
const m = createmodel(modelname)
const ϕns = ϕnames(m)
const ϕn = length(ϕns)
const objfnname = lowercase(ARGS[2])
const isitlvd = occursin("itlvd", objfnname)
const basefilename = (useratid ? "$(modelname)_$(objfnname)_rat$(ratid)"
                               : "$(modelname)_$(objfnname)")
const fitfilename = "$(datapath)/fit_$(basefilename).csv"
if !isfile(fitfilename)
    error("Fitting data $fitfilename not found")
end
println("Loading fitting data from $fitfilename")
const d = CSV.read(fitfilename)
(maxllh, maxi) = findmax(d[:llh])
if maxi != size(d,1) && maxllh > d[end,end]
    warn("Maximum llh $maxllh found in line $maxi instead of last line")
end
const ϕ = Vector(d[maxi, ϕns])
#ϕ[1] *= 3.0  # raise bound 3-fold
println("Simulating model $modelname with parameters")
println(foldl((a, b) -> a * ", " * b,
        [@sprintf("%s=%.2f", ϕns[i], ϕ[i]) for i in 1:ϕn]))

# helper functions to convert vectors to ones useable for matwrite
dfarr2mat(X::Array{Float64,1}) = X
dfarr2mat(X::Array{Int64,1}) = X
dfarr2mat(X::Array{Bool,1}) = X
dfarr2mat(X::BitArray{1}) = convert(Vector{Bool}, X)

if isitlvd
    # interleaved trials only
    if extendedstats
        error("Extended statistics not supported for the inverleaved condition")
    end

    Random.seed!(0)
    const evidence, ids = catitlvdevidence(itlvd_cs, 10, 10^5)
    const m_itlvd = simitlvdmodel(m, ϕ, evidence, ids)

    # write simulations to file (the below failes with non-ascii column names)
    d_itlvd = Dict{String,Any}([string(s) => dfarr2mat(m_itlvd[s]) for s in names(m_itlvd)])
    matwrite("$(datapath)/$(basefilename)_itlvd.mat", d_itlvd)
    println("Data written to $(datapath)/$(basefilename)_itlvd.mat")
else
    # separate categ/ident trials
    Random.seed!(0)
    const evidence, ids = catevidence(ident_cs, categ_cs, 10, 10^5);
    const m_ident, m_categ = (extendedstats ?
        simmodelwithstats(m, ϕ, evidence, ids) : simmodel(m, ϕ, evidence, ids))

    # write simulations to file (the below failes with non-ascii column names)
    d_ident = Dict{String,Any}([string(s) => dfarr2mat(m_ident[s]) for s in names(m_ident)])
    d_categ = Dict{String,Any}([string(s) => dfarr2mat(m_categ[s]) for s in names(m_categ)])
    matwrite("$(datapath)/$(basefilename)_ident.mat", d_ident)
    println("Data written to $(datapath)/$(basefilename)_ident.mat")
    matwrite("$(datapath)/$(basefilename)_categ.mat", d_categ)
    println("Data written to $(datapath)/$(basefilename)_categ.mat")
end


text:
function typestable(@nospecialize(f), @nospecialize(t); checkonlyany=false)
    v = code_typed(f, t)
    stable = true
    for vi in v
        for (name, ty) in zip(vi[1].slotnames, vi[1].slottypes)
            !(ty isa Type) && continue
            if (checkonlyany && ty === Any) ||
                (!checkonlyany && (!Base.isdispatchelem(ty) || ty == Core.Box))
                stable = false
                println("Type instability is detected! the variable is $(name) ::$ty")
            end
        end
    end
    return stable
end


text:
export rule

@rule AR(:x, Marginalisation) (m_y::NormalDistributionsFamily, q_θ::NormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta) = begin
    mθ, Vθ = mean(q_θ), cov(q_θ)
    my, Vy = mean(m_y), cov(m_y)

    mγ = mean(q_γ)

    mA = as_companion_matrix(mθ)
    mV = ar_transition(getvform(meta), getorder(meta), mγ)

    D = mA'*inv(Vy + mV)*mA + mγ*Vθ
    Vx = inv(D)
    mx = inv(D)*mA'*inv(Vy + mV)*my

    return convert(promote_variate_type(getvform(meta), NormalMeanVariance), mx, Vx)
end


text:
using Distributions
using QuadGK
using Plots
using HCubature

unicodeplots()

function gaussian(x; σ=1, μ=0)
   y = @. exp(-(x-μ)^2/2)/sqrt(2π)/σ
   return y
end

function transform(f)
    function g(t)
        return f(t/(1-t^2))*(1+t^2)/(1-t^2)/(1-t^2)
    end
    return g
end

function integral(f, a, b)
    if a == -Inf && b == Inf
        println("Use inf case")
        g = transform(f)
        return quadgk(g, -1, 1)
    else
        return quadgk(f, a, b)
    end
end

function plotgaussian()
    xs = -3:0.01:3
    ys = gaussian(xs)
    plot(xs, ys)
end

function integraloverinf()
    integral(gaussian, -Inf, Inf)
end

struct Parameter
    a
    b
    Parameter(; a,b)=new(a,b)
end

function p(x; w)
    t -> (1-w.a)*gaussian(t) + a * gaussian(x-w.b)
end

function post_distribution(p, φ, w, X; beta=1.)
    f(w)=x->φ(w)*prod(p(X[i], w=Parameter(a=w[1], b=w[2]))^β for i in 1:length(X))
    res=HCubature(f,[0,-5],[1,5])
end

function main()
x    w = Parameter(0.5, 3.0)
    p(x, w=w)
end

#main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Exemplo
using Test

#import Exemplo
#Exemplo.minha_funcao #não funcionou, mas acho que não é necessário

@testset "Exemplo.jl" begin
    @test soma(2,3) == 5
    @test soma(9,1) == 10
    @test multiplicação(2,3) == 6
end


text:
include("mesh.jl")

function findPolyMeshNData(fileLines)
    lineCounter = 1
    nData = 0
    startLine = 0
    for line in fileLines
        if isNumber(strip(line))
            nData = parse(Int64, line)
            startLine = lineCounter + 2
            break
        end
        lineCounter += 1
    end
    return nData, startLine
end

function isNumber(str)
    re = r"^[+-]?([0-9]+([.][0-9]*)?|[.][0-9]+)$"
    return occursin(re, str)
end

function findInLines(str, lines, startLine)
    nLines = size(lines, 1)
    for i in startLine:nLines
        if occursin(str, lines[i])
            return i
        end
    end
    return -1
end

function readPolyMeshPoints(filePath)
    f = open(filePath, "r")
    fileLines = readlines(f)
    nPoints, startLine = findPolyMeshNData(fileLines)

    points = Vector{Float64}[]
    for iLine in startLine:startLine + nPoints - 1
        noWhiteSpace = strip(fileLines[iLine])
        coordLine = split(noWhiteSpace[2:end - 1])
        coord = Float64[]
        for i in 1:3
            push!(coord, parse(Float64, coordLine[i]))
        end
        push!(points, coord)
    end
    return points
end

function readPolyMeshFaces(filePath)
    f = open(filePath, "r")
    fileLines = readlines(f)
    nFaces, startLine = findPolyMeshNData(fileLines)

    faces = Vector{Int64}[]
    for iLine in startLine:startLine + nFaces - 1
        noWhiteSpace = strip(fileLines[iLine])
        faceLine = split(noWhiteSpace[3:end - 1])
        facePoints = Int64[]
        for iPoint in 1:length(faceLine)
            push!(facePoints, parse(Int64, faceLine[iPoint]) + 1)
        end
        push!(faces, facePoints)
    end
    return faces
end

function readPolyMeshOwner(filePath)
    f = open(filePath, "r")
    fileLines = readlines(f)
    nOwner, startLine = findPolyMeshNData(fileLines)

    owner = Int64[]
    for iLine in startLine:startLine + nOwner - 1
        noWhiteSpace = strip(fileLines[iLine])
        ownerLine = noWhiteSpace
        push!(owner, parse(Int64, ownerLine) + 1)
    end
    return owner
end

function readPolyMeshNeighbour(filePath)
    f = open(filePath, "r")
    fileLines = readlines(f)
    nNeighbour, startLine = findPolyMeshNData(fileLines)

    neighbour = Int64[]
    for iLine in startLine:startLine + nNeighbour - 1
        noWhiteSpace = strip(fileLines[iLine])
        neighbourLine = noWhiteSpace
        push!(neighbour, parse(Int64, neighbourLine) + 1)
    end
    return neighbour
end

function readPolyMeshBoundary(filePath)
    f = open(filePath, "r")
    fileLines = readlines(f)
    nBoundaries, startLine = findPolyMeshNData(fileLines)

    boundaries = Dict[]
    for iBoundary in 1:nBoundaries
        boundaryDict = Dict{String, Any}("name" => "", "type" => "", "nFaces" => 0, "startFace" => 0)

        nameLine = findInLines("{", fileLines, startLine) - 1
        boundaryDict["name"] = strip(fileLines[nameLine])

        typeLine = findInLines("type", fileLines, startLine)
        boundaryDict["type"] = split(strip(fileLines[typeLine]))[2][1:end - 1]

        nFacesLine = findInLines("nFaces", fileLines, startLine)
        boundaryDict["nFaces"] = parse(Int64, split(strip(fileLines[nFacesLine]))[2][1:end - 1])

        startFaceLine = findInLines("startFace", fileLines, startLine)
        boundaryDict["startFace"] = parse(Int64, split(strip(fileLines[startFaceLine]))[2][1:end - 1]) + 1

        push!(boundaries, boundaryDict)

        startLine = findInLines("}", fileLines, startLine) + 1
    end
    return boundaries
end

function readOpenFoamMesh(polyMeshPath)
    println("Reading OpenFOAM mesh in $polyMeshPath")
    points = readPolyMeshPoints("$polyMeshPath/points")
    faces = readPolyMeshFaces("$polyMeshPath/faces")
    owner = readPolyMeshOwner("$polyMeshPath/owner")
    neighbour = readPolyMeshNeighbour("$polyMeshPath/neighbour")
    boundary= readPolyMeshBoundary("$polyMeshPath/boundary")

    return processOpenFoamMesh(points, faces, owner, neighbour, boundary)
end


text:
"""
Creates a sampler from distribution(s)
"""
function make_sample_prior(distr::D) where D <: Distributions.Distribution
    return function sample_prior()
        return rand(distr)
    end
end

function make_sample_prior(distr::Array{D, 1}) where D <: Distributions.Distribution
    return function sample_prior()
        return rand.(distr)
    end
end

function make_sample_prior(param::T, distr::D) where {T <: Tuple, D <: UnionAll}
    return function sample_prior()
        return rand(distr(param...))
    end
end

function make_sample_prior(param::Array{T, 1}, distr::D) where {T <: Tuple, D <: UnionAll}
    d = [distr(p...) for p in param]
    return function sample_prior()
        return rand.(d)
    end
end

function make_sample_prior(param::Array{T, 1}, distr::Array{D, 1}) where {T <: Tuple, D <: UnionAll}
    dA = [d(p...) for (p, d) in zip(param, distr)]
    return function sample_prior()
        return rand.(dA)
    end
end

"""
Creates a density from distribution or set of distributions
"""
function make_density_prior(distr::D) where D <: Distributions.Distribution
    return function density_prior(x::Float64)
        return Distributions.pdf(distr, x)
    end
end

function make_density_prior(distr::Array{D, 1}) where D <: Distributions.Distribution
    return function density_prior(x::Array{Float64, 1})
        tot = 0.0
        for (d, x) in zip(distr, x)
            tot += Distributions.logpdf(d, x)
        end
        return exp(tot)
    end
end

"""
Simultaneously generate sampler and density function for distribution and parameters
"""
function make_model_prior(d::D) where D <: Distributions.Distribution
    return make_sample_prior(d), make_density_prior(d)
end

function make_model_prior(d::Array{D, 1}) where D <: Distributions.Distribution
    return make_sample_prior(d), make_density_prior(d)
end

function make_model_prior(param::T, distr::D) where {T <: Tuple, D <: UnionAll}
    d = distr(param...)
    return make_sample_prior(d), make_density_prior(d)
end

function make_model_prior(param::Array{T, 1}, distr::D) where {T <: Tuple, D <: UnionAll}
    dA = [distr(p...) for p in param]
    return make_sample_prior(dA), make_density_prior(dA)
end

function make_model_prior(param::Array{T, 1}, d::Array{D, 1}) where {T <: Tuple, D <: UnionAll}
    dA = [d(p...) for (p, d) in zip(param, distr)]
    return make_sample_prior(dA), make_density_prior(dA)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using Base.Test

@testset "DCTI field" begin

    # scalar field
    a = DCTI(1)
    @test interpolate(a, 0.0) == 1
    @test a[1] == 1
    @test length(a) == length(1)
    @test size(a) == size(1)
    update!(a, 2)
    @test a == 2

    # vector field
    b = DCTI([1,2])
    @test interpolate(b, 0.0) == [1,2]
    @test b[1] == [1,2]
    @test length(b) == length([1,2])
    @test size(b) == size([1,2])
    update!(b, [2,3])
    @test b == [2,3]

    # tensor field
    c = DCTI([1 2; 3 4])
    @test interpolate(c, 0.0) == [1 2; 3 4]
    @test c[1] == [1 2; 3 4]
    @test length(c) == length([1 2; 3 4])
    @test size(c) == size([1 2; 3 4])
    update!(c, [2 3; 4 5])
    @test c == [2 3; 4 5]

end

@testset "DVTI field" begin

    # scalar field
    a = DVTI((1, 2))
    @test a[1] == 1
    @test a[2] == 2
    @test interpolate(a, 0.0) == (1, 2)
    update!(a, (2,3))
    @test a == (2,3)
    @test (2,3) == a

    # vector field
    b = DVTI(([1,2], [2,3]))
    @test b[1] == [1,2]
    @test b[2] == [2,3]
    @test interpolate(b, 0.0) == ([1,2], [2,3])
    update!(b, ([2,3], [4,5]))
    @test b == ([2,3], [4,5])

    # tensor field
    c = DVTI(([1 2; 3 4], [2 3; 4 5]))
    @test c[1] == [1 2; 3 4]
    @test c[2] == [2 3; 4 5]
    @test interpolate(c, 0.0) == ([1 2; 3 4], [2 3; 4 5])
    update!(c, ([2 3; 4 5], [5 6; 7 8]))
    @test c == ([2 3; 4 5], [5 6; 7 8])

    d = DVTI(2, 3)
    @test a == d
end

@testset "DCTV field" begin

    # scalar field
    a = DCTV(0.0 => 0.0, 1.0 => 1.0)
    @test isapprox(interpolate(a, -1.0), 0.0)
    @test isapprox(interpolate(a, 0.0), 0.0)
    @test isapprox(interpolate(a, 0.5), 0.5)
    @test isapprox(interpolate(a, 1.0), 1.0)
    update!(a, 1.0 => 2.0)
    @test isapprox(interpolate(a, 0.5), 1.0)
    update!(a, 2.0 => 1.0)
    @test isapprox(interpolate(a, 1.5), 1.5)

    # vector field
    b = DCTV(0.0 => [1.0, 2.0], 1.0 => [2.0, 3.0])
    @test isapprox(interpolate(b, 0.5), [1.5, 2.5])

    # tensor field
    c = DCTV(0.0 => [1.0 2.0; 3.0 4.0], 1.0 => [2.0 3.0; 4.0 5.0])
    @test isapprox(interpolate(c, 0.5), [1.5 2.5; 3.5 4.5])
end

@testset "DVTV field" begin
    # scalar field
    a = DVTV(0.0 => (0.0, 1.0), 1.0 => (1.0, 0.0))
    update!(a, 2.0 => (2.0, 0.0))
    r = interpolate(a, 0.5)
    @test isapprox(r[1], 0.5)
    @test isapprox(r[2], 0.5)
    update!(a, 2.0 => (4.0, 0.0))
end

@testset "CVTV field" begin
    f = CVTV((xi,t) -> xi[1]*xi[2]*t)
    @test isapprox(f([1.0,2.0],3.0), 6.0)
end

@testset "Dictionary fields" begin
    X = Dict(1=>[0.0,0.0], 1000=>[1.0,0.0], 100000=>[1.0,1.0])
    G = DVTId(X)
    @test isapprox(G[1], X[1])
    @test isapprox(G[1000], X[1000])
    @test isapprox(G[100000], X[100000])
    Y = Dict(1=>[2.0,2.0], 1000=>[3.0,2.0], 100000=>[3.0,3.0])
    F = DVTVd(0.0 => X, 1.0 => Y)
    @test isapprox(interpolate(F, 0.5)[100000], [2.0,2.0])
end

@testset "update dictionary field" begin
    f1 = Dict(1=>1.0, 2=>2.0, 3=>3.0)
    f2 = Dict(1=>2.0, 2=>3.0, 3=>4.0)
    fld = DVTVd(0.0 => f1)
    update!(fld, 1.0 => f2)
    @test isapprox(interpolate(fld, 0.5)[1], 1.5)
    update!(fld, 1.0 => f1)
    @test isapprox(interpolate(fld, 0.5)[1], 1.0)
end

@testset "use of common constructor field" begin
    @test isa(field(1.0), DCTI)
    @test isa(field(1.0 => 1.0), DCTV)
    @test isa(field((1.0,2.0)), DVTI)
    @test isa(field(1, 2), DVTI)
    @test isa(field(1.0 => (1.0,2.0)), DVTV)
    @test isa(field((xi,t) -> xi[1]*t), CVTV)
    @test isa(field(1 => [1.0, 2.0], 10 => [2.0, 3.0]), DVTId)
    @test isa(field(0.0 => (1=>1.0,10=>2.0), 1.0 => (1=>2.0,10=>3.0)), DVTVd)
    X = Dict(1 => [0.0,0.0], 2 => [1.0,0.0])
    X1 = field(X)
    X2 = field(0.0 => X)
    @test isa(X1, DVTId)
    @test isa(X2, DVTVd)
end

@testset "general interpolation" begin
    a = [1, 2, 3]
    b = (2, 3, 4)
    @test interpolate(a, b) == 2+6+12
    a = (1, 2)
    b = (2, 3, 4)
    @test interpolate(a, b) == 2+6
    @test_throws AssertionError interpolate(b, a)
end