text:
using ReinforcementLearning
using Flux
using CUDA
using Random, Distributions
using Setfield: @set
using Zygote: ignore
using LinearAlgebra

export DDPGPolicyDemo, update_critic!

#This is a modified version of JuliaReinforcementLearning's ddpg policy. It has been adapted to handle all the various training scenarios used during the project.

mutable struct DDPGPolicyDemo{
    BA<:NeuralNetworkApproximator,
    BC<:NeuralNetworkApproximator,
    TA<:NeuralNetworkApproximator,
    TC<:NeuralNetworkApproximator,
    P,
    R<:AbstractRNG,
} <: AbstractPolicy

    behavior_actor::BA
    behavior_critic::BC
    target_actor::TA
    target_critic::TC
    γ::Float32
    ρ::Float32
    batch_size::Int
    start_steps::Int
    start_policy::P
    update_after::Int
    update_every::Int
    act_limit::Float64
    act_noise::Float64
    step::Int
    rng::R
    # for logging
    actor_loss::Float32
    critic_loss::Float32
end

Flux.functor(x::DDPGPolicyDemo) = (
    ba = x.behavior_actor,
    bc = x.behavior_critic,
    ta = x.target_actor,
    tc = x.target_critic,
),
y -> begin
    x = @set x.behavior_actor = y.ba
    x = @set x.behavior_critic = y.bc
    x = @set x.target_actor = y.ta
    x = @set x.target_critic = y.tc
    x
end

"""
    DDPGPolicy(;kwargs...)
# Keyword arguments
- `behavior_actor`,
- `behavior_critic`,
- `target_actor`,
- `target_critic`,
- `start_policy`,
- `γ = 0.99f0`,
- `ρ = 0.995f0`,
- `batch_size = 32`,
- `start_steps = 10000`,
- `update_after = 1000`,
- `update_every = 50`,
- `act_limit = 1.0`,
- `act_noise = 0.1`,
- `step = 0`,
- `rng = Random.GLOBAL_RNG`,
"""
function DDPGPolicyDemo(;
    behavior_actor,
    behavior_critic,
    target_actor,
    target_critic,
    start_policy,
    γ = 0.99f0,
    ρ = 0.995f0,
    batch_size = 32,
    start_steps = 10000,
    update_after = 1000,
    update_every = 50,
    act_limit = 1.0,
    act_noise = 0.1,
    step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(behavior_actor, target_actor)  # force sync
    copyto!(behavior_critic, target_critic)  # force sync
    DDPGPolicyDemo(
        behavior_actor,
        behavior_critic,
        target_actor,
        target_critic,
        γ,
        ρ,
        batch_size,
        start_steps,
        start_policy,
        update_after,
        update_every,
        act_limit,
        act_noise,
        step,
        rng,
        0.0f0,
        0.0f0,
    )
end

function (p::DDPGPolicyDemo)(env::AbstractEnv; training = true)
    p.step += 1

    if p.step <= p.start_steps
        p.start_policy(env)
    else
        D = device(p.behavior_actor)
        s = state(env)
        s = Flux.unsqueeze(s, ndims(s) + 1)
        a = p.behavior_actor(send_to_device(D, s)) |> vec |> send_to_host
        clamp(a[] + randn(p.rng) * p.act_noise * training, -p.act_limit, p.act_limit)
    end
end

function action(p::DDPGPolicyDemo, s; training = true)
    p.step += 1

    if p.step <= p.start_steps
        p.start_policy(s)
    else
        D = device(p.behavior_actor)
        s = Flux.unsqueeze(s, ndims(s) + 1)
        a = p.behavior_actor(send_to_device(D, s)) |> vec |> send_to_host
        clamp(a[] + randn(p.rng) * p.act_noise * training, -p.act_limit, p.act_limit)
    end
end

function RLBase.update!(p::DDPGPolicyDemo, traj::CircularArraySARTSTrajectory, ::AbstractEnv, ::PreActStage)
    length(traj) > p.update_after || return
    p.step % p.update_every == 0 || return
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(p.batch_size))
    update!(p, batch)
end

function update_critic!(p::DDPGPolicyDemo, traj::CircularArraySARTSTrajectory, demo_traj::CircularArraySARTSTrajectory)
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(ceil(Int, 3*p.batch_size/4)))
    d_inds, d_batch = sample(p.rng, traj, BatchSampler{SARTS}(ceil(Int, p.batch_size/4)))
    _update_critic(p, batch, d_batch=d_batch)
end

function update_critic!(p::DDPGPolicyDemo, traj::CircularArraySARTSTrajectory)
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(ceil(Int, p.batch_size)))
    _update_critic(p, batch)
end

function RLBase.update!(p::DDPGPolicyDemo, traj::CircularArraySARTSTrajectory)
    length(traj) > p.update_after || return
    p.step % p.update_every == 0 || return
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(p.batch_size))
    update!(p, batch)
end

function RLBase.update!(p::DDPGPolicyDemo, traj::CircularArraySARTSTrajectory, ::AbstractEnv, stage::PreTrainStage)
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(p.batch_size))
    update!(p, batch, stage)
end

function RLBase.update!(p::DDPGPolicyDemo, traj::CircularArraySARTSTrajectory, demo_traj::CircularArraySARTSTrajectory)
    length(traj) > p.update_after || return
    p.step % p.update_every == 0 || return
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(ceil(Int, p.batch_size*0.75)))
    demo_inds, demo_batch = sample(p.rng, demo_traj, BatchSampler{SARTS}(ceil(Int, p.batch_size*0.25)))
    update!(p, batch, demo_batch)
end

function RLBase.update!(p::DDPGPolicyDemo, batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(p), batch)

    A = p.behavior_actor
    C = p.behavior_critic
    Aₜ = p.target_actor
    Cₜ = p.target_critic

    γ = p.γ
    ρ = p.ρ

    a′ = Aₜ(s′)
    qₜ = Cₜ(vcat(s′, a′)) |> vec
    y = r .+ γ .* (1 .- t) .* qₜ
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(C)) do
        q = C(vcat(s, a)) |> vec
        loss = mean((y .- q) .^ 2)
        ignore() do
            p.critic_loss = loss
        end
        loss
    end

    update!(C, gs1)

    gs2 = gradient(Flux.params(A)) do
        loss = -mean(C(vcat(s, A(s))))
        ignore() do
            p.actor_loss = loss
        end
        loss
    end

    update!(A, gs2)

    # polyak averaging
    for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
        dest .= ρ .* dest .+ (1 - ρ) .* src
    end
end

function RLBase.update!(p::DDPGPolicyDemo, batch::NamedTuple{SARTS}, ::PreTrainStage)
    s, a, r, t, s′ = send_to_device(device(p), batch)

    A = p.behavior_actor
    C = p.behavior_critic
    Aₜ = p.target_actor
    Cₜ = p.target_critic

    γ = p.γ
    ρ = p.ρ

    a′ = Aₜ(s′)
    qₜ = Cₜ(vcat(s′, a′)) |> vec
    y = r .+ γ .* (1 .- t) .* qₜ
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(C)) do
        q = C(vcat(s, a)) |> vec
        loss = mean((y .- q) .^ 2) + sum(norm, Flux.params(C))/2
        ignore() do
            p.critic_loss = loss
        end
        loss
    end

    update!(C, gs1)

    gs2 = gradient(Flux.params(A)) do
        loss = -mean(C(vcat(s, A(s)))) + 1/2*mean((A(s) .- a) .^ 2) + sum(norm, Flux.params(A))/4
        ignore() do
            p.actor_loss = loss
        end
        loss
    end

    update!(A, gs2)

    # polyak averaging
    for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
        dest .= ρ .* dest .+ (1 - ρ) .* src
    end
end

function RLBase.update!(p::DDPGPolicyDemo, batch::NamedTuple{SARTS}, demo_batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(p), batch)
    sᵈ, aᵈ, rᵈ, tᵈ, s′ᵈ = send_to_device(device(p), demo_batch)

    A = p.behavior_actor
    C = p.behavior_critic
    Aₜ = p.target_actor
    Cₜ = p.target_critic

    γ = p.γ
    ρ = p.ρ

    a′ = Aₜ(s′)
    qₜ = Cₜ(vcat(s′, a′)) |> vec
    y = r .+ γ .* (1 .- t) .* qₜ
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(C)) do
        q = C(vcat(s, a)) |> vec
        loss = mean((y .- q) .^ 2) + sum(norm, Flux.params(C))/2
        ignore() do
            p.critic_loss = loss
        end
        loss
    end

    update!(C, gs1)

    aᵈ = Flux.unsqueeze(aᵈ, 1)
    Qcomp = Cₜ(vcat(sᵈ, aᵈ)) .> Cₜ(vcat(sᵈ, A(sᵈ)))
    gs2 = gradient(Flux.params(A)) do
        loss = -mean(C(vcat(s, A(s)))) .+ 1/2*mean(((A(sᵈ) .- aᵈ) .^ 2) .* Qcomp) .+ sum(norm, Flux.params(A)) ./ 4
        ignore() do
            p.actor_loss = loss
        end
        loss
    end

    update!(A, gs2)

    # polyak averaging
    for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
        dest .= ρ .* dest .+ (1 - ρ) .* src
    end
end

# gradient gives error: scalar getindex is disallowed if x is a CuArray with Float64, not FLoat32.
# Solved by calling |> gpu at the end, converting CuArray{Float64,n} to CuArray{Float32,n}
function _update_critic(p::DDPGPolicyDemo, batch::NamedTuple{SARTS}; d_batch = nothing)
    s, a, r, t, s′ = send_to_device(device(p), batch)
    if !isnothing(d_batch)
        sᵈ, aᵈ, rᵈ, tᵈ, s′ᵈ = send_to_device(device(p), d_batch)

        s = hcat(s,sᵈ)
        a = vcat(a,aᵈ)
        r = vcat(r,rᵈ)
        t = vcat(t,tᵈ)
        s′ = hcat(s′,s′ᵈ)
    end

    A = p.behavior_actor
    C = p.behavior_critic
    Aₜ = p.target_actor
    Cₜ = p.target_critic

    γ = p.γ
    ρ = p.ρ

    a′ = Aₜ(s′)
    qₜ = Cₜ(vcat(s′, a′)) |> vec
    y = r .+ γ .* (1 .- t) .* qₜ
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(C)) do
        q = C(vcat(s, a)) |> vec
        loss = mean((y .- q) .^ 2)
        ignore() do
            p.critic_loss = loss
        end
        loss
    end

    update!(C, gs1)

    for (dest, src) in zip(Flux.params(Cₜ), Flux.params(C))
        dest .= ρ .* dest .+ (1 - ρ) .* src
    end
end

function _pgd(A, C, x, ϵ)
    x0 = copy(x)
    α = 2.5*ϵ/10
    for k in 1:10
        grad = gradient(Flux.params(x)) do
            -mean(C(vcat(x, A(x))))
        end

        x = x + α*sign.(grad[x])
        x = clamp.(x, x0.-ϵ, x0.+ϵ)
        x = send_to_device(device(A), x)
    end
    x
end

max_stars_repo_name: polhager/FurutaRL

max_issues_repo_name: polhager/FurutaRL

max_forks_repo_name: polhager/FurutaRL


text:
include("params.jl")
include("run_rtaa.jl")
include("run_cmax.jl")
include("run_finite_model_class.jl")
include("run_model_search.jl")
include("run_agnostic_sysid.jl")
include("run_online_model_search.jl")

function mountaincar_all_main()
    rtaa_steps = mountaincar_rtaa_main()
    cmax_steps = mountaincar_cmax_main()
    true_steps = mountaincar_true_main()
    # finite_model_class_steps = mountaincar_finite_model_class_main()
    # true_finite_model_class_steps = mountaincar_finite_model_class_main(true_model = true)
    # local_finite_model_class_steps = mountaincar_finite_model_class_main(local_agent = true)
    return_based_model_search_steps = mountaincar_return_based_model_search_main()
    bellman_based_model_search_steps = mountaincar_bellman_based_model_search_main()

    plot(range_of_values, rtaa_steps, lw = 3, label = "RTAA*", legend = :topleft)
    plot!(range_of_values, cmax_steps, lw = 3, label = "CMAX")
    plot!(range_of_values, true_steps, lw = 3, label = "True")
    # plot!(range_of_values, finite_model_class_steps, lw = 3, label = "Finite Model Class")
    # plot!(
    #     range_of_values,
    #     true_finite_model_class_steps,
    #     lw = 3,
    #     label = "Finite Model Class with true model",
    # )
    # plot!(
    #     range_of_values,
    #     local_finite_model_class_steps,
    #     lw = 3,
    #     label = "Finite Model Class with Local Data",
    # )
    plot!(range_of_values, return_based_model_search_steps, lw = 3, label = "RBMS")
    plot!(range_of_values, bellman_based_model_search_steps, lw = 3, label = "BBMS")
    xlabel!("Misspecification")
    ylabel!("Number of steps to reach goal")
end

ext: jl


text:
@testset "Testing matrix_eqn.jl" begin
    rough_kwargs = Dict(:atol => 1e-7, :rtol => 1e-7)

    @testset "simple test where X is all zero" begin
        A = Matrix(I, 4, 4) .* .95
        B = zeros(4, 4)

        X = solve_discrete_lyapunov(A, B)

        @test X == (B)
    end

    @testset "simple test where X is same as B" begin
        A = fill(0.5, 2, 2)
        B = [.5 -.5; -.5 .5]

        X = solve_discrete_lyapunov(A, B)

        @test X == (B)
    end

    @testset "simple_ones lyap" begin
        A = zeros(4, 4)
        B = ones(4, 4)

        sol = solve_discrete_lyapunov(A, B)

        @test sol == (B)
    end

    @testset "scalars lyap" begin
        a, b = 0.5, 0.75
        sol = solve_discrete_lyapunov(a, b)

        @test sol == ones(1, 1)
    end

    @testset "testing ricatti golden_num_float" begin
        val = solve_discrete_riccati(1.0, 1.0, 1.0, 1.0)
        gold_ratio = (1 + sqrt(5)) / 2.
        @test isapprox(val[1], gold_ratio; rough_kwargs...)
    end

    @testset "testing ricatti golden_num_2d" begin
        If64 = Matrix{Float64}(I, 2, 2)
        A, B, R, Q = If64, If64, If64, If64
        gold_diag = If64 .* (1 + sqrt(5)) ./ 2.
        val = solve_discrete_riccati(A, B, Q, R)
        @test isapprox(val, gold_diag; rough_kwargs...)
    end

    @testset "test tjm 1" begin
        A = [0.0 0.1 0.0
             0.0 0.0 0.1
             0.0 0.0 0.0]
        B = [1.0 0.0
             0.0 0.0
             0.0 1.0]
        Q = [10^5 0.0 0.0
             0.0 10^3 0.0
             0.0 0.0 -10.0]
        R = [0.0 0.0
             0.0 1.0]
        X = solve_discrete_riccati(A, B, Q, R)
        Y = diagm(0 => [1e5, 1e3, 0.0])

        @test isapprox(X, Y; rough_kwargs...)
    end

    @testset "test tjm 2" begin
        A = [0.0 -1.0
             0.0 2.0]
        B = [1.0 0.0
             1.0 1.0]
        Q = [1.0 0.0
             0.0 0.0]
        R = [4.0 2.0
             2.0 1.0]
        X = solve_discrete_riccati(A, B, Q, R)
        Y = zeros(2, 2)
        Y[1, 1] = 1
        @test isapprox(X, Y; rough_kwargs...)
    end

    @testset "test tjm 3" begin
        r = 0.5
        Im = Matrix{Float64}(I, 2, 2)
        A = [2.0+r^2 0.0
             0.0     0.0]
        B = Im
        R = [1.0 r
             r   r*r]
        Q = Im - A' * A + A' * ((R + Im) \ A)
        X = solve_discrete_riccati(A, B, Q, R)
        Y = Im
        @test isapprox(X, Y; rough_kwargs...)
    end

end  # @testset


text:
include("imports.jl")

@testset ExtendedTestSet "getcallback" begin
    learner = testlearner(usedefaultcallbacks = true)
    @test getcallback(learner, Metrics) isa Metrics
    @test getcallback(learner, Checkpointer) isa Nothing
end

@testset ExtendedTestSet "removecallback!" begin
    learner = testlearner(usedefaultcallbacks = true)
    @test removecallback!(learner, Metrics) isa Metrics
    @test getcallback(learner, Metrics) isa Nothing

end

@testset ExtendedTestSet "replacecallback!" begin
    learner = testlearner(usedefaultcallbacks = true)
    cb = getcallback(learner, Metrics)
    cbtaken = replacecallback!(learner, Metrics(accuracy))
    @test cb === cbtaken
    removecallback!(learner, Metrics)
    @test replacecallback!(learner, cb) isa Nothing

end

@testset ExtendedTestSet "addcallback!" begin
    learner = testlearner(usedefaultcallbacks = true)
    cb = Checkpointer(mktempdir())
    @test getcallback(learner, Checkpointer) isa Nothing
    @test_nowarn addcallback!(learner, cb)
    @test getcallback(learner, Checkpointer) isa Checkpointer
end

@testset ExtendedTestSet "setcallbacks!" begin
    learner = testlearner(usedefaultcallbacks = true)
    @test length(learner.callbacks.cbs) != 1
    @test_nowarn setcallbacks!(learner, [Metrics(accuracy)])
    @test length(learner.callbacks.cbs) == 1
end


text:
# taken from: http://www.sunshine2k.de/coding/java/TriangleRasterization/TriangleRasterization.html#algo2

function _draw_hline!(img, xi::Int64, xf::Int64, y::Int64, color)
    y < 1 && return
    y > size(img)[1] && return
    if xi>xf
        xi,xf = xf,xi
    end
    for x in max(xi,1):min(size(img)[2],xf)
        img[y,x,:] .= color
    end
end

function _fill_bottom_flat_triangle!(img, v1x::Int64,v1y::Int64, v2x::Int64,v2y::Int64, v3x::Int64,v3y::Int64, color)
    invslope1::Float64 = (v2x - v1x) / (v2y - v1y);
    invslope2::Float64 = (v3x - v1x) / (v3y - v1y);
    curx1::Float64 = v1x;
    curx2::Float64 = v1x;
    scanlineY::Int64 = v1y;
    while scanlineY <= v2y
        _draw_hline!(img, floor(Int64, curx1), floor(Int64, curx2), scanlineY, color);
        curx1 += invslope1;
        curx2 += invslope2;
        scanlineY += 1
    end
end

function _fill_top_flat_triangle!(img, v1x::Int64,v1y::Int64, v2x::Int64,v2y::Int64, v3x::Int64,v3y::Int64, color)
    invslope1::Float64 = (v3x - v1x) / (v3y - v1y);
    invslope2::Float64 = (v3x - v2x) / (v3y - v2y);
    curx1::Float64 = v3x;
    curx2::Float64 = v3x;
    scanlineY::Int64 = v3y;
    while scanlineY > v1y
        _draw_hline!(img, floor(Int64, curx1), floor(Int64, curx2), scanlineY, color);
        curx1 -= invslope1;
        curx2 -= invslope2;
        scanlineY -= 1
    end
end

function draw_triangle!(img, v1x::Int64,v1y::Int64, v2x::Int64,v2y::Int64, v3x::Int64,v3y::Int64, color)
    # at first sort the three vertices by y-coordinate ascending so v1 is the topmost vertice
    if v1y > v2y
        v1x,v2x = v2x,v1x
        v1y,v2y = v2y,v1y
    end
    if v2y > v3y
        v2x,v3x = v3x,v2x
        v2y,v3y = v3y,v2y
    end
    if v1y > v2y
        v1x,v2x = v2x,v1x
        v1y,v2y = v2y,v1y
    end
    v1y > size(img)[1] && return
    v3y < 1 && return
    max(v1x,v2x,v3x) < 1 && return
    min(v1x,v2x,v3x) > size(img)[2] && return

    # here we know that v1.y <= v2.y <= v3.y
    # check for trivial case of bottom-flat triangle
    if v2y == v3y
        _fill_bottom_flat_triangle!(img, v1x,v1y, v2x,v2y, v3x,v3y, color);
    elseif v1y == v2y
        # check for trivial case of top-flat triangle
        _fill_top_flat_triangle!(img, v1x,v1y, v2x,v2y, v3x,v3y, color);
    else
        # general case - split the triangle in a topflat and bottom-flat one
        v4x::Int64 = floor(Int64, (v1x + ((v2y - v1y) / (v3y - v1y)) * (v3x - v1x)))
        v4y::Int64 = v2y
        _fill_bottom_flat_triangle!(img, v1x,v1y, v2x,v2y, v4x,v4y, color);
        _fill_top_flat_triangle!(img, v2x,v2y, v4x,v4y, v3x,v3y, color);
    end
end

# line from: https://github.com/JuliaImages/ImageDraw.jl/blob/master/src/line2d.jl
function line!(img, x0::Int, y0::Int, x1::Int, y1::Int, color)
    dx = abs(x1 - x0)
    dy = abs(y1 - y0)

    sx = x0 < x1 ? 1 : -1
    sy = y0 < y1 ? 1 : -1;

    err = (dx > dy ? dx : -dy) / 2

    while true
        if y0 >= 1 && y0 <= size(img)[1] && x0 >= 1 && x0 <= size(img)[2]
            img[y0, x0] = color
        end
        (x0 != x1 || y0 != y1) || break
        e2 = err
        if e2 > -dx
            err -= dy
            x0 += sx
        end
        if e2 < dy
            err += dx
            y0 += sy
        end
    end

    img
end


text:
module ARMANets

export ARConv

include("layers.jl")

end


text:
"""
Case 29:
This case study a three bus system with 1 Generic Renewable Model (REPCA, REECB, REGCA) and an infinite source.
The fault drop the connection between buses 1 and 3, eliminating the direct connection between the load at bus 1
and the generator located in bus 3. The infinite generator is located at bus 2.
"""

##################################################
############### LOAD DATA ########################
##################################################

include(joinpath(TEST_FILES_DIR, "data_tests/test29.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

names = ["RENA: No Flags", "RENA: Freq Flag"]
F_flags = [0, 1]

csv_files = (
    joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/TEST_RENA_DEFAULT_FLAGS.csv"),
    joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/TEST_RENA_FREQ_FLAG.csv"),
)

init_conditions = [test29_x0_init, test29_x0_Fflag_init]

eigs_values = [test29_eigvals, test29_eigvals_fflag]

# time span
tspan = (0.0, 5.0);

function test_renA_implicit(csv_file, init_cond, eigs_value, F_Flag)
    path = (joinpath(pwd(), "test-psse-renA"))
    !isdir(path) && mkdir(path)
    try
        sys = System(threebus_file_dir)
        add_source_to_ref(sys)

        for g in get_components(Generator, sys)
            case_gen = inv_generic_renewable(g, F_Flag)
            add_component!(sys, case_gen, g)
        end

        Ybus_change = BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2")

        sim = Simulation(ResidualModel, sys, path, tspan, Ybus_change)

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in init_cond
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2

        # Solve problem
        @test execute!(
            sim,
            IDA(),
            dtmax = 0.005,
            saveat = 0.005,
            abstol = 1e-9,
            reltol = 1e-9,
        ) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for generator
        t, voltage = get_voltage_magnitude_series(results, 103)
        _, power = get_activepower_series(results, "generator-103-1")
        _, rpower = get_reactivepower_series(results, "generator-103-1")

        M = get_csv_data(csv_file)
        M_t = M[:, 1]
        M_p = M[:, 2]
        M_q = M[:, 3]
        M_v = M[:, 4]
        t_psse, p_psse = clean_extra_timestep!(M_t, M_p)
        _, q_psse = clean_extra_timestep!(M_t, M_q)
        _, v_psse = clean_extra_timestep!(M_t, M_v)

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(voltage - v_psse, 2) <= 1e-2
        @test LinearAlgebra.norm(power - p_psse, 2) <= 1e-2
        @test LinearAlgebra.norm(rpower - q_psse, 2) <= 1e-2
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

function test_renA_mass_matrix(csv_file, init_cond, eigs_value, F_Flag)
    path = (joinpath(pwd(), "test-psse-renA"))
    !isdir(path) && mkdir(path)
    try
        sys = System(threebus_file_dir)
        add_source_to_ref(sys)

        for g in get_components(Generator, sys)
            case_gen = inv_generic_renewable(g, F_Flag)
            add_component!(sys, case_gen, g)
        end

        Ybus_change = BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2")

        sim = Simulation(MassMatrixModel, sys, path, tspan, Ybus_change)

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in init_cond
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        # Test Eigenvalues
        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2

        # Solve problem
        @test execute!(
            sim,
            Rodas4(),
            dtmax = 0.005,
            saveat = 0.005,
            abstol = 1e-6,
            reltol = 1e-6,
        ) == PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for generator
        t, voltage = get_voltage_magnitude_series(results, 103)
        _, power = get_activepower_series(results, "generator-103-1")
        _, rpower = get_reactivepower_series(results, "generator-103-1")

        M = get_csv_data(csv_file)
        M_t = M[:, 1]
        M_p = M[:, 2]
        M_q = M[:, 3]
        M_v = M[:, 4]
        t_psse, p_psse = clean_extra_timestep!(M_t, M_p)
        _, q_psse = clean_extra_timestep!(M_t, M_q)
        _, v_psse = clean_extra_timestep!(M_t, M_v)

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(voltage - v_psse, 2) <= 1e-2
        @test LinearAlgebra.norm(power - p_psse, 2) <= 1e-2
        @test LinearAlgebra.norm(rpower - q_psse, 2) <= 1e-2
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 29 RENA ResidualModel" begin
    for (ix, name) in enumerate(names)
        @testset "$(name)" begin
            #dyr_file = dyr_files[ix]
            csv_file = csv_files[ix]
            F_flag = F_flags[ix]
            init_cond = init_conditions[ix]
            eigs_value = eigs_values[ix]
            test_renA_implicit(csv_file, init_cond, eigs_value, F_flag)
        end
    end
end

@testset "Test 29 RENA MassMatrixModel" begin
    for (ix, name) in enumerate(names)
        @testset "$(name)" begin
            #dyr_file = dyr_files[ix]
            csv_file = csv_files[ix]
            F_flag = F_flags[ix]
            init_cond = init_conditions[ix]
            eigs_value = eigs_values[ix]
            test_renA_mass_matrix(csv_file, init_cond, eigs_value, F_flag)
        end
    end
end

"""
Case 29: Test with dyr files
This case study a three bus system with 1 Generic Renewable Model (REPCA, REECB, REGCA) and an infinite source.
The fault drop the connection between buses 1 and 3, eliminating the direct connection between the load at bus 1
and the generator located in bus 3. The infinite generator is located at bus 2.
"""

raw_file_dir = joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/ThreeBusRenewable.raw")

names_dyr = [
    "RENA: No F_Flag, Ref_Flag = 1, V_Flag = 1, Q_Flag = 0",
    "RENA: No F_Flag, Ref_Flag = 1, V_Flag = 0, Q_Flag = 0",
    "RENA: No F_Flag, Ref_Flag = 1, V_Flag = 1, Q_Flag = 1",
]

csv_files_dyr = [
    joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/TEST_RENA_NOFREQ_REF_FLAG.csv"),
    joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/TEST_RENA_NOFREQ_REF_FLAG.csv"),
    joinpath(TEST_FILES_DIR, "benchmarks/psse/RENA/TEST_RENA_NOFREQ_REF_FLAG_Q_FLAG.csv"),
]

dyr_files = [
    joinpath(
        TEST_FILES_DIR,
        "benchmarks/psse/RENA/ThreeBus_REN_A_NOFREQFLAG_with_REF_FLAG.dyr",
    ),
    joinpath(
        TEST_FILES_DIR,
        "benchmarks/psse/RENA/ThreeBus_REN_A_NOFREQFLAG_with_REF_FLAG_no_V_FLAG.dyr",
    ),
    joinpath(
        TEST_FILES_DIR,
        "benchmarks/psse/RENA/ThreeBus_REN_A_NOFREQFLAG_with_REF_FLAG_Q_FLAG.dyr",
    ),
]

init_conditions_dyr =
    [test29_x0_init_Rflag, test29_x0_init_Rflag_no_Vflag, test29_x0_init_Rflag_Qflag]

eigs_values_dyr =
    [test29_eigvals_Rflag, test29_eigvals_Rflag_no_Vflag, test29_eigvals_Rflag_Qflag]

tspans = [(0.0, 10.0), (0.0, 10.0), (0.0, 20.0)]

function test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan)
    path = (joinpath(pwd(), "test-psse-renA_dyr"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file_dir, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            ResidualModel,
            sys, #system
            path,
            tspan, #time span
            BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2"), #Type of Fault
        ) #Type of Fault

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in init_cond
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Test Eigenvalues
        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2

        # Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for generator
        t, voltage = get_voltage_magnitude_series(results, 103)
        _, angl = get_voltage_angle_series(results, 103)
        _, angl_ref = get_voltage_angle_series(results, 102)

        # Obtain data from PSS/E
        M, _ = readdlm(csv_file, ',', header = true)

        M_t = M[:, 1]
        M_θ = M[:, 2]
        M_v = M[:, 3]

        t_psse, v_psse = clean_extra_timestep!(M_t, M_v)
        _, θ_psse = clean_extra_timestep!(M_t, M_θ)

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(voltage - v_psse, 2) <= 1e-2
        @test LinearAlgebra.norm((angl - angl_ref) - θ_psse .* pi / 180, 2) <= 1e-2
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

function test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan)
    path = (joinpath(pwd(), "test-psse-renA_dyr"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file_dir, dyr_file)

        # Define Simulation Problem
        sim = Simulation!(
            MassMatrixModel,
            sys, #system
            path,
            tspan, #time span
            BranchTrip(1.0, Line, "BUS 1-BUS 3-i_2"), #Type of Fault
        ) #Type of Fault

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in init_cond
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Test Eigenvalues
        @test LinearAlgebra.norm(eigs - eigs_value) < 1e-2

        # Solve problem
        @test execute!(sim, Rodas4(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # Obtain data for generator
        t, voltage = get_voltage_magnitude_series(results, 103)
        _, angl = get_voltage_angle_series(results, 103)
        _, angl_ref = get_voltage_angle_series(results, 102)

        # Obtain data from PSS/E
        M, _ = readdlm(csv_file, ',', header = true)

        M_t = M[:, 1]
        M_θ = M[:, 2]
        M_v = M[:, 3]

        t_psse, v_psse = clean_extra_timestep!(M_t, M_v)
        _, θ_psse = clean_extra_timestep!(M_t, M_θ)

        # Test Transient Simulation Results
        # PSSE results are in Degrees
        @test LinearAlgebra.norm(voltage - v_psse, 2) <= 1e-2
        @test LinearAlgebra.norm((angl - angl_ref) - θ_psse .* pi / 180, 2) <= 1e-2
        @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 29 RENA ResidualModel with dyr" begin
    for (ix, name) in enumerate(names_dyr)
        @testset "$(name)" begin
            dyr_file = dyr_files[ix]
            csv_file = csv_files_dyr[ix]
            init_cond = init_conditions_dyr[ix]
            eigs_value = eigs_values_dyr[ix]
            tspan = tspans[ix]
            test_renA_implicit_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan)
        end
    end
end

@testset "Test 29 RENA MassMatrixModel with dyr" begin
    for (ix, name) in enumerate(names_dyr)
        @testset "$(name)" begin
            dyr_file = dyr_files[ix]
            csv_file = csv_files_dyr[ix]
            init_cond = init_conditions_dyr[ix]
            eigs_value = eigs_values_dyr[ix]
            tspan = tspans[ix]
            test_renA_massmatrix_dyr(dyr_file, csv_file, init_cond, eigs_value, tspan)
        end
    end
end


text:
module Example

using Comm, Par, Refs

function main()
    println("Setting up $(nprocs()) matrices...")
    # @remote performs a calculation remotely, and leaves the result
    # on the remote process. The calculation proceeds in the
    # background, while the local process continues.
    As = [@remote p rand(1000, 1000) for p in 1:nprocs()]

    println("Squaring each matrix...")
    # To access a remote object, use get() on the process where the
    # object lives
    Bs = [@remote A get(A)^2 for A in As]

    println("Calculating and collecting matrix norms...")
    # @rpar performs a remote calculation, but returns the result to
    # the caller instead of leaving it on the remote system. The
    # calculation and collecting still proceeds in the background.
    ns = [@rpar B norm(get(B)) for B in Bs]

    println("Outputting norms:")
    # This call to get() implicitly waits for all the previous
    # calculations to finish.
    println([get(n) for n in ns])

    println("Done.")
end

run_main(main)

end


text:
using OffsetArrays
using SpecialFunctions: gamma, erfcx
using PyPlot
using Polynomials
using Printf
using LinearAlgebra

Γ(z) = gamma(z)

function pade_matrix(α::T, β::T,
                     m::Integer, n::Integer) where T <: AbstractFloat
    if mod(m+n,2) == 0
        throw(ArgumentError("m+n must be odd"))
    end
    r = div(m+n-1,2)
    a = OffsetArray{T}(undef, 0:m-1)
    b = OffsetArray{T}(undef, 0:n-1)
    pow = one(T)
    for k = 0:m-1
        a[k] = pow / Γ(β+k*α)
        pow = -pow
    end
    b[0] = zero(T)
    pow = one(T)
    for k = 1:n-1
        pow = -pow
        b[k] = pow * sinpi(k*α-β) * Γ(1+k*α-β) / π
    end
    A = zeros(T, 2r, 2r+1)
    if m ≥ r+1
        row = 0
        for k = 0:r-1
            row += 1
            A[row,k+1] = one(T)
            for j = 0:k
                A[row,r+1+j] = -a[k-j]
            end
        end
        for k = r:m-1
            row += 1
            for j = 0:r
                A[row,r+1+j] = -a[k-j]
            end
        end
        for k = r-n+1:r-1
            row += 1
            A[row,k+1] = one(T)
            for j = k+1:r
                A[row,r+1+j] = -b[j-k]
            end
        end
    else
        row = 0
        for k = 0:m-1
            row += 1
            A[row,k+1] = one(T)
            for j = 0:k
                A[row,r+1+j] = -a[k-j]
            end
        end
        for k = -(r-m):-1
            row += 1
            for j = 0:r
                A[row,r+1+j] = -b[j-k]
            end
        end
        for k = 0:r-1
            row += 1
            A[row,k+1] = one(T)
            for j = k+1:r
                A[row,r+1+j] = -b[j-k]
            end
        end
    end
    return A
end

function pade_svd(α::T, β::T, m::Integer, n::Integer) where T <: AbstractFloat
    A = pade_matrix(α, β, m, n)
    p, q = pade_svd(A)
    return p, q
end

function pade_svd(A::Matrix{T}) where T <: AbstractFloat
    r = div(size(A,1), 2)
    F = svd(A, full=true)
    p = F.V[1:r,2r+1]
    q = F.V[r+1:2r+1,2r+1]
    q1 = q[1]
    p ./= q1
    q ./= q1
    return p, q
end

function pade_lu(A::Matrix{T}) where T <: AbstractFloat
    r = div(size(A,1), 2)
    F = lu(A)
    y = zeros(T, 2r+1)
    y[2r+1] = one(T)
    for i = 2r:-1:1
        s = zero(T)
        for j = i+1:2r+1
            s += F.U[i,j] * y[j]
        end
        y[i] = -s / F.U[i,i]
    end
    p = y[1:r]
    q = y[r+1:2r+1]
    μ = q[1]
    for k = 1:r
        p[k] /= μ
    end
    for k = 1:r+1
        q[k] /= μ
    end
    return p, q, y
end

#T = BigFloat
T = Float64
α = 1 / parse(T, "2")
β = parse(T, "1")

x = range(zero(T), parse(T, "10"), length=201)
y = erfcx.(x)
figure(1)
r = 3
for m = r:5
    n = 2r + 1 - m
    A = pade_matrix(α, β, m, n)
    p, q = pade_lu(A)
    P = Polynomial(p)
    Q = Polynomial(q)
    err = P.(x) ./ Q.(x) - y
    plot(x, err, label="m, n = $m, $n")
end
legend()
grid(true)

figure(2)
for (r, m, n) in [(4, 5, 4), (8, 9, 8), (12, 13, 12)]
    A = pade_matrix(α, β, m, n)
    p, q = pade_lu(A)
    P = Polynomial(p)
    Q = Polynomial(q)
    err = P.(x) ./ Q.(x) - y
    semilogy(x, abs.(err), label="m, n = $m, $n")
end
legend()
grid(true)

figure(3)
p, q = pade_svd(1.85, 1.0, 12, 11)
P = Polynomial(p)
Q = Polynomial(q)
R(x) = P(x) / Q(x)
plot(x, R.(x))
grid(true)


text:
using JSON

#Twitter Database Type
abstract TwitterDB

#Twitter Stream Object
include("TwStream.jl")

function stream_tweets(auth::TwAuth, db::TwitterDB, api_options::Dict;
                      max_tweets::Int = 10,
                      max_time::Int = 60,
                      alert_every::Int = 100,
                      max_err_attemps::Int = 3)

  twstrm = TwStream(auth, api_options)
  streamdb(twstrm, db, max_tweets = max_tweets,
                       max_time = max_time,
                       alert_every = alert_every,
                       max_err_attemps = max_err_attemps)
end

#Use default auth
function stream_tweets(db::TwitterDB, api_options::Dict;
                      max_tweets::Int = 10,
                      max_time::Int = 60,
                      alert_every::Int = 100,
                      max_err_attemps::Int = 3)

  twstrm = TwStream(JuliaTwitter.auth, api_options)
  streamdb(twstrm, db, max_tweets = max_tweets,
                       max_time = max_time,
                       alert_every = alert_every,
                       max_err_attemps = max_err_attemps)
end

#Move tweets from the stream into the database
function streamdb(tw_stream::TwStream, db::TwitterDB;
                  max_tweets::Int = 10,
                  max_time::Int = 60,
                  alert_every::Int = 100,
                  max_err_attemps::Int = 3)

  #Create a stream buffer
  buf = start(eachline(tw_stream.twstream))
  tweet_count = 0
  err_attempts = 0

  #Calculate the ending time (convert  mins to secs to ms)
  end_time = time() +  60 * max_time

  while tweet_count < max_tweets
    #Exit the loop if we have passed the time limit
    if time() > end_time
      println("Time limit reached. Collected $tweet_count tweets")
      break
    end

    #Capture the next line and save the buf state
    line,buf = next(eachline(tw_stream.twstream), buf)

    try
      #Parse JSON
      dict = JSON.parse(line)
      add_entry(db, dict)
      tweet_count+=1

      #reset err_attempts counter
      err_attempts = 0
    catch stream_error
      #Twitter will return an empty line if there hasn't been a tweet in awhile
      #     that should not count as an error

      #If the line is not empty, alert the user
      if length(line) > 1
        err_attempts+=1

        println("[WARNING] Error processing entry: $line")
        println("          Continuing execution...")

        #Check err attempts
        if err_attempts >= max_err_attemps
          println("[ERROR] Maximum number of consecutive errors reached. Exiting...")

          #Try closing the database
          try
            close_db(db)
          catch
          end

          return false
        end

      end
    end

    #Update user
    if tweet_count % alert_every == 0 && tweet_count != 0
      println(STDOUT, "Wrote tweet $tweet_count...")
    end
  end

  close_db(db)

  true
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
# Sampler
# In: parameter space
# Out: parameter sets
#
#
# For step in budget:
#
#   For rep in reps:
#
#       Model
#       In: parameter set[step]
#       Out: model data
#
#       Evaluation
#       In: model data & either true data or stylised facts
#       Out: labelled data (for one run)
#
#   Result averaging
#   In: Labels for each run
#   Out: Aggregated labels (one for each criterion)
#
#   If step is big enough to run surrogate with:
#
#       Surrogate
#       In: labelled data
#       Out: predictions for unlabelled data
#
#       Selector
#       In: predictions for unlabelled data
#       Out: unlabelled parameter set


text:
export findcallers

function get_typed_instances!(srcs, @nospecialize(tt), method::Method, world, interp)
    # This is essentially taken from code_typed_by_type
    matches = Base._methods_by_ftype(tt, -1, world)
    if matches === false
        error("signature $(item.specTypes) does not correspond to a generic function")
    end
    for match in matches
        match.method == method || continue
        meth = Base.func_for_method_checked(match.method, tt, match.sparams)
        (src, ty) = isdefined(Core.Compiler, :NativeInterpreter) ?
            Core.Compiler.typeinf_code(interp, meth, match.spec_types, match.sparams, false) :
            Core.Compiler.typeinf_code(meth, match.spec_types, match.sparams, false, interp)
        push!(srcs, (src, match.sparams))
    end
    return srcs
end
get_typed_instances!(srcs, mi::MethodInstance, world, interp) = get_typed_instances!(srcs, mi.specTypes, mi.def, world, interp)

defaultinterp(world) = isdefined(Core.Compiler, :NativeInterpreter) ?
                       Core.Compiler.NativeInterpreter(world) :
                       Core.Compiler.Params(world)

function get_typed_instances(mi::MethodInstance; world=typemax(UInt), interp=defaultinterp(world))
    return get_typed_instances!(Tuple{CodeInfo,Core.SimpleVector}[], mi, world, interp)
end
function get_typed_instances(@nospecialize(tt), method::Method; world=typemax(UInt), interp=defaultinterp(world))
    return get_typed_instances!(Tuple{CodeInfo,Core.SimpleVector}[], tt, method, world, interp)
end

"""
    CallMatch

A structure to summarize a "matching" caller/callee pair. The fields are:

- `mi`: the `MethodInstance` for the caller
- `src`: its corresponding `CodeInfo`
- `sparams`: the type parameters for the caller, given `mi`'s signature (alternatively use `mi.sparam_vals`)
- `line`: the statement number (in SSAValue sense) on which the call occurs
- `argtypes`: the caller's inferred types passed as arguments to the callee
"""
struct CallMatch
    mi::MethodInstance
    src::CodeInfo
    sparams::SimpleVector
    line::Int
    argtypes::Vector{Any}
end

function Base.show(io::IO, cm::CallMatch)
    print(io, "CallMatch from ", cm.mi, " on statement ", cm.line, "/", length(cm.src.code), " with inferred argument types ", cm.argtypes)
end

"""
    callers = findcallers(f, argmatch::Union{Nothing,Function}, mis; callhead=:call | :iterate)

Find "static" callers of a function `f` with *inferred* argument types for which `argmatch(types)` returns true.
Optionally pass `nothing` for `argmatch` to allow any calls to `f`.
`mis` is the list of `MethodInstance`s you want to check, for example obtained from [`methodinstances`](@ref).

`callhead` controls whether you're looking for an ordinary (`:call`) or a splatted (varargs) call (`:iterate`).

`callers` is a vector of [`CallMatch`](@ref) objects.

# Examples

```julia
f(x) = rand()
function g()
    a = f(1.0)
    z = Any[7]
    b = f(z[1]::Integer)
    c = f(z...)
    return nothing
end

g()
mis = union(methodinstances(f), methodinstances(g))

# callhead = :call is the default
callers1 = findcallers(f, argtyps->length(argtyps) == 1 && argtyps[1] === Float64, mis)

# Get the partial-inference case with known number of arguments (this is definitely :call)
callers2 = findcallers(f, argtyps->length(argtyps) == 1 && argtyps[1] === Integer, mis; callhead=:call)

# Get the splat call
callers3 = findcallers(f, argtyps->length(argtyps) == 1 && argtyps[1] === Vector{Any}, mis; callhead=:iterate)
```

!!! compat
    `findcallers` is available on Julia 1.6 and higher

!!! warn
    `findcallers` is not guaranteed to find all calls. Calls can be "obfuscated" by many mechanisms,
    including calls from top level, calls where the function is a runtime variable, etc.
"""
function findcallers(f, argmatch::Union{Function,Nothing}, mis::AbstractVector{Core.MethodInstance};
                     callhead::Symbol=:call, world=typemax(UInt), interp=defaultinterp(world))
    callhead === :call || callhead === :invoke || callhead === :iterate || error(":call and :invoke are supported, got ", callhead)
    # Check that f is not a type with specified parameters
    if f isa DataType && !isempty(f.parameters)
        throw(ArgumentError("findcallers does not support constructor calls with type parameters, got $f. Pass in the type without parameters."))
    end
    # Construct a GlobalRef version of `f`
    ref = GlobalRef(parentmodule(f), nameof(f))
    callers = CallMatch[]
    srcs = Tuple{CodeInfo,Core.SimpleVector}[]
    for item in mis
        empty!(srcs)
        try
            get_typed_instances!(srcs, item, world, interp)
        catch err
            # println("skipping ", item, ": ", err)
            continue
        end
        for (src, sparams) in srcs
            for (i, stmt) in enumerate(src.code)
                if isa(stmt, Expr)
                    if stmt.head === :(=)
                        stmt = stmt.args[2]   # call must come from right hand side of assignment
                        isa(stmt, Expr) || continue
                    end
                    stmt = stmt::Expr
                    callee = nothing
                    if stmt.head === callhead
                        callee = stmt.args[1]
                    elseif callhead === :iterate && stmt.head === :call && isglobalref(stmt.args[1], Core, :_apply_iterate)
                        callee = stmt.args[3]
                    end
                    callee === nothing && continue
                    if isa(callee, Core.SSAValue)
                        callee = try
                            eval_ssa(src, sparams, callee.id)
                        catch err
                            @show src i
                            throw(err)
                        end
                    end
                    matches = false
                    if callee === f
                        matches = true
                    elseif isa(callee, GlobalRef)
                        if isglobalref(callee, ref.mod, ref.name) ||
                           (isdefined(callee.mod, callee.name) && getfield(callee.mod, callee.name) === f)
                            matches = true
                        end
                    elseif isa(f, Type) && isa(callee, DataType) && !isempty(callee.parameters)
                        # This came from a `typeof(a)(b)` call
                        matches = callee <: f
                    end
                    matches || continue
                    # Collect the arg types
                    argtypes = []
                    for i = (callhead === :iterate ? 4 : 2):length(stmt.args)
                        a = stmt.args[i]
                        if isa(a, Core.SSAValue)
                            push!(argtypes, extract(src.ssavaluetypes[a.id], sparams))
                        elseif isa(a, Core.SlotNumber)
                            push!(argtypes, extract(src.slottypes[a.id], sparams))
                        elseif isa(a, GlobalRef)
                            push!(argtypes, Core.Typeof(getfield(a.mod, a.name)))
                        elseif isexpr(a, :static_parameter)
                            a = a::Expr
                            push!(argtypes, Type{sparams[a.args[1]::Int]})
                        else
                            push!(argtypes, extract(a, sparams))
                        end
                    end
                    # Check whether the types match
                    if argmatch === nothing
                        push!(callers, CallMatch(item, src, sparams, i, argtypes))
                        break
                    else
                        if argmatch(argtypes)
                            push!(callers, CallMatch(item, src, sparams, i, argtypes))
                            break
                        end
                    end
                end
            end
        end
    end
    return callers
end

isglobalref(@nospecialize(g), mod::Module, name::Symbol) = isa(g, GlobalRef) && g.mod === mod && g.name === name

extract(a, sparams) = isa(a, Core.Const) ? Core.Typeof(a.val) :
                      isa(a, Core.PartialStruct) ? (a.typ <: Tuple{Any} ? a.typ.parameters[1] : a.typ) :
                      isa(a,Type) ? a :
                      isexpr(a, :static_parameter) ? sparams[(a::Expr).args[1]] :
                      Core.Typeof(a)

# This is deliberately simple to prevent performance from tanking
function eval_ssa(src, sparams, id)
    stmt = src.code[id]
    if isa(stmt, Expr)
        if stmt.head === :call
            callee = stmt.args[1]
            if isglobalref(callee, Core, :apply_type)
                return stmt.args[2]
            elseif isglobalref(callee, Base, :getproperty)
                modg, objq = stmt.args[2], stmt.args[3]
                if isa(modg, Core.SSAValue)
                    modg = eval_ssa(src, sparams, modg.id)
                end
                if isa(modg, GlobalRef)
                    isdefined(modg.mod, modg.name) || return nothing
                    mod = getfield(modg.mod, modg.name)
                    objname = objq.value::Symbol
                    return GlobalRef(mod, objname)
                end
            elseif isglobalref(callee, Core, :kwfunc)
                return stmt.args[2]
            elseif isa(callee, GlobalRef) && callee.name === :typeof && length(stmt.args) == 2
                isa(stmt.args[2], Core.SlotNumber) && return src.slottypes[(stmt.args[2]::Core.SlotNumber).id]
                isa(stmt.args[2], Core.SSAValue) && return eval_ssa(src, sparams, (stmt.args[2]::Core.SSAValue).id)
            end
        end
    end
    isa(stmt, GlobalRef) && return stmt
    isa(stmt, QuoteNode) && return stmt.value
    return src.ssavaluetypes[id]
end


text:
using NonstationaryProcesses
using SafeTestsets
using Test

@testset "NonstationaryProcesses.jl" begin
    @test 1 == 1
end


text:
module ICARTTUnits

using Unitful;

using ..ICARTTExceptions: ICARTTUnitsException;

# The unit aliases dictionary defines aliases for existing Unitful units that
# might be written differently in ICARTT files. The key should be the Unitful
# abbreviation string, and the value must be a tuple of alternate aliases that
# that unit may go by in the ICARTT files. It must be a tuple, even if there is
# only one alias. Note that the aliases will be matched in order of decreasing
# length, no matter what order they are listed in.
const unit_aliases = Dict{String, AbstractArray{AbstractString,1}}();

# Treat "unitless" and "none" slightly differently: "unitless" implies that it
# is a physical quantity that has no dimensionality; while "none" implies that
# it is not a physical quantity, e.g. an index of some sort.
@unit(unitless, "unitless", unitless, Quantity(1, Unitful.NoUnits), false);
@unit(none, "no_units", none, Quantity(1, Unitful.NoUnits), false);
# mach numbers are technically unitless b/c the represent the speed of the craft
# relative to the speed of sound _in that medium_ so they do not map to an
# absolute speed
@unit(mach, "mach", MachNumber, Quantity(1, Unitful.NoUnits), false);

# units relating to amount of matter
@unit(molec, "molec.", molec, Quantity(1/Unitful.Na.val, Unitful.mol), false); # can't use // because Na is a float
@unit(DU, "DU", DobsonUnit, Quantity(0.4462, Unitful.mmol), false);

# the UnitfulMoles package was not working on 5 Dec 2018, so I'm defining
# mixing ratios here
@unit(ppm, "ppm", parts_per_million, Quantity(1, Unitful.μmol / Unitful.mol), false);
@unit(ppmv, "ppmv", parts_per_million_volume, Quantity(1, Unitful.μL / Unitful.L), false);
@unit(ppb, "ppb", parts_per_billion, Quantity(1, Unitful.nmol / Unitful.mol), false);
@unit(ppbv, "ppb", parts_per_billion_volume, Quantity(1, Unitful.nL / Unitful.L), false);
@unit(ppt, "ppt", parts_per_trillion, Quantity(1, Unitful.pmol / Unitful.mol), false);
@unit(pptv, "pptv", parts_per_trillion_volume, Quantity(1, Unitful.pL / Unitful.L), false);

# time units
@unit(days, "days", days, Quantity(24, Unitful.hr), false);

# volume units
@unit(std_m, "std m", m_at_stp, Quantity(1, Unitful.m), false);

# Recommended by http://ajkeller34.github.io/Unitful.jl/stable/extending/
# if a package gets precompiled

const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(ICARTTUnits)

    # this part will initialize the units aliases dict
    merge!(unit_aliases, _read_alias_config(joinpath(@__DIR__, "common_data", "standard_unit_aliases.txt")));
end

function list_default_unit_aliases()
    for (key,val) in unit_aliases
        println("\"$key\" will be substituted for: \"$(join(val, "\", \""))\"");
    end
end

"""
    parse_unit_string(ustr; aliases_dict=nothing)

Given a string describing a unit or combination of units, convert it into a
Unitful.Units instance. This uses `sanitize_raw_unit_strings` to preformat the
string into a format that Unitful is more likely to understand.

`aliases_dict` is passed through to `sanitize_raw_unit_strings`.
"""
function parse_unit_string(ustr; aliases_dict=nothing)
    # Since I could not find a version of the Unitful @u_str macro that was a
    # function, I'm replicating the internals of @u_str here
    if aliases_dict == nothing
        ustr = sanitize_raw_unit_strings(ustr);
    else
        ustr = sanitize_raw_unit_strings(ustr, aliases_dict)
    end
    # The way the @u_str macro works in Unitful is that, given an expression or
    # symbol, tries to replace each symbol (whether standalone or in the expression)
    # with a Units instance. Then the final expression is automatically evaluated
    # by the macro before returning, so e.g. "m * s^-1" becomes:
    #   FreeUnit(m) / FreeUnit(s)
    # which gets evaluated to
    #   FreeUnit(m / s)
    # When reading from a file, I don't think there's any way to replicate that
    # without doing a run-time eval. This should be safe, as Unitful.replace_value
    # will error if a symbol is not recognized, and replaces all symbols with the
    # apporpriate Unitful.FreeUnit instances.
    try
        expr = Unitful.replace_value(Meta.parse(ustr));
        return eval(expr); # do not need the `esc` call since not returning from macro
    catch err
        # Since ICARTT files are ASCII encoded, the "micro" prefix will be represent
        # most often by a "u", but we can't just replace "u" with "μ" in all cases,
        # that would break e.g. "unitless". The only way I can figure out that doesn't
        # involve manually entering _every_ unit with the prefix "μ" that we want
        # to handle is to try the parsing and, if it fails because a symbol starting
        # with "u" isn't recognized, replace that with "μ" and try again.
        if ~(:msg in fieldnames(typeof(err)))
            # if no message, then not the right error
            rethrow(err);
        end
        m = match(r"(?<=Symbol )u[a-zA-Z]+", err.msg);
        if m === nothing
            rethrow(err);
        else
            orig_str = ustr
            new_str = replace(ustr, Regex("$(m.match)") => SubstitutionString("μ" * m.match[2:end]), count=1);
            try
                return parse_unit_string(new_str);
            catch err2
                msg = "Tried replacing 'u' prefix with 'μ' ($orig_str -> $new_str) but this failed: $(err2)";
                throw(ICARTTUnitsException(msg));
            end
        end
    end
end

"""
    sanitize_raw_unit_strings(ustr)
    sanitize_raw_unit_strings(ustr, aliases_dict=nothing)

Preprocesses a unit string (`ustr`) read in from ICARTT files into a form that
can be understood by Unitful. Does several things:

    1. Replaces any substrings defined by `unit_aliases` with their key value.
       This helps standardize the units used.
    2. Replaces blank space between an alphanumeric character and a letter with
       a "*" e.g. "m s-1" becomes "m * s-1"
    3. Inserts a "^" between units and their exponents; specifically, if a letter
       is followed immediately by a number, +, or -, a "^" is inserted.

The second, two argument form, allows you to specify the dictionary of aliases
to use.
"""
function sanitize_raw_unit_strings(ustr::AbstractString)
    return sanitize_raw_unit_strings(ustr, unit_aliases);
end

function sanitize_raw_unit_strings(ustr::AbstractString, aliases_dict::AbstractDict)
    #print("Sanitizing '$ustr', ")
    # The first one is the most complicated because we need to look for any of
    # the aliases defined in the unit_aliases dictionary, but they need to match
    # a whole word, or be prefixed. That makes this a bit of a mess. We require
    # either:
    #   A) the string to substitute (S) is preceeded by a start of the string, or
    #   B) S is preceeded by a whitspace character, or
    #   C) S is preceeded by one of the defined metric prefixes, which is itself
    #       preceeded by a start-of-string or whitespace
    # and
    #   D) S is followed by a non-letter character, or an end-of-string

    # This is the first step, we need to construct a look-ahead pattern that matches
    # on start-of-string (\A), whitespace (\s), or one of the prefixes preceeded
    # by a start-of-string or whitespace.
    #
    # The joins will correctly put the | and either \A or \s in front of each
    # prefix except the first one, so the first "\\A" adds it for the first prefix
    # and the "|\\s" adds it for the first prefix of the second group. We shouldn't
    # need to explicitly include \A and \s without a prefix because the prefixdict
    # includes an empty string for no prefix.
    #
    # This part relies heavily on ICARTT files being ASCII encoded; it's not compatible
    # with Unicode encodings. This has to happen first in case we want to treat
    # a string with spaces, e.g. "std m" specially.
    prefixes = "\\A" * join(values(Unitful.prefixdict), "|\\A") * "|\\s" * join(values(Unitful.prefixdict), "|\\s")
    for (key, val) in pairs(aliases_dict)
        # Force the aliases to be searched in order of decreasing length (longest
        # first); this ensures that if a shorter alias is a subset of a longer one
        # that the entirety of the longer one gets matched. E.g. if searching
        # for "deg" or "degrees" and "deg" came up first, it would replace only
        # the first three letters of "degrees" with the proper ° symbol, i.e. we'd
        # get "°rees", which wouldn't match any known unit.\
        # Also there's no `sort` method for tuples, so we have to convert to a
        # temporary array
        sorted_aliases = sort([val...], by=length, rev=true)
        aliases = join(val, "|");
        sub = SubstitutionString(key)
        # Look for any of the aliases preceeded by any of the metric prefixes and
        # either the start-of-string or whitespace, and succeeded by any non-letter
        # character or an end-of-string
        re = Regex("(?<=$(prefixes))($(aliases))(?=[^a-zA-Z]|\\Z)");
        #println("sub = $sub, re = $re")
        ustr = replace(ustr, re => sub)
    end

    # Replace any spaces between a letter/number and letter with *, so e.g.
    # m s-1 -> m * s-1 and m2 s-1 -> m2 * s-1. User groups to keep the last
    # character of the preceeding unit and first letter of the succeeding unit
    # in the result
    ustr = replace(ustr, r"([a-zA-Z0-9])\s+([a-zA-Z])" => s"\1 * \2")

    # Insert a caret between units and their exponents. Assumes units will never
    # contain digits. Gets the last non-digit character of the unit as group 1
    # and the exponent which is the numbers and optionally a leading + or - as
    # group 2, and insert the ^ between them. Right now, there cannot be spaces
    # between the unit and exponent.
    ustr = replace(ustr, r"([a-zA-Z])([+\-]?\d+)" => s"\1^\2")

    #println(" result: '$ustr'")
    return ustr
end

"""
    _read_alias_config(config_file; verbose=0)

Read a unit alias configuration file `config_file`. Returns the aliases as a
dictionary where the key is the string to replace with and the values are arrays
of strings to replace.

`verbose` sets the level of logging, change to > 0 to increase logging or < 0 to
suppress warning messages.

See `ReadICARTT.read_icartt_file` for information about the formatting of the
configuration file.
"""
function _read_alias_config(config_file; verbose=0)
    alias_dict = Dict{String,Array}();
    return _read_alias_config!(alias_dict, config_file; verbose=verbose);
end

"""
    _read_alias_config!(alias_dict, config_file; verbose=0)

Reads the unit alias configuration file, placing the results into `alias_dict`.
The modified `alias_dict` is also returned. Values in `alias_dict` may be added
to or overwritten depending on the configuration file.

`verbose` sets the level of logging, change to > 0 to increase logging or < 0 to
suppress warning messages.

See `ReadICARTT.read_icartt_file` for information about the formatting of the
configuration file.
"""
function _read_alias_config!(alias_dict, config_file; verbose=0)
    # the config file needs to be formatted as
    #  Unitful abbrev: alias1, alias2, ... [:{append,replace}]
    open(config_file, "r") do io
        for (ln,line) in enumerate(eachline(io))
            line_chunks = split(line, ":");
            unit = strip(line_chunks[1]);
            aliases = [strip(a) for a in split(line_chunks[2], ",")];
            op_mode = length(line_chunks) > 2 ? lowercase(strip(line_chunks[3])) : "append";

            if op_mode == "replace"
                if unit in keys(alias_dict) && verbose >= 0
                    old_aliases = join(alias_dict[unit], ", ");
                    new_aliases = join(aliases, ", ");
                    @warn "Replacing existing aliases for unit \"$unit\" ($old_aliases) with
                     those defined on line $ln of $config_file ($new_aliases)"
                end
                alias_dict[unit] = aliases;
            elseif op_mode == "append"
                if !(unit in keys(alias_dict))
                    alias_dict[unit] = aliases;
                else
                    append!(alias_dict[unit], aliases);
                end
            end
        end
    end
    return alias_dict;
end

end

max_stars_repo_name: joshua-laughner/Julia-ICARTT

max_issues_repo_name: joshua-laughner/Julia-ICARTT

max_forks_repo_name: joshua-laughner/Julia-ICARTT


text:
struct FieldDataset{F, M, P}
    fields :: F
  metadata :: M
  filepath :: P
end

"""
    FieldDataset(filepath;
                 architecture=CPU(), grid=nothing, backend=InMemory(), metadata_paths=["metadata"])

Returns a `Dict` containing a `FieldTimeSeries` for each field in the JLD2 file located
at `filepath`. Note that model output **must** have been saved with halos.

Keyword arguments
=================
- `backend`: Either `InMemory()` (default) or `OnDisk()`. The `InMemory` backend will
load the data fully in memory as a 4D multi-dimensional array while the `OnDisk()`
backend will lazily load field time snapshots when the `FieldTimeSeries` is indexed
linearly.

- `metadata_paths`: A list of JLD2 paths to look for metadata. By default it looks in
  `file["metadata"]`.

- `grid`: May be specified to override the grid used in the JLD2 file.
"""
function FieldDataset(filepath;
                    architecture=CPU(), grid=nothing, backend=InMemory(), metadata_paths=["metadata"])

  file = jldopen(filepath)

  field_names = keys(file["timeseries"])
  filter!(k -> k != "t", field_names)  # Time is not a field.

  ds = Dict{String, FieldTimeSeries}(
      name => FieldTimeSeries(filepath, name; architecture, backend, grid)
      for name in field_names
  )

  metadata = Dict(
      k => file["$mp/$k"]
      for mp in metadata_paths if haskey(file, mp)
      for k in keys(file["$mp"])
  )

  close(file)

  return FieldDataset(ds, metadata, abspath(filepath))
end

Base.getindex(fds::FieldDataset, inds...) = Base.getindex(fds.fields, inds...)

Base.show(io::IO, fds::FieldDataset) where {X, Y, Z, K, A} =
  print(io, "FieldDataset with $(length(fds.fields)) fields and $(length(fds.metadata)) metadata entries.")


text:
module PokerHandEvaluator

using LinearAlgebra: dot
using PlayingCards
using Combinatorics

export CompactHandEval, FullHandEval
export hand_rank, hand_type, best_cards, all_cards

abstract type AbstractHandEvaluation end

"""
    hand_rank(::AbstractHandEvaluation)

The hand rank, between 1 (straight royal flush)
and 7462 (7 high with lowest kickers):

    (A♡,K♡,Q♡,J♡,10♡) -> 1
    ...
    (7♡,5♢,4♣,3♠,2♡) -> 7462
"""
function hand_rank end

"""
    hand_type(::AbstractHandEvaluation)

A Symbol corresponding to the type of hand.
Concretely, `hand_type` returns one of:
 - `:straight_flush`
 - `:quads`
 - `:full_house`
 - `:flush`
 - `:straight`
 - `:trips`
 - `:two_pair`
 - `:one_pair`
 - `:high_card`
"""
function hand_type end

"""
    best_cards(::FullHandEval)

The set of cards with the lowest (best) rank.
"""
function best_cards end

"""
    all_cards(::FullHandEval)

All of the input cards to `FullHandEval`.
"""
function all_cards end

"""
    FullHandEval(cards::Card...)
    FullHandEval(::Card,::Card,::Card,::Card,::Card[,::Card,::Card])
    FullHandEval(::Tuple{Card,Card,Card,Card,Card[,::Card,::Card]})

A compact hand evaluation, with supporting methods:
[`hand_rank`](@ref), [`hand_type`](@ref), [`best_cards`](@ref), [`all_cards`](@ref).

# Examples

```julia
using PlayingCards, PokerHandEvaluator
table_cards = (J♡,J♣,2♣,3♢,5♣)
player_cards = (
  (A♠,2♠,table_cards...),
  (J♠,T♣,table_cards...),
);
fhe = FullHandEval.(player_cards);
winner_id = argmin(hand_rank.(fhe)) # = 2
winning_hand = hand_type(fhe[winner_id]) # = :trips
winning_rank = hand_rank(fhe[winner_id]) # = 1842
winning_cards = best_cards(fhe[winner_id]) # = (J♠, T♣, J♡, J♣, 5♣)
allcards = all_cards(fhe[winner_id]) # = (J♠, T♣, J♡, J♣, 2♣, 3♢, 5♣)
```
"""
struct FullHandEval{AC <: Tuple, TC <: Tuple} <: AbstractHandEvaluation
    hand_type::Symbol
    rank::Int
    all_cards::AC
    best_cards::TC
    FullHandEval(player_cards::Tuple, table_cards::Tuple) =
        FullHandEval((player_cards..., table_cards...))
    FullHandEval(all_cards...) = FullHandEval(all_cards)
    function FullHandEval(all_cards::Tuple)
        @assert 5 ≤ length(all_cards) ≤ 7
        rank, hand_type, best_cards = evaluate_full(all_cards)

        return new{
            typeof(all_cards),
            typeof(best_cards)}(hand_type, rank, all_cards, best_cards)
    end
end

hand_rank(th::FullHandEval) = th.rank
hand_type(th::FullHandEval) = th.hand_type
best_cards(th::FullHandEval) = th.best_cards
all_cards(th::FullHandEval) = th.all_cards

"""
    CompactHandEval(cards::Card...)
    CompactHandEval(::Card,::Card,::Card,::Card,::Card[,::Card,::Card])
    CompactHandEval(::Tuple{Card,Card,Card,Card,Card[,::Card,::Card]})

A compact hand evaluation, with supporting methods:
[`hand_rank`](@ref), [`hand_type`](@ref).

# Examples

```julia
using PlayingCards
using PokerHandEvaluator
table_cards = (J♡,J♣,A♣,A♢);
player_cards = (
  (A♠,2♠,table_cards...),
  (J♠,J♣,table_cards...),
);
che = CompactHandEval.(player_cards);
winner_id = argmin(hand_rank.(che)) # = 2
winning_hand = hand_type(che[winner_id]) # = :quads
```
"""
struct CompactHandEval <: AbstractHandEvaluation
    rank::Int
    CompactHandEval(player_cards::Tuple, table_cards::Tuple) =
        CompactHandEval((player_cards..., table_cards...))
    CompactHandEval(all_cards...) = CompactHandEval(all_cards)
    function CompactHandEval(all_cards::Tuple)
        @assert 5 ≤ length(all_cards) ≤ 7
        return new(evaluate_compact(all_cards))
    end
end

hand_type(th::CompactHandEval) =
    hand_type_binary_search(hand_rank(th))
hand_rank(th::CompactHandEval) = th.rank

"""
    evaluate(cards::Card...)
    evaluate(::Card,::Card,::Card,::Card,::Card[,::Card,::Card])
    evaluate(::Tuple{Card,Card,Card,Card,Card[,Card,Card]})

Evaluates 5, 6, and 7-card hands using
[`evaluate5`](@ref). This is done by using
`Combinatorics.combinations` to evaluate
all possible 5, 6, and 7-card hand combinations.
"""
function evaluate end

function hand_type_binary_search(rank::Int)
    @assert 1 ≤ rank ≤ 7462
    if rank ≤ 1609
        if rank ≤ 166
            rank ≤ 10 && return :straight_flush
            return :quads
        else
            rank ≤ 322 && return :full_house
            rank ≤ 1599 && return :flush
            return :straight
        end
    else
        if rank ≤ 3325
            rank ≤ 2467 && return :trips
            return :two_pair
        else
            rank ≤ 6185 && return :one_pair
            return :high_card
        end
    end
end

#####
##### evaluate_full
#####

# Wrapper for 5 card hands:
function evaluate_full(t::Tuple{<:Card,<:Card,<:Card,<:Card,<:Card})
    hand_rank = evaluate5(t)
    return (hand_rank, hand_type_binary_search(hand_rank), t)
end

# Wrapper for 6 card hands:
function evaluate_full(t::Tuple{<:Card,<:Card,<:Card,<:Card,<:Card,<:Card})
    return best_hand_rank_from_6_cards_full(t)
end
# Wrapper for 7 card hands:
function evaluate_full(t::Tuple{<:Card,<:Card,<:Card,<:Card,<:Card,<:Card,<:Card})
    return best_hand_rank_from_7_cards_full(t)
end

#####
##### evaluate_compact
#####

# Wrapper for 5 card hands:
evaluate_compact(t::Tuple{<:Card,<:Card,<:Card,<:Card,<:Card}) = evaluate5(t)

# Wrapper for 6 card hands:
function evaluate_compact(t::Tuple{<:Card,<:Card,<:Card,<:Card,<:Card,<:Card})
    return best_hand_rank_from_6_cards_compact(t)
end
# Wrapper for 7 card hands:
function evaluate_compact(t::Tuple{<:Card,<:Card,<:Card,<:Card,<:Card,<:Card,<:Card})
    return best_hand_rank_from_7_cards_compact(t)
end

include("best_hand_rank_from_n_cards.jl")
include("HandCombinations.jl")
include("evaluate5.jl")

end # module


text:
function gensys_uncertain_altpol(m::AbstractDSGEModel, prob_vec::AbstractVector{S},
                                 altpolicies::Vector{AltPolicy} = [default_policy()];
                                 regime_switching::Bool = false, regimes::Vector{Int} = Int[1],
                                 TTT::Matrix{S} = Matrix{S}(undef, 0, 0),
                                 Γ0s::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0),
                                 Γ1s::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0),
                                 Cs::Vector{Vector{S}} = Vector{Vector{S}}(undef, 0),
                                 Ψs::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0),
                                 Πs::Vector{Matrix{S}} = Vector{Matrix{S}}(undef, 0)) where {S <: Real}

    @assert sum(prob_vec) == 1. "The vector of probabilities must sum to 1"

    @assert length(altpolicies) == length(prob_vec) - 1

    altpolicies = altpolicies[prob_vec[2:end] .> 0.0]
    prob_vec = vcat(prob_vec[1], prob_vec[2:end][prob_vec[2:end] .> 0.0])
    if regime_switching
        if length(regimes) == 1
            Γ0, Γ1, C, Ψ, Π = haskey(get_settings(m), :regime_eqcond_info) &&
                haskey(get_setting(m, :regime_eqcond_info), regimes[1]) ?
                get_setting(m, :regime_eqcond_info)[regimes[1]].alternative_policy.eqcond(m, regimes[1]) :
                eqcond(m, regimes[1])

            if isempty(TTT)
                assert_cond = haskey(get_settings(m), :uncertain_altpolicy) ? !get_setting(m, :uncertain_altpolicy) : true
                @assert assert_cond "The Setting :uncertain_altpolicy must be false if TTT is empty."
                TTT, _, _ = solve(m; regime_switching = regime_switching, regimes = regimes)
            end

            Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til = gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π)

            inds = 1:n_states(m) # Don't want augmented states
            Td = Vector{Matrix{S}}(undef, length(prob_vec))
            Cd = Vector{Vector{S}}(undef, length(prob_vec))
            Td[1] = TTT[inds, inds]
            Cd[1] = C[inds]

            # Turn off uncertain altpolicy b/c when computing the state space system under alternative policies,
            # we want to treat these alternative policies as perfectly credible.
            is_uncertain_altpol = haskey(get_settings(m), :uncertain_altpolicy) && get_setting(m, :uncertain_altpolicy)
            if is_uncertain_altpol
                m <= Setting(:uncertain_altpolicy, false)
            end

            # Solve each alternative policy
            for (i, altpolicy) in enumerate(altpolicies)
                j = i + 1
                tmpT, _, tmpC = if altpolicy.solve == solve      # if solve is being called for gensys_uncertain_altpol
                    altpolicy.solve(m; regime_switching = false) # then it is interpreted as the "default" policy, and
                else                                             # it is further assumed that we can safely ignore regime-switching
                    altpolicy.solve(m; regime_switching = regime_switching, regimes = regimes)
                end
                Td[j] = tmpT[inds, inds]
                Cd[j] = tmpC[inds]
            end

            # Turn uncertain_altpolicy back on
            if is_uncertain_altpol
                m <= Setting(:uncertain_altpolicy, true)
            end

            T̄ = sum([p .* Td[i] for (i, p) in enumerate(prob_vec)])
            C̄ = sum([p .* Cd[i] for (i, p) in enumerate(prob_vec)])
            Lmat = (Γ2_til * T̄ + Γ0_til)
            Tcal = Lmat \ Γ1_til
            Rcal = Lmat \ Ψ_til
            Ccal = Lmat \ (C_til - Γ2_til * C̄)
            return Tcal, Rcal, Ccal
        else
            error("Regime switching for multiple regimes has not been completed for gensys_uncertain_altpol")
            if isempty(Γ0s) || isempty(Γ1s) || isempty(Cs) || isempty(Ψs) || isempty(Πs)
                # If any of these are empty, we recompute the relevant matrices
                for fcast_reg in regimes
                    Γ0s[fcast_reg], Γ1s[fcast_reg], Cs[fcast_reg], Ψs[fcast_reg], Πs[fcast_reg] =
                        haskey(get_settings(m), :regime_eqcond_info) ?
                        (haskey(get_setting(m, :regime_eqcond_info), fcast_reg) ?
                         get_setting(m, :regime_eqcond_info)[fcast_reg].alternative_policy.eqcond(m, fcast_reg) :
                         eqcond(m, fcast_reg)) : eqcond(m, fcast_reg)
                end
            end

            # Solve model for creating the "default" matrices in each period
            TTTs_gensys, RRRs_gensys, CCCs_gensys, eus =
                map((Γ0, Γ1, C, Ψ, Π) -> gensys(Γ0, Γ1, C, Ψ, Π, 1 + 1e-6, verbose = verbose), Γ0s, Γ1s, Cs, Ψs, Πs)

            # Check for LAPACK exception, existence and uniqueness
            for eu in eus
                if eu[1] != 1 || eu[2] != 1
                    throw(GensysError())
                end
            end

            TTTs_gensys = map(TTT_gensys -> real(TTT_gensys), TTTs_gensys)
            RRRs_gensys = map(RRR_gensys -> real(RRR_gensys), RRRs_gensys)
            CCCs_gensys = map(CCC_gensys -> real(CCC_gensys), CCCs_gensys)

            Γ0_tils, Γ1_tils, Γ2_tils, C_tils, Ψ_tils =
                map((Γ0, Γ1, C, Ψ, Π) -> gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π), Γ0s, Γ1s, Cs, Ψs, Πs)

            # Rows are different altpolicies in each regime, columns are different regimes
            Td = Matrix{Matrix{S}}(undef, length(prob_vec))
            Cd = Matrix{Vector{S}}(undef, length(prob_vec))
            Td[1, :] = TTTs_gensys
            Cd[1, :] = CCCs_gensys
            for fcast_reg in regimes
                for (i, altpolicy) in enumerate(altpolicies)
                    j = i + 1
                    tmpT, _, tmpC = altpolicy.solve(m; regime_switching = regime_switching, regimes = regimes)
                    Td[j, fcast_reg] = tmpT
                    Cd[j, fcast_reg] = tmpC
                end
            end

            T̄s = map(reg -> sum([p .* Td_reg[j, reg] for (j, p) in enumerate(prob_vec)]), regimes)
            C̄s = map(reg -> sum([p .* Cd_reg[j, reg] for (j, p) in enumerate(prob_vec)]), regimes)

            Tcal = Vector{Matrix{S}}(undef, length(regimes))
            Rcal = Vector{Matrix{S}}(undef, length(regimes))
            Ccal = Vector{Vector{S}}(undef, length(regimes))

            for i in length(regimes)
                Lmat = (Γ2_tils[i] * T̄s[i] + Γ0_tils[i])
                Tcals[i] = Lmat \ Γ1_tils[i]
                Rcals[i] = Lmat \ Ψ_tils[i]
                Ccals[i] = Lmat \ (C_tils[i] - Γ2_tils[i] * C̄s[i])
            end

            return Tcals, Rcals, Ccals
        end
    else
        Γ0, Γ1, C, Ψ, Π = haskey(get_settings(m), :regime_eqcond_info) ?
            (haskey(get_setting(m, :regime_eqcond_info), regimes[1]) ?
             get_setting(m, :regime_eqcond_info)[regimes[1]].alternative_policy.eqcond(m, regimes[1]) :
             eqcond(m, regimes[1])) : eqcond(m, regimes[1])

        if isempty(TTT)
            assert_cond = haskey(get_settings(m), :uncertain_altpolicy) ? !get_setting(m, :uncertain_altpolicy) : true
            @assert assert_cond "The Setting :uncertain_altpolicy must be false if TTT is empty."
            TTT, _, _ = solve(m)
        end

        Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til = gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π)

        # Turn off uncertain altpolicy b/c when computing the state space system under alternative policies,
        # we want to treat these alternative policies as perfectly credible.
        is_uncertain_altpol = haskey(get_settings(m), :uncertain_altpolicy) && get_setting(m, :uncertain_altpolicy)
        if is_uncertain_altpol
            m <= Setting(:uncertain_altpolicy, false)
        end

        inds = 1:n_states(m) # Don't want augmented states
        Td = Vector{Matrix{S}}(undef, length(prob_vec))
        Cd = Vector{Vector{S}}(undef, length(prob_vec))
        Td[1] = TTT[inds, inds]
        Cd[1] = C[inds]
        for (i, altpolicy) in enumerate(altpolicies)
            j = i + 1
            tmpT, _, tmpC = altpolicy.solve(m) # regime_switching is false, so it is assumed that it is unnecessary here
            Td[j] = tmpT[inds, inds]
            Cd[j] = tmpC[inds]
        end

        # Turn uncertain_altpolicy back on
        if is_uncertain_altpol
            m <= Setting(:uncertain_altpolicy, true)
        end

        T̄ = sum([p .* Td[i] for (i, p) in enumerate(prob_vec)])
        C̄ = sum([p .* Cd[i] for (i, p) in enumerate(prob_vec)])

        Lmat = (Γ2_til * T̄ + Γ0_til)
        Tcal = Lmat \ Γ1_til
        Rcal = Lmat \ Ψ_til
        Ccal = Lmat \ (C_til - Γ2_til * C̄)

        return Tcal, Rcal, Ccal
    end
end

# systems should be a vector of perfectly credible alternative policies, where
# the first element is the implemented one, and the others are alternative policies
# which agents believe may occur.
function gensys_uncertain_altpol(m::AbstractDSGEModel, prob_vec::AbstractVector{S}, regime::Int,
                                 systems::Vector{Union{System, RegimeSwitchingSystem}},
                                 is_altpol::Vector{Bool}) where {S <: Real}

    @assert sum(prob_vec) == 1. "The vector of probabilities must sum to 1"

    @assert length(systems) == length(prob_vec)

    has_pos_prob = findall(x -> x > 0., (@view prob_vec[2:end])) # Note that has_pos_prob has 1 less length than prob_vec

    Γ0, Γ1, C, Ψ, Π = haskey(get_setting(m, :regime_eqcond_info), regime) ?
        get_setting(m, :regime_eqcond_info)[regime].alternative_policy.eqcond(m, regime) :
            eqcond(m, regime)

    Γ0_til, Γ1_til, Γ2_til, C_til, Ψ_til = gensys_to_predictable_form(Γ0, Γ1, C, Ψ, Π)

    # Use systems[1] (b/c first element is perfect cred implemented policy)
    # to infer the initial values of T̅, C̅
    inds = 1:n_states(m)
    T̅, C̅ = (prob_vec[1] == 0.) ? (zeros(size(Γ0_til)), zeros(size(C))) :
        (prob_vec[1] .* (@view systems[1][regime, :TTT][inds, inds]), prob_vec[1] .* (@view systems[1][regime, :CCC][inds]))
    for i in has_pos_prob
        if is_altpol[i]
            T̅ .+= prob_vec[i + 1] * (@view systems[i + 1][:TTT][inds, inds])
            C̅ .+= prob_vec[i + 1] * (@view systems[i + 1][:CCC][inds])
        else
            T̅ .+= prob_vec[i + 1] * (@view systems[i + 1][min(length(systems[i+1].transitions), regime), :TTT][inds, inds])
            C̅ .+= prob_vec[i + 1] * (@view systems[i + 1][min(length(systems[i+1].transitions), regime), :CCC][inds])
        end
    end

    Lmat = (Γ2_til * T̅ + Γ0_til)
    Tcal = Lmat \ Γ1_til
    Rcal = Lmat \ Ψ_til
    Ccal = Lmat \ (C_til - Γ2_til * C̅)

    return Tcal, Rcal, Ccal
end


text:
using SymbolicUtils: Sym, FnType, Term, Add, Mul, Pow, symtype, operation, arguments
using SymbolicUtils
using IfElse: ifelse
using Test

@syms w::Real z::Real a::Real b::Real x

addr = addr = Add(Real, 10, Dict(w=>1, z=>1))

substitute(addr, Dict(w=>1, z=>3.0)) |> display

fieldnames(typeof(addr)) |> display

addr.coeff |> display


text:
module MultiAgentPathFinding

# stdlib
using Random
using DataStructures

# others
using Parameters
using Graphs

# Abstract types
""" MAPFState

Abstract base type for tracking the state of an agent in a MAPF problem
"""
abstract type MAPFState end

""" MAPFAction

Abstract base type for the action that an agent can take in a MAPF problem.
"""
abstract type MAPFAction end

""" MAPFConflict

Abstract base type for a path-path conflict in a MAPF problem.
"""
abstract type MAPFConflict end

""" MAPFConstraints

Abstract base type for the constraints on the low-level path search induced by a conflict in a MAPF problem.
"""
abstract type MAPFConstraints end

""" MAPFEnvironment

Abstract base type for the environment in a MAPF problem.
"""
abstract type MAPFEnvironment end

# Types
export
    MAPFState,
    MAPFAction,
    MAPFConflict,
    MAPFConstraints,
    MAPFEnvironment

# Utils
export
    PlanResult,
    get_mapf_state_from_idx,
    get_mapf_action,
    get_plan_result_from_astar

# High Level Cost elements
export
    HighLevelCost,
    compute_cost,
    accumulate_cost,
    deaccumulate_cost,
    SumOfCosts,
    Makespan

# CBS elements
export
    get_empty_constraint,
    set_low_level_context!,
    get_first_conflict,
    create_constraints_from_conflict,
    overlap_between_constraints,
    add_constraint,
    low_level_search!,
    CBSSolver,
    search!,
    focal_heuristic,
    ECBSSolver

# Grid 2D Types
export
    Grid2DState,
    Grid2DAction,
    Grid2DConflict,
    Grid2DConstraints,
    Grid2DEnvironment,
    Grid2DLocation,
    AStarGrid2DEnvironment

include("utils.jl")
include("high_level_cost.jl")
include("cbs.jl")
include("ecbs.jl")
include("domains/grid2d/types.jl")
include("domains/grid2d/cbs_grid2d.jl")
include("domains/grid2d/ecbs_grid2d.jl")
include("domains/grid2d/a_star_epsilon_grid2d.jl")

end # module


text:
# Admission probability functions

abstract type AbstractAdmProbFctSwitches{F1 <: Real} <: ModelSwitches end
abstract type AbstractAdmProbFct{F1 <: Real} <: ModelObject end

"""
	$(SIGNATURES)

Initialize admission probability functions for all colleges (as a set).
Still need to make the function for each college.
"""
function init_admprob_fct(switches :: AbstractAdmProbFctSwitches{F1}) where F1 end

n_colleges(switches :: AbstractAdmProbFctSwitches{F1}) where F1 =
    switches.nColleges;

Lazy.@forward AbstractAdmProbFct.switches (
    ModelParams.has_pvector, n_colleges, n_open_colleges,
    by_college, by_college!, college_index, param_names
    );

Base.show(io :: IO, af :: AbstractAdmProbFct{F1}) where F1 =
    print(io, typeof(af));
Base.show(io :: IO, af :: AbstractAdmProbFctSwitches{F1}) where F1 =
    print(io, typeof(af));

ModelParams.has_pvector(switches :: AbstractAdmProbFctSwitches{F1}) where F1 = true;

ModelParams.get_object_id(switches :: AbstractAdmProbFctSwitches{F1}) where F1 =
    get_object_id(switches.pvec);

# function prob_admit(admProbFct :: AF1,
#     iCollege :: Integer, hsGpa) where {AF1 <: AbstractAdmProbFct{<: Real}, I1, F1}

#     # This is expensive, but hard to avoid. Need to ensure that current parameters are used when admission prob fct is constructed.
#     prob_fct = make_admprob_function(admProbFct, iCollege);
#     return prob_fct.(hsGpa)
# end

function make_admprob_fct_switches(objId :: ObjectId, switchType :: DataType,
    nc, nOpen, byCollegeV, calibratedV)

    pvecV = Vector{Param}();
    for pName in param_names(switchType)
        nColl = pvec_length(nc, nOpen, byCollegeV, pName);
        isCalibrated = any(isequal.(pName, calibratedV));
        push!(pvecV, init_param(switchType, pName, nColl, isCalibrated));
    end
    pvec = ParamVector(objId, pvecV);

    switches = switchType(pvec, nc, nOpen);
    @assert validate_admprob_fct_switches(switches);
    return switches
end

# Returns the no of colleges with parameters for a parameter.
# Exammple: `pvec_length(4, 1, [:mV], :bV) == 1`.
function pvec_length(nc, nOpen, byCollegeV, pName)
    if any(isequal.(byCollegeV, pName))
        return nc - nOpen;
    else
        return 1;
    end
end

## ------------  Access

# Index into parameters that vary by college.
function college_index(switches :: AbstractAdmProbFctSwitches{F1}, iCollege) where F1
    return iCollege .- n_open_colleges(switches);
end

"""
	$(SIGNATURES)

Does the parameter `pName` vary by college?
"""
function by_college(switches :: AbstractAdmProbFctSwitches{F1}, pName) where F1
    v = param_value(switches, pName);
    @assert !isnothing(v)  "Not found: $pName";
    return length(v) > 1
end

"""
	$(SIGNATURES)

Switch a parameter to vary by college.
"""
function by_college!(switches :: AbstractAdmProbFctSwitches{F1}, pName) where F1
    set_by_college!(switches, pName, true);
end

"""
	$(SIGNATURES)

Switch a parameter NOT to vary by college.
"""
function not_by_college!(switches :: AbstractAdmProbFctSwitches{F1}, pName) where F1
    set_by_college!(switches, pName, false);
end

function set_by_college!(switches :: AbstractAdmProbFctSwitches{F1},
    pName, byCollege :: Bool) where F1

    # Check if we need to make a change
    (by_college(switches, pName) == byCollege)  &&  (return nothing);
    if byCollege
        nColl = n_colleges(switches) - n_open_colleges(switches);
    else
        nColl = 1;
    end
    isCalibrated = is_calibrated(switches, pName);
    p = init_param(typeof(switches), pName, nColl, isCalibrated);
    ModelParams.replace!(get_pvector(switches), p);
end

n_open_colleges(switches :: AbstractAdmProbFctSwitches{F1}) where F1 =
    switches.nOpenColleges;

"""
	$(SIGNATURES)

Get a parameter that may or may not vary by college for college `iCollege`.
For college 1, throw an error.
"""
function get_param(af :: AbstractAdmProbFct{F1}, pName,
    iCollege :: Integer) where F1

    @assert (iCollege > n_open_colleges(af))  "No parameters for college $iCollege";
    if by_college(af, pName)
        # There are no parameters for the first college
        idx = college_index(af, iCollege);
    else
        idx = 1;
    end
    return getproperty(af, pName)[idx]
end

get_params(af :: AbstractAdmProbFct{F1}, iCollege :: Integer) where F1 =
    [get_param(af, varName, iCollege)  for varName in param_names(af)];

include("admprob_fct_open.jl");
include("admprob_fct_logistic.jl");
include("admprob_fct_linear.jl");
include("admprob_fct_step.jl");

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


text:
__precompile__(true)

module SampleDistributions

using Distributions
using StatsBase
using NamedTuples

include("utils.jl")
include("sampledistribution.jl")
include("distributionsapi.jl")
include("informationtheory.jl")

end # module


text:
addprocs(3)

using TransformUtils, Gadfly
using KernelDensityEstimate
using IncrementalInference
using RoME
using Caesar
using Distributions
# using Base.Test

# import RoME: ⊕

# Parametric solution with hand checked solution as ground truth
# XYZ
gt = Dict{Symbol, Tuple{Symbol,Vector{Float64}}}()
gt[:l1] =(:XYZ ,  [2.07618; -2.85304; 2.41059] )
gt[:l2] =(:XYZ ,  [2.27909; -3.53382; 1.5933] )
gt[:l3] =(:XYZ ,  [1.83487; -3.26179; 2.38189] )
gt[:l4] =(:XYZ ,  [1.85991; -3.26293; 2.77598] )
gt[:l5] =(:XYZ ,  [2.11136; -2.86976; 2.72259] )
gt[:l6] =(:XYZ ,  [1.79027; -3.22242; 2.15042] )
gt[:l7] =(:XYZ ,  [2.02859; -2.82187; 2.16534] )

# XYZ qWqXqYqZ
gt[:x1]= (:XYZqWXYZ , [0.000229458; 0.299991; -0.00230474;   0.696163;  0.0546845;  -0.123042;  0.705144] )
gt[:x2]= (:XYZqWXYZ , [-0.314997; -1.1901; 1.73293;          0.943023; -0.0847915;  -0.102545; -0.304964] )
gt[:x3]= (:XYZqWXYZ , [-0.363944; -1.18799; 1.74864;         0.937102; -0.0847783;  -0.108046; -0.320903] )
gt[:x4]= (:XYZqWXYZ , [-0.262489; -0.865555; 1.82557;        0.935268; -0.0797656; -0.0966424; -0.331015] )
gt[:x5]= (:XYZqWXYZ , [-0.680204; -1.37827; 1.73447;         0.949729; -0.0644836;  -0.120014; -0.281874] )
gt[:x6]= (:XYZqWXYZ , [-0.614131; -1.29263; 1.75964;         0.951792; -0.0638253;  -0.117314; -0.276146] )

"""
    gt2fg!(gt)

Insert into existing or new FactorGraph type from ground truth data in gt.

Example

  gt[:l7] = (:XYZ ,  [2.02859; -2.82187; 2.16534] )
  gt[:x1] = (:XYZqWXYZ , [0.000229458; 0.299991; -0.00230474;   0.696163;  0.0546845;  -0.123042;  0.705144] )
"""
function gt2fg!(gt::Dict{Symbol, Tuple{Symbol,Vector{Float64}}};
    fg=nothing,
    session="GROUNDTRUTH" )
  #
  fg = fg != nothing ? fg : initfg(sessionname=session)

  for (key, val) in fg
    #
  end

  warn("Under construction")
  fg
end

# setup visualization process and default drawings
vc = startdefaultvisualization()

N = 200
fg = Caesar.initfg()

sqrtinv = [[10;0;0;0;0;0]';
[0;10;0;0;0;0]';
[0;0;10;0;0;0]';
[0;0;0;57.2958;0;0]';
[0;0;0;0;57.2958;0]';
[0;0;0;0;0;57.2958]'];

initCov = inv(sqrtinv^2)
odoCov = deepcopy(initCov)

sqrtinv = [[286.479;0]';
[0;200]']
brCov = inv(sqrtinv^2)

println("Adding PriorPose3 to graph...")
init = MvNormal([0.0;0.0;0.0;-0.00772052; 0.0; -0.0992321],initCov)
v1 = addNode!(fg, :x1,  zeros(6,N) ⊕ Pose3Pose3(init ),  N=N)
initPosePrior = PriorPose3(init)
f1  = addFactor!(fg,[v1], initPosePrior)

# Odometry pose 0 pose 1: (-0.380711, -1.02585, 1.77348; -2.19796, -0.151721, -0.0929671)
# odo = SE3([-0.380711, -1.02585, 1.77348], Euler(-0.0929671, -0.151721, -2.19796) )
odo = MvNormal([-0.380711; -1.02585; 1.77348;-0.0929671; -0.151721; -2.19796], odoCov )
addOdoFG!(fg, Pose3Pose3(deepcopy(odo)) )

# Odometry pose 1 pose 2: (-0.00138117, -0.0600476, 0.0204065; -0.038278, 0.0186151, 0.00606331)
odo = MvNormal([-0.00138117; -0.0600476; 0.0204065; 0.00606331; 0.0186151; -0.038278], odoCov )
addOdoFG!(fg, Pose3Pose3(odo) )

# Odometry pose 2 pose 3: (-0.298704, -0.113974, 0.0244868; -0.00805163, -0.00425057, -0.0275746)
odo = MvNormal([-0.298704; -0.113974; 0.0244868; -0.0275746; -0.00425057; -0.00805163], odoCov )
addOdoFG!(fg, Pose3Pose3(deepcopy(odo)) )

# Odometry pose 3 pose 4: (0.670492, 0.0251882, -0.0705555; 0.109113, 0.0228966, 0.039246)
odo = MvNormal([0.670492; 0.0251882; -0.0705555; 0.039246; 0.0228966; 0.109113], odoCov )
addOdoFG!(fg, Pose3Pose3(deepcopy(odo)) )

# Odometry pose 4 pose 5: (-0.0615969, 0.00849829, -0.0170747; 0.0315662, -0.0115416, -0.00605423)
odo = MvNormal([-0.0615969; 0.00849829; -0.0170747;  -0.00605423;-0.0115416;0.0315662], odoCov )
addOdoFG!(fg, Pose3Pose3(deepcopy(odo)) )

println("Adding landmarks to graph...")
# Bearing Range pose 1 landmark 0: (0.0103329, 2.99139)
# Bearing Range pose 2 landmark 0: (0.0456434, 3.02145)
# Bearing Range pose 3 landmark 0: (-0.0277112, 3.12405)
# Bearing Range pose 4 landmark 0: (0.0778371, 3.20003)
# Bearing Range pose 5 landmark 0: (0.03676, 3.17675)

rangecov, bearingcov=3e-5, 3e-5

# vl1 = addNode!(fg, :l1,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (2.9914, 0.010333), :x2, :l1, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.02145, 0.04564), :x3, :l1, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.12405, -0.027711), :x4, :l1, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.2, 0.0778371), :x5, :l1, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.17675, 0.03676), :x6, :l1, rangecov=rangecov,bearingcov=bearingcov)

# Bearing Range pose 1 landmark 1: (-0.0971847, 3.49785)
# Bearing Range pose 2 landmark 1: (-0.0522769, 3.5344)
# Bearing Range pose 3 landmark 1: (-0.115672, 3.69064)
# Bearing Range pose 4 landmark 1: (-0.0471636, 3.66387)
# Bearing Range pose 5 landmark 1: (-0.0887547, 3.66507)
# addNode!(fg, :l2,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (3.49785, -0.0971847), :x2, :l2, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.5344, -0.0522769), :x3, :l2, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.69064, -0.115672), :x4, :l2, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.66387, -0.0471636), :x5, :l2, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.66507, -0.0887547), :x6, :l2, rangecov=rangecov,bearingcov=bearingcov)

# Bearing Range pose 1 landmark 2: (-0.149587, 3.06005)
# Bearing Range pose 2 landmark 2: (-0.0960879, 3.07682)
# Bearing Range pose 3 landmark 2: (-0.176641, 3.23283)
# Bearing Range pose 4 landmark 2: (-0.065336, 3.21081)
# Bearing Range pose 5 landmark 2: (-0.115751, 3.20317)
# addNode!(fg, :l3,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (3.06005, -0.149587), :x2, :l3, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.07682, -0.0960879), :x3, :l3, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.23283, -0.176641), :x4, :l3, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.21081, -0.065336), :x5, :l3, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.20317, -0.115751), :x6, :l3, rangecov=rangecov,bearingcov=bearingcov)

# Bearing Range pose 1 landmark 3: (-0.151731, 3.18027)
# Bearing Range pose 2 landmark 3: (-0.099368, 3.20953)
# Bearing Range pose 3 landmark 3: (-0.174175, 3.33787)
# Bearing Range pose 4 landmark 3: (-0.0662606, 3.32899)
# Bearing Range pose 5 landmark 3: (-0.114886, 3.32137)
# addNode!(fg, :l4,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (3.18027, -0.151731),:x2, :l4, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.20953, -0.099368), :x3, :l4, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.33787, -0.174175),:x4, :l4, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.32899, -0.0662606),:x5, :l4, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.32137, -0.114886),:x6, :l4, rangecov=rangecov,bearingcov=bearingcov)

# Bearing Range pose 1 landmark 4: (0.00322788, 3.10878)
# Bearing Range pose 2 landmark 4: (0.0370399, 3.15066)
# Bearing Range pose 3 landmark 4: (-0.030071, 3.22893)
# Bearing Range pose 4 landmark 4: (0.0683279, 3.31492)
# Bearing Range pose 5 landmark 4: (0.0323075, 3.29366)
# addNode!(fg, :l5,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (3.10878, 0.00322788),:x2, :l5, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.15066, 0.0370399), :x3, :l5, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.22893, -0.030071), :x4, :l5, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.31492, 0.0683279), :x5, :l5, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.29366, 0.0323075), :x6, :l5, rangecov=rangecov,bearingcov=bearingcov)

# Bearing Range pose 1 landmark 5: (-0.146491, 2.95272)
# Bearing Range pose 2 landmark 5: (-0.090717, 2.99605)
# Bearing Range pose 3 landmark 5: (-0.178344, 3.14063)
# Bearing Range pose 4 landmark 5: (-0.0603386, 3.10463)
# Bearing Range pose 5 landmark 5: (-0.111392, 3.10912)
# addNode!(fg, :l6,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (2.95272, -0.146491), :x2, :l6, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (2.99605, -0.090717), :x3, :l6, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.14063, -0.178344), :x4, :l6, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.10463, -0.0603386), :x5, :l6, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.10912, -0.111392), :x6, :l6, rangecov=rangecov,bearingcov=bearingcov)

# Bearing Range pose 1 landmark 6: (0.0179481, 2.8873)
# Bearing Range pose 2 landmark 6: (0.0543992, 2.9317)
# Bearing Range pose 3 landmark 6: (-0.0250711, 3.03224)
# Bearing Range pose 4 landmark 6: (0.0878667, 3.09855)
# Bearing Range pose 5 landmark 6: (0.041406, 3.07297)
# addNode!(fg, :l7,  0.1*randn(3,N),  N=N)
addLinearArrayConstraint(fg, (2.8873, 0.0179481), :x2, :l7, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (2.9317, 0.0543992), :x3, :l7, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.03224, -0.0250711), :x4, :l7, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.09855, 0.0878667), :x5, :l7, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (3.07297, 0.041406), :x6, :l7, rangecov=rangecov,bearingcov=bearingcov)

# solving = [true]
# @async begin
#   while solving[1]
#     println("visualizing")
#     visualizeallposes!(vc, fg)
#     visualizeDensityMesh!(vc, fg, :l1)
#     sleep(2)
#   end
# end
#
# tree = wipeBuildNewTree!(fg,drawpdf=false);
# @time inferOverTree!(fg, tree)
# solving[1]=false;

batchSolve(fg)
visualize(fg, vc, N=N, drawtype=:fit, gt=gt)

visualizeallposes!(vc, fg, drawtype=:fit, gt=gt)

# ls(fg)

visualizeDensityMesh!(vc, fg, :l6)

[inferOverTree!(fg, tree) for i in 1:3]; # should not be required

draw(PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
 plotKDE( getVertKDE(fg,:l1), dimLbls=["x";"y";"z";"phi";"the";"psi"]) )

 draw(PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
  plotKDE( pl1, dimLbls=["x";"y";"z";"phi";"the";"psi"]) )

draw(PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
 plotKDE( getVertKDE(fg,:x6), dimLbls=["x";"y";"z";"phi";"the";"psi"]) )

XX = String["x$i" for i in 1:6];
LL = String["l$i" for i in 1:7];
[draw(PDF("/home/dehann/Desktop/imgs/$(i).pdf",30cm,20cm),
  plotKDE( getVertKDE(fg,i), dimLbls=["x";"y";"z";"phi";"the";"psi"]) ) for i in union(XX,LL)];

# # new install
# Pkg.add("KernelDensityEstimate")
#
# # or update
# Pkg.update()
#
# # draw this
# using KernelDensityEstimate, Gadfly
#
# p, q = kde!(randn(2,100)), kde!(randn(2,100)+2);
#
# plotKDE(p)
# plotKDE([p;q],c=["red";"blue"],levels=3)
#
# p, q = kde!(randn(4,100)), kde!(randn(4,100)+2);
# draw(PDF("test.pdf",30cm,20cm),
#         plotKDE( [p;q],c=["red";"blue"],levels=3,
# 				dimLbls=["x";"y";"z";"phi";"the";"psi"]) )

#