text:
function Base.show(io::IO, sol::FODESolution)
    println("Time span $(typeof(sol.t))")
    println(sol.t)
    println("Solution $(typeof(sol.u))")
    print(sol.u)
end

function Base.show(io::IO, prob::SingleTermFODEProblem)
    printstyled("SingleTermFODEProblem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.α)", color=:red)
    println()
    println("timespan: $(prob.tspan)")
    println("u0: $(prob.u0)")
end

function Base.show(io::IO, prob::MultiTermsFODEProblem)
    printstyled("MultiTermsFODEProblem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.orders)", color=:red)
    println()
    println("timespan: $(prob.tspan)")
    println("u0: $(prob.u0)")
end

function Base.show(io::IO, prob::FODESystem)
    printstyled("FODESystem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.α)", color=:red)
    println()
    println("timespan: $(prob.tspan)")
    println("u0: $(prob.u0)")
end

function Base.show(io::IO, prob::FDDEProblem)
    printstyled("FDDEProblem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.α)", color=:red)
    println()
    println("timespan: $(prob.tspan)")
    println("history function: $(prob.ϕ)")
end

function Base.show(io::IO, prob::FDDESystem)
    printstyled("FDDESystem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.α)", color=:red)
    println()
    println("timespan: $(prob.T)")
    println("history function: $(prob.ϕ)")
end

function Base.show(io::IO, prob::DODEProblem)
    printstyled("DODEProblem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.orders)", color=:red)
    println()
    println("timespan: $(prob.tspan)")
end

function Base.show(io::IO, prob::FractionalDifferenceProblem)
    printstyled("FractionalDifferenceProblem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.α)", color=:red)
    println()
    println("u0: $(prob.u0)")
end

function Base.show(io::IO, prob::FIEProblem)
    printstyled("FIEProblem", color=:light_blue)
    printstyled(" with order ")
    printstyled("$(prob.orders)", color=:red)
    println()
    println("timespan: $(prob.tspan)")
end

"""
Fractional differential equation solutions visulization hooks.
"""
@recipe f(sol::FODESolution) = sol.t, sol.u

@recipe f(sol::FDifferenceSolution) = sol.t, sol.u

@recipe f(sol::FIESolution) = sol.t, sol.u

@recipe f(sol::DODESolution) = sol.t, sol.u


text:
function champs_3(t, fft_ubar, fft_vbar, fft_ug, fft_vg,
                  A1, A2, matr, conjmatr, sigma, llambda,
                  Ktaubis, epsilon, Ntaumm)

    h1u, h1v = ftau(t, fft_ubar, fft_vbar,
                    A1, A2, matr, conjmatr,
                    sigma, llambda)

    h1u_fft = fft(h1u, dims=1)
    h1v_fft = fft(h1v, dims=1)

    dtu0u = h1u_fft[1, :] / Ntaumm
    dtu0v = h1v_fft[1, :] / Ntaumm

    h1u_fft[1, :] .= 0
    h1v_fft[1, :] .= 0

    h1u = epsilon * ifft(h1u_fft ./ (1im * Ktaubis), dims=1)
    h1v = epsilon * ifft(h1v_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
			      fft_ubar, fft_vbar,
			      dtu0u, dtu0v,
                              A1, A2,
			      matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar, fft_vbar,
			      A1, A2,
			      matr, conjmatr,
                              sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2
    champv = champv1 .+ champv2

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0.0
    champv_fft[1, :] .= 0.0

    dth1u = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth1v = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    h2u, h2v = ftau(t,
		    fft_ubar .+ h1u,
                    fft_vbar .+ h1v,
                    A1, A2, matr,
                    conjmatr, sigma, llambda)

    h2u = h2u .- dth1u
    h2v = h2v .- dth1v

    h2u_fft = fft(h2u, dims=1)
    h2v_fft = fft(h2v, dims=1)

    h2u_fft[1, :] .= 0.0
    h2v_fft[1, :] .= 0.0

    dtauh2u = ifft(h2u_fft, dims=1)
    dtauh2v = ifft(h2v_fft, dims=1)

    h2u = epsilon * ifft(h2u_fft ./ (1im * Ktaubis), dims=1)
    h2v = epsilon * ifft(h2v_fft ./ (1im * Ktaubis), dims=1)

    C2u = fft_ubar .+ h2u
    C2v = fft_vbar .+ h2v

    ffu, ffv = ftau(t,
		    C2u .+ fft_ug,
		    C2v .+ fft_vg,
                    A1, A2,
		    matr, conjmatr,
                    sigma, llambda)

    champu, champv = ftau(t,
			  C2u, C2v,
			  A1, A2,
			  matr, conjmatr,
			  sigma, llambda)

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champubaru = champu_fft[1, :] / Ntaumm
    champubarv = champv_fft[1, :] / Ntaumm

    champu1, champv1 = duftau(t, fft_ubar, fft_vbar,
                              champubaru, champubarv,
                              A1, A2, matr, conjmatr,
                              sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t, fft_ubar, fft_vbar,
                              A1, A2, matr, conjmatr,
                              sigma, llambda, Ntaumm)

    champu_fft = fft(champu1 .+ champu2, dims=1)
    champv_fft = fft(champv1 .+ champv2, dims=1)

    dttu0baru = champu_fft[1, :] / Ntaumm
    dttu0barv = champv_fft[1, :] / Ntaumm

    champu_fft[1, :] .= 0.0
    champv_fft[1, :] .= 0.0

    dth1baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dth1barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = d2uftau(t, fft_ubar, fft_vbar, dtu0u, dtu0v,
                               champubaru, champubarv, A1, A2, matr,
                               conjmatr, sigma, llambda, Ntaumm)

    champu2, champv2 = duftau(t, fft_ubar, fft_vbar, dttu0baru, dttu0barv,
                              A1, A2, matr, conjmatr, sigma, llambda,
                              Ntaumm)

    champu3, champv3 = dtduftau(t, fft_ubar, fft_vbar, dtu0u, dtu0v,
                                A1, A2, matr, conjmatr, sigma, llambda,
				Ntaumm)

    champu4, champv4 = dtduftau(t, fft_ubar, fft_vbar,
				champubaru, champubarv,
                                A1, A2, matr, conjmatr, sigma, llambda,
                                Ntaumm)

    champu5, champv5 = d2tftau(t, fft_ubar, fft_vbar, A1, A2,
			       matr, conjmatr,
                               sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2 .+ champu3 .+ champu4 .+ champu5
    champv = champv1 .+ champv2 .+ champv3 .+ champv4 .+ champv5

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0
    champv_fft[1, :] .= 0

    dtth1baru = epsilon * ifft(champu_fft ./ (1im * Ktaubis), dims=1)
    dtth1barv = epsilon * ifft(champv_fft ./ (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
                              champubaru .+ dth1baru,
			      champubarv .+ dth1barv,
                              A1, A2, matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t,
			      fft_ubar .+ h1u,
			      fft_vbar .+ h1v,
			      A1, A2,
                              matr, conjmatr,
			      sigma, llambda, Ntaumm)

    champu = champu1 .+ champu2 .- dtth1baru
    champv = champv1 .+ champv2 .- dtth1barv

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champu_fft[1, :] .= 0.0
    champv_fft[1, :] .= 0.0

    dtC2u = champubaru .+ epsilon * ifft(champu_fft ./ (1im * Ktaubis),
					 dims=1)
    dtC2v = champubarv .+ epsilon * ifft(champv_fft ./ (1im * Ktaubis),
					 dims=1)

    champgu = ffu .- dtauh2u .- dtC2u
    champgv = ffv .- dtauh2v .- dtC2v

    champgu_fft = fft(champgu, dims=1)
    champgv_fft = fft(champgv, dims=1)

    champmoyu = champgu_fft[1, :] / Ntaumm
    champmoyv = champgv_fft[1, :] / Ntaumm

    champgu_fft[1, :] .= 0.0
    champgv_fft[1, :] .= 0.0

    ichampgu = ifft(champgu_fft ./ (1im * Ktaubis), dims=1)
    ichampgv = ifft(champgv_fft ./ (1im * Ktaubis), dims=1)

    champubaru, champubarv, ichampgu, ichampgv, champmoyu, champmoyv

end

ext: jl


text:
verbose = true
 if verbose && !isdefined(Main,:RvSpectML)   println("# Loading RvSpecML")    end
 using RvSpectML

# USER: You must create a data_paths.jl file in one of the default_paths_to_search listed below. It need only contain one line:
# expres_data_path = "/path/to/EXPRES/data/not/including/target_subdir"
target_subdir = "101501"   # USER: Replace with directory of your choice
 fits_target_str = "101501"
 output_dir = "examples/output/"
 default_paths_to_search = [pwd(),"examples",joinpath(pkgdir(RvSpectML),"examples"),"/gpfs/group/ebf11/default/ebf11/expres/inputs"]
 # NOTE: make_manifest does not update its paths_to_search when default_paths_to_search is defined here, so if you change the line above, you must also include "paths_to_search=default_paths_to_search" in the make_manifest() function call below
 pipeline_plan = PipelinePlan()
 dont_make_plot!(pipeline_plan, :movie)

reset_all_needs!(pipeline_plan)
 if need_to(pipeline_plan,:read_spectra)
   if verbose println("# Finding what data files are avaliable.")  end
   df_files = make_manifest(target_subdir, EXPRES )

   if verbose println("# Reading in customized parameters from param.jl.")  end
   eval(code_to_include_param_jl())

   if verbose println("# Reading in ", size(df_files_use,1), " FITS files.")  end
   @time all_spectra = map(EXPRES.read_data,eachrow(df_files_use))
   RvSpectML.discard_blaze(all_spectra)
   RvSpectML.discard_continuum(all_spectra)
   dont_need_to!(pipeline_plan,:read_spectra)
   #all_spectra
 end


text:
passwordMatcher = r"(\d*)-(\d*) (.): (\w*)"

validPasswords = 0
for line in eachline()
    result = match(passwordMatcher, line)
    min = parse(Int32, result[1])
    max = parse(Int32, result[2])

    count = 0
    for match in eachmatch(Regex(result[3]), result[4])
        count += 1
    end

    if count >= min && count <= max
        global validPasswords += 1
    end
end

println(validPasswords)


text:
if !@isdefined(TestEdge)
mutable struct TestEdge <: AEdge
    src::Int
    dst::Int
end
Erdos.src(e::TestEdge) = e.src
Erdos.dst(e::TestEdge) = e.dst
end

@testset "$TEST $G" begin

# EdgeMap matrix
n = 4
m = randn(n, n)
m = (m + m'); m -= Diagonal(m)

g = G(m)
em = EdgeMap(g, m)
@test length(em) == n*n
@test em[TestEdge(1,2)] == em[1,2]
@test values(em) === m
@test valtype(em) == Float64
# @test sprint(show, em) == "EdgeMap{Float64}" #TODO

d = [1 0 1;
     1 0 2;
     3 0 0]
g = G(d, upper=true)
m = EdgeMap(g, d)
@test m[1,2] == m[2,1] == 0
@test m[1,3] == m[3,1] == 1
@test m[1,1] == 1
@test m[2,2] == 0
m[3, 1] = -1
@test m[1,3] == m[3,1] == -1
@test d[1,3] == -1
@test d[3,1] == 3

d = [1 0 1;
     1 0 2;
     3 0 0]
g = DG(d)
m = EdgeMap(g, d)
@test m[1,2] != m[2,1]
@test m[1,3] != m[3,1]
@test m[1,1] == 1
@test m[2,2] == 0
m[3, 1] = -1
@test m[1,3]  == 1
@test m[3,1]  == -1
@test d[1,3] == 1
@test d[3,1] == -1

@test get(m, edge(g, 1, 2), -100) == 0

m = EdgeMap(g, sparse(d))
@test values(m) == nonzeros(sparse(d))
@test get(m, edge(g, 1, 2), -100) == 0

# ConstEdgeMap
em = ConstEdgeMap(g, 1)
@test valtype(em) == Int
@test length(em) == typemax(Int)
# @test_throws ErrorException em[1,2] = 18 #TODO
# @test_throws ErrorException em[Edge(1,2)] = 18
# @test_throws ErrorException em[E(1,2)] = 18
em[1,2] = 18
em[Edge(1,2)] = 18
em[E(1,2)] = 18

@test em[1,2] == 1
@test values(em) == [1]
# @test sprint(show, em) == "EdgeMap{Int64}"

n = 10
g = G(n)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3)

MAPS_LIST = [  EdgeMap(g, Int), EdgeMap(g, spzeros(Int, n, n)), EdgeMap(g, zeros(Int, n, n)),
            EdgeMap(g, Dict{Edge{Int}, Int}()), EdgeMap(g, Dict{E, Int}()) ]
E <: IndexedEdge && append!(MAPS_LIST, [EdgeMap(g, zeros(Int, ne(g))), EdgeMap(g, Dict{Int,Int}())])

for m in MAPS_LIST
    for (k, e)  in enumerate(edges(g))
        i, j = src(e), dst(e)
        m[e] = k
        @test k == m[e]
        @test k == m[i, j]
        @test get(m, e, -1) == k
        @test haskey(m, e)
        @test haskey(m, i, j)
    end
    m[1,2] = 10
    @test m[1,2] == 10
    # @test sprint(show, m) == "EdgeMap: $(m.data)"
    if typeof(m.data) <: AbstractMatrix
        @test haskey(m, 3, 4)
    else
        @test !haskey(m, 3, 4)
    end
    @test !haskey(m, 100,100)
end

g = G(10,10)
m = EdgeMap(g, rand(10,10))
adj = adjacency_list(g)
madj = edgemap2adjlist(m)
@test length(madj) == length(adj)
for i=1:nv(g)
    @test length(madj[i]) == length(adj[i])
    for (k, j) in enumerate(adj[i])
        @test madj[i][k] == m[i, j]
    end
end

g = G(10,20)
emap = EdgeMap(g, rand(10,20))
m = Matrix(emap)
@test typeof(m) == Matrix{Float64}
sp = sparse(emap)
@test typeof(sp) <: SparseMatrixCSC{Float64}

@test size(m) == size(sp) == (nv(g), nv(g))
@test count(!iszero,m) == count(!iszero,sp) == 2ne(g)

for e in edges(g)
    u, v = src(e), dst(e)
    @test m[u,v] == m[v,u] == sp[u, v] == sp[v, u] == emap[e]
end

@test size(emap) == (nv(g), nv(g))
@test size(emap, 1) == nv(g)
@test size(emap, 2) == nv(g)
@test_throws ErrorException size(emap, 3)

g = DG(10,20)
emap = EdgeMap(g, rand(10,20))
m = Matrix(emap)
@test typeof(m) == Matrix{Float64}
sp = sparse(emap)
@test typeof(sp) <: SparseMatrixCSC{Float64}

@test size(m) == size(sp) == (nv(g), nv(g))
@test count(!iszero,m) == count(!iszero,sp) == ne(g)

for e in edges(g)
    u, v = src(e), dst(e)
    @test m[u,v] == sp[u, v] == emap[e]
end
@test valtype(emap) == Float64

@test size(emap) == (nv(g), nv(g))

g = G(10,20)
m = EdgeMap(g, e -> src(e) + dst(e))
if E <: IndexedEdge
    @test typeof(m.data) == Vector{V}
else
    @test typeof(m.data) == Dict{E, V}
end
for e in edges(g)
    @test m[e] == src(e) + dst(e)
end
@test valtype(m) == V

m = EdgeMap(g, e -> norm(rand(2)))
@test valtype(m) == Float64
m = EdgeMap(g, e -> rand(2))
@test valtype(m) == Vector{Float64}

end # testset


text:
#=
Author : Guy Farmer
generates a distribution of ambulances for various ambulance counts and stores to a jld file
=#
using AmbulanceDeployment

#using AmbulanceDeployment
using DataFrames, Winston, JLD, CSV, Gurobi, JuMP, GLPK, JLD, JSON

# #I CAN'T for the life of me see the difference.
adjacent_nbhd = CSV.File("C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy/test/austin_data/adjacent_nbhd.csv") |> DataFrame;
coverage = CSV.File("C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy/test/austin_data/coverage.csv") |> DataFrame;
#coverage = convert(Array{Bool, 2}, coverage[:, :])
hourly_calls = CSV.File("C:/Users/Owner/Documents/Austin/AmbulanceDeployment/AmbulanceDeployment.jl-legacy/test/austin_data/Full_WeekdayCalls.csv") |> DataFrame;

num_regions = 209;
regions = collect(1:209);
locations = collect(1:size(coverage,2));
coverage = Matrix(coverage);
coverage = convert(Array{Bool, 2}, coverage[:, :]);
adjacent = Matrix(adjacent_nbhd[!,2:end]);
adjacent = adjacent[regions,regions] .> 0.5;
demand = Matrix(hourly_calls[!,5:end-1]); #we were trimming first 6 collumns but it seems to me it only the first contain ex

#= previously there was an error that didnt account for all the regions because it
 it only included regions for incidents. the regions without incidents are changed
 to 0=#
# for x in incidents[!,:neighborhood]
#     if(!haskey(regions2index,x))
#         regions2index[x] = 0
#     end
# end

# We focus on emergency calls during the "peak period" (8AM - 8PM),
# with the emergency calls from the first 3 month as our training set,
# and the subsequent emergency calls from the remaining months as our test set
peak_period = (hourly_calls[!,:hour] .>= 8) .* (hourly_calls[!,:hour] .<= 20);
indices = 1:DataFrames.nrow(hourly_calls);
train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 3);
test_filter  = .~train_filter;
train_indices = indices[train_filter];
test_indices = indices[test_filter];

p = DeploymentProblem(30, length(locations), length(regions), demand, train_indices,
                      test_indices, coverage[regions,:], Array{Bool,2}(adjacent));

amb_deployment = Dict{Symbol, Dict{Int, Vector{Int}}}();

(deployment_model, name) = (next_dp -> StochasticDeployment(next_dp, nperiods=500), :Stochastic)
println("$name: ")
amb_deployment[name] = Dict{Int, Vector{Int}}()
namb = 40
println("$namb ")
p.nambulances = namb
println("time for model generation - stochastic $namb ambulances")
## make sure to include @time begin
@time model = deployment_model(p);
#@time model2 = StochasticDeployment_hyp(p)
set_optimizer(model.m, Gurobi.Optimizer);
println("time for model solution - stochastic $namb ambulances")
@time AmbulanceDeployment.optimize!(model, p)
#@time AmbulanceDeployment.StochasticDeployment_hyp!(model, p, extra_amb = 1)

#deployment(m::StochasticDeployment) = [round(Int,x) for x in JuMP.getValue(m.x)]
#deployment2(m::StochasticDeployment_hyp) = [round(Int,x) for x in JuMP.getValue(m.x)]
amb_deployment[name][namb] = deployment(model)

#[0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0, -0.0, -0.0, 0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.0, 0.0]

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
# ---
# title: 920. Number of Music Playlists
# id: problem920
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-music-playlists/description/>
# hidden: true
# ---
#
# Your music player contains `N` different songs and she wants to listen to `L`
# **** (not necessarily different) songs during your trip.  You create a
# playlist so that:
#
#   * Every song is played at least once
#   * A song can only be played again only if `K` other songs have been played
#
# Return the number of possible playlists.  **As the answer can be very large,
# return it modulo`10^9 + 7`**.
#
#
#
# **Example 1:**
#
#
#
#     Input: N = 3, L = 3, K = 1
#     Output: 6
#     **Explanation** : There are 6 possible playlists. [1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1].
#
#
# **Example 2:**
#
#
#
#     Input: N = 2, L = 3, K = 0
#     Output: 6
#     **Explanation** : There are 6 possible playlists. [1, 1, 2], [1, 2, 1], [2, 1, 1], [2, 2, 1], [2, 1, 2], [1, 2, 2]
#
#
# **Example 3:**
#
#
#
#     Input: N = 2, L = 3, K = 1
#     Output: 2
#     **Explanation** : There are 2 possible playlists. [1, 2, 1], [2, 1, 2]
#
#
#
#
# **Note:**
#
#   1. `0 <= K < N <= L <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "Fast scanning" begin
  # 2-D array
  img = zeros(Float64, 9, 9)
  img[4:6,:] .= 10
  img[:,4:6] .= 11

  seg = fast_scanning(img, 2)
  expected = fill(4,(9,9))
  expected[1:3,1:3] .= 1
  expected[1:3,7:9] .= 5
  expected[7:9,1:3] .= 3
  expected[7:9,7:9] .= 6
  expected_labels = [1,3,4,5,6]
  expected_means = Dict(1=>0.0,3=>0.0,4=>10.6,5=>0.0,6=>0.0)
  expected_count = Dict(1=>9,3=>9,4=>45,5=>9,6=>9)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  # Offset array
  img = centered(img)

  seg = fast_scanning(img, 2)
  expected = centered(expected)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  # 3-D array
  img = zeros(Float64, (7,7,7))
  img[3:5,2:6,1:3] .= 2
  img[1:4,5:7,2:6] .= 3
  img[3:5,3:5,3:5] .= 8

  seg = fast_scanning(img, 1.5)
  expected = fill(5, (7,7,7))
  expected[3:5,2:6,1:3] .= 3
  expected[1:4,5:7,2:6] .= 3
  expected[3:5,3:5,3:5] .= 4
  expected_labels = [3,4,5]
  expected_means = Dict(3=>222/84,4=>8.0,5=>0.0)
  expected_count = Dict(3=>84,4=>27,5=>232)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  # custom diff_fn
  img = zeros(RGB{N0f8}, (3,3))
  img[2,:] .= RGB{N0f8}(0.2, 0.4, 0.4)
  img[3,:] .= RGB{N0f8}(0.2, 1.0, 1.0)

  seg = fast_scanning(img, 0.22)
  expected = ones(Int, (3,3))
  expected[2,:] .= 2
  expected[3,:] .= 3
  expected_labels = [1,2,3]
  expected_means = Dict(1=>RGB{Float64}(0.0,0.0,0.0),2=>RGB{Float64}(0.2,0.4,0.4),3=>RGB{Float64}(0.2,1.0,1.0))
  expected_count = Dict(1=>3,2=>3,3=>3)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  seg = fast_scanning(img, 0.22, (i,j)->(abs(i.r-j.r)))
  expected = ones(Int, (3,3))
  expected_labels = [1]
  expected_means = Dict(1=>RGB{Float64}(1.2/9,4.2/9,4.2/9))
  expected_count = Dict(1=>9)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  # Adaptive Thresholding
  img = [ 0.0 0.1 0.7 0.9;
          0.0 0.1 0.7 0.9;
          0.2 0.3 0.7 0.9;
          0.2 0.3 0.7 0.9; ]

  seg = fast_scanning(img, [0.2 0.1; 0.1 0.3])
  expected =  [ 1 1 4 4;
                1 1 4 4;
                2 2 4 4;
                2 2 4 4; ]
  expected_labels = [1,2,4];
  expected_means = Dict(1=>0.05,2=>0.25,4=>0.8)
  expected_count = Dict(1=>4,2=>4,4=>8)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected

  seg = fast_scanning(img, (2,2))
  expected =  [ 1 3 5 6;
                1 3 5 6;
                2 4 5 6;
                2 4 5 6; ]
  expected_labels = [1,2,3,4,5,6]
  expected_means = Dict(1=>0.0,2=>0.2,3=>0.1,4=>0.3,5=>0.7,6=>0.9)
  expected_count = Dict(1=>2,2=>2,3=>2,4=>2,5=>4,6=>4)

  @test all(label->(label in expected_labels), seg.segment_labels)
  @test all(label->(label in seg.segment_labels), expected_labels)
  @test expected_count == seg.segment_pixel_count
  @test all(label->(expected_means[label] ≈ seg.segment_means[label]), seg.segment_labels)
  @test seg.image_indexmap == expected
end


text:
# Traits used for statistical models

"""
    implicit_intercept(T::Type)
    implicit_intercept(x::T) = implicit_intercept(T)

Return `true` if models of type `T` should include an implicit intercept even
if none is specified in the formula.  Is `true` by default for all
`T<:StatisticalModel`, and `false` for others.  To specify that a model type `T`
includes an intercept even if one is not specified explicitly in the formula,
overload this function for the corresponding type:
`implicit_intercept(::Type{<:T}) = true`

If a model has an implicit intercept, it can be explicitly excluded by using `0`
in the formula, which generates [`InterceptTerm{false}`](@ref InterceptTerm) with
[`apply_schema`](@ref).
"""
implicit_intercept(x::T) where {T} = implicit_intercept(T)
implicit_intercept(::Type{T}) where {T} = false
implicit_intercept(::Type{<:StatisticalModel}) = true

"""
    drop_intercept(T::Type)
    drop_intercept(x::T) = drop_intercept(T)

Define whether a given model automatically drops the intercept. Return `false` by default.
To specify that a model type `T` drops the intercept, overload this function for the
corresponding type: `drop_intercept(::Type{<:T}) = true`

Models that drop the intercept will be fitted without one: the intercept term will be
removed even if explicitly provided by the user. Categorical variables will be expanded
in the rank-reduced form (contrasts for `n` levels will only produce `n-1` columns).
"""
drop_intercept(x::T) where {T} = drop_intercept(T)
drop_intercept(::Type) = false


text:
# test DrakeVisualizer

# for more low level interaction with the DraekVisualizer/Director, please see
# www.github.com/rdeits/DrakeVisualizer.jl

using Caesar, CoordinateTransformations, Rotations
using RoME, Distributions
using TransformUtils
using CloudGraphs
using IncrementalInference
# for illustration only
using GeometryTypes
using DrakeVisualizer
import DrakeVisualizer: Triad
import Caesar: prepcolordepthcloud!

 # @load "usercfg.jld"
include(joinpath(dirname(@__FILE__),"..","database","blandauthremote.jl"))
user_config = addrdict
user_config["session"] = "SESSHAUVDEV1"
backend_config, user_config = standardcloudgraphsetup(addrdict=user_config)

vis = startdefaultvisualization(draworigin=true)
# Launch the viewer application if it isn't running already:
# DrakeVisualizer.any_open_windows() || DrakeVisualizer.new_window();

# Mongo keys with BSONpointclouds, take from SESSHAUVDEV5
mkpc = Vector{AbstractString}(5)
mkpc[1] = "58f268f7b2ca1516ff509b82"
mkpc[2] = "58f268feb2ca1516ff509b84"
mkpc[3] = "58f26900b2ca1516ff509b86"
mkpc[4] = "58f26902b2ca1516ff509b88"
mkpc[5] = "58f26904b2ca1516ff509b8a"

wTx = Vector{AffineMap}(5)
wTx[1] = Translation(8.7481, 7.51, 4.5662) ∘ LinearMap(Quat(0.9939, -0.0280019, -0.00824962, -0.106355))
wTx[2] = Translation(9.43, 11.0879, 4.9) ∘ LinearMap(Quat(0.992387, 0.0575732, -0.00698004, -0.108645))
wTx[3] = Translation(9.94, 15.1484, 4.8792) ∘ LinearMap(Quat(0.997597, 0.0563353, -0.0115003, -0.0386678))
wTx[4] = Translation(9.93104, 14.5831, 5.8) ∘ LinearMap(Quat(0.997626, 0.0336346, -0.011705, -0.0589371))
wTx[5] = Translation(9.12, 10.1615, 5.84927) ∘ LinearMap(Quat(0.998351, 0.0362801, -0.00973871, -0.0434122))

for i in 1:5
  # set pose
  sym = Symbol("x$(i)")
  setgeometry!(vis[:debughauv][sym], Triad())
  settransform!(vis[:debughauv][sym],   wTx[i])
  # and a point cloud
  data = read_MongoData(backend_config, mkpc[i])
  pcarr = getPointCloudFromBSON(data)
  pc = prepcolordepthcloud!(i, pcarr)
  setgeometry!(vis[:debughauv][sym][:pointcloud], pc)
end

# >>> view.camera().SetViewUp([0,0,-1])
# >>> view.render()

# fetch and study the graph as built by server.jl

# setBackendWorkingSet!(fg.cg.neo4j.connection, user_config["session"])
fg = initfg(cloudgraph=backend_config, sessionname=user_config["session"])
fullLocalGraphCopy!(fg, reqbackendset=true)
writeGraphPdf(fg)
run(`evince fg.pdf`)

tree = wipeBuildNewTree!(fg,drawpdf=true)
@async run(`evince bt.pdf`)

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl", 144)

ret = inferOverTreeR!(fg, tree, N=100, dbg=true)
visualizeallposes!(vis, fg)

spyCliqMat(tree, :x2)
plotUpMsgsAtCliq(tree, :x1, :x2, marg=[6;3])

plotPriorsAtCliq(tree, :x1, :x5, dims=[1;2])
plotPriorsAtCliq(tree, :x1, :x5, dims=[6;3])
plotPriorsAtCliq(tree, :x1, :x5, dims=[4;5])

drawAllBinaryFactorEdges!(vis, backend_config, user_config["session"])

# drawLineBetween!(vis, fg,:x1,:x2, api=localapi , color=RGBA(0,1.0,1.0,0.5), scale=0.045 )
# drawAllOdometryEdges!(vis, fg, api=localapi)
# drawAllBinaryFactorEdges!(vis, fg)

X1 = getBelief(fg, :x1)
plotKDE(X1, dims=[1;2])
X2 = getBelief(fg, :x2)
plotKDE(X2, dims=[1;2])

pp, = localProduct(fg, :x2)
plotKDE(pp, dims=[1;2])

## Recalculate XYH
# change toolbox
wTx1 = convert(SE3, wTx[1])
wTx2 = convert(SE3, wTx[2])
wTx3 = convert(SE3, wTx[3])
wTx4 = convert(SE3, wTx[4])
wTx5 = convert(SE3, wTx[5])
# get rotation to world frame for local level
wEx1 = convert(Euler, wTx1.R);  wEx1.Y = 0.0;
wRlx1 = SE3(zeros(3), wEx1)
wEx2 = convert(Euler, wTx2.R);  wEx2.Y = 0.0;
wRlx2 = SE3(zeros(3), wEx2)
wEx3 = convert(Euler, wTx3.R);  wEx3.Y = 0.0;
wRlx3 = SE3(zeros(3), wEx3)
wEx4 = convert(Euler, wTx4.R);  wEx4.Y = 0.0;
wRlx4 = SE3(zeros(3), wEx4)
wEx5 = convert(Euler, wTx5.R);  wEx5.Y = 0.0;
wRlx5 = SE3(zeros(3), wEx5)
# Odometries
wRlx1Tx2 = wRlx1 * (wTx1\wTx2)
wRlx2Tx3 = wRlx2 * (wTx2\wTx3)
wRlx3Tx4 = wRlx3 * (wTx3\wTx4)
wRlx4Tx5 = wRlx4 * (wTx4\wTx5)

vEx1_2 = veeEuler(wRlx1Tx2)
vEx2_3 = veeEuler(wRlx2Tx3)
vEx3_4 = veeEuler(wRlx3Tx4)
vEx4_5 = veeEuler(wRlx4Tx5)

XYH1_2 = [vEx1_2[1], vEx1_2[2], vEx1_2[6]]
XYH2_3 = [vEx2_3[1], vEx2_3[2], vEx2_3[6]]
XYH3_4 = [vEx3_4[1], vEx3_4[2], vEx3_4[6]]
XYH4_5 = [vEx4_5[1], vEx4_5[2], vEx4_5[6]]

odo_xyh_12 = getfnctype(getVert(fg, fg.fIDs[:x1x2]))
odo_xyh_23 = getfnctype(getVert(fg, fg.fIDs[:x2x3]))
odo_xyh_34 = getfnctype(getVert(fg, fg.fIDs[:x3x4]))
odo_xyh_45 = getfnctype(getVert(fg, fg.fIDs[:x4x5]))

zpr2 = veeEuler(wTx2)[[3;4;5]]
zpr3 = veeEuler(wTx3)[[3;4;5]]
zpr4 = veeEuler(wTx4)[[3;4;5]]
zpr5 = veeEuler(wTx5)[[3;4;5]]

# compare differences, these should all be zero
@show odo_xyh_12.xyy.μ - XYH1_2
@show odo_xyh_23.xyy.μ - XYH2_3
@show odo_xyh_34.xyy.μ - XYH3_4
@show odo_xyh_45.xyy.μ - XYH4_5

# build a new factor graph from this info
N = 100
fgnew = initfg(cloudgraph=backend_config, sessionname="SESSHAUVCODED")

addVariable!(fgnew,:x1,N=N,tags=[:POSE;], dims=6)
addFactor!(fgnew,[:x1], PriorPose3(MvNormal(veeEuler(wTx1), 0.0001*eye(6)) ) )
initVariable!(fgnew, :x1)
plotKDE(fgnew, :x1, dims=[1;2], title="$(veeEuler(wTx1))")

addVariable!(fgnew,:x2,N=N,tags=[:POSE;], dims=6)
a,b = MvNormal(zpr2[2:3],0.001*eye(2)), Distributions.Normal(zpr2[1], 0.01)
addFactor!(fgnew,[:x2], PartialPriorRollPitchZ( a, b ) )
addFactor!(fgnew,[:x1,:x2], PartialPose3XYYaw(MvNormal(XYH1_2,0.001*eye(3)))  )
initVariable!(fgnew, :x2)
plotKDE(fgnew, :x2, dims=[6;3], title="$(veeEuler(wTx2))")

visualizeallposes!(vis, fgnew)

pp, parr, par, lb, infdim = localProduct(fgnew, :x2)
plotKDE(pp, dims=[6;3], title="$(veeEuler(wTx2))")
# plotLocalProduct(fgnew, :x2, dims=[6;3], api=localapi)

# addVariable!(fgnew,:x2,N=N,tags=[:POSE;], dims=6)
# a,b = MvNormal(zpr2[2:3],0.001*eye(2)), Distributions.Normal(zpr2[1], 0.1)
# addFactor!(fgnew,[:x2], PartialPriorRollPitchZ( a, b ) )
# addFactor!(fgnew,[:x1,:x2], deepcopy(odo_xyh_12)  )
# initVariable!(fgnew, :x2)
# plotKDE(fgnew, :x2, dims=[1;2], title="$(veeEuler(wTx2))")

# addVariable!(fgnew,:x2,N=N,tags=[:POSE;], dims=6)
# a,b = MvNormal(zpr2[2:3],0.001*eye(2)), Distributions.Normal(zpr2[1], 0.1)
# addFactor!(fgnew,[:x2], PartialPriorRollPitchZ( a, b ) )
# odo = MvNormal( veeEuler(wTx1\wTx2),0.001*eye(6))
# addFactor!(fgnew,[:x1,:x2], Pose3Pose3(odo)  )
# initVariable!(fgnew, :x2)
# plotKDE(fgnew, :x2, dims=[1;2], title="$(veeEuler(wTx2))")

addVariable!(fgnew,:x3,N=N,tags=[:POSE;], dims=6)
a,b = MvNormal(zpr3[2:3],0.001*eye(2)), Distributions.Normal(zpr3[1], 0.1)
addFactor!(fgnew,[:x3], PartialPriorRollPitchZ( a, b ) )
addFactor!(fgnew,[:x2,:x3], PartialPose3XYYaw(MvNormal(XYH2_3,0.001*eye(3)))  )
initVariable!(fgnew, :x3)
plotKDE(fgnew, :x3, dims=[1;2], title="$(veeEuler(wTx3))")

visualizeallposes!(vis, fgnew)

addVariable!(fgnew,:x4,N=N,tags=[:POSE;], dims=6)
a,b = MvNormal(zpr4[2:3],0.001*eye(2)), Distributions.Normal(zpr4[1], 0.1)
addFactor!(fgnew,[:x4], PartialPriorRollPitchZ( a, b ) )
addFactor!(fgnew,[:x3,:x4], PartialPose3XYYaw(MvNormal(XYH3_4,0.001*eye(3)))  )
initVariable!(fgnew, :x4)
plotKDE(fgnew, :x4, dims=[1;2], title="$(veeEuler(wTx4))")

addVariable!(fgnew,:x5,N=N,tags=[:POSE;], dims=6)
a,b = MvNormal(zpr5[2:3],0.001*eye(2)), Distributions.Normal(zpr5[1], 0.1)
addFactor!(fgnew,[:x5], PartialPriorRollPitchZ( a, b ) )
addFactor!(fgnew,[:x4,:x5], PartialPose3XYYaw(MvNormal(XYH4_5,0.001*eye(3)))  )
initVariable!(fgnew, :x5)
plotKDE(fgnew, :x5, dims=[1;2], title="$(veeEuler(wTx5))")

visualizeallposes!(vis, fgnew)
solve(fgnew)

treenew = wipeBuildNewTree!(fgnew, drawpdf=true)
inferOverTreeR!(fgnew, treenew, N=100, dbg=true)

plotUpMsgsAtCliq(treenew, :x1, :x2, marg=[6;3])

plotPriorsAtCliq(treenew, :x1, :x5, dims=[1;2])
plotPriorsAtCliq(treenew, :x1, :x5, dims=[6;3])
plotPriorsAtCliq(treenew, :x1, :x5, dims=[4;5])

#


text:
# The code to generate Fig 8
using PyPlot
using SurNoR_2020
using Statistics
using MAT
using HypothesisTests

PyPlot.svg(true)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["svg.fonttype"] = "none"
rcParams["pdf.fonttype"] = 42

# WARNING: put n_set = 1 for the left panel, put n_set = 2 for the middle,
#          and put n_set = 3 for the right panel
n_set = 3
# ------------------------------------------------------------------------------
# Load simulated data
# ------------------------------------------------------------------------------
Sim_Data = load("data/Simulated_data/Simulated_data.jld")
Inputs = Sim_Data["Inputs"]

Sub_Num = 12
N_start = (n_set-1)*Sub_Num
BehavData = Inputs[N_start .+ (1:Sub_Num),:,:]

# ------------------------------------------------------------------------------
# Training Algorithms
# ------------------------------------------------------------------------------
Param_Path = string("data/Fitted_parameters_simulated/CV_3Fold/nset_", string(n_set), "/")
Fold_Set = [[9,10,11,12],[5,6,7,8],[1,2,3,4]]
for Part_ind = 1:4
    Agents_MFNS    = SurNoR_2020.Func_MFNS_CV_train(BehavData=BehavData,
                                                    Param_Path=Param_Path,
                                                    Fold_Set=Fold_Set,
                                                    Part_ind=Part_ind)
    Agents_HybSExp = SurNoR_2020.Func_HybSExp_CV_train(BehavData=BehavData,
                                                       Param_Path=Param_Path,
                                                       Fold_Set=Fold_Set,
                                                       Part_ind=Part_ind)
    Agents_HybN    = SurNoR_2020.Func_HybN_CV_train(BehavData=BehavData,
                                                    Param_Path=Param_Path,
                                                    Fold_Set=Fold_Set,
                                                    Part_ind=Part_ind)
    Agents_SurNoR  = SurNoR_2020.Func_SurNoR_CV_train(BehavData=BehavData,
                                                      Param_Path=Param_Path,
                                                      Fold_Set=Fold_Set,
                                                      Part_ind=Part_ind)

    if Part_ind == 1
        global Agents_Tot_Alg = cat(Agents_MFNS,Agents_HybSExp, Agents_HybN, Agents_SurNoR, dims=4)

    else
        Agents_Tot_Alg_temp = cat(Agents_MFNS,Agents_HybSExp, Agents_HybN, Agents_SurNoR, dims=4)

        Agents_Tot_Alg = cat(Agents_Tot_Alg, Agents_Tot_Alg_temp, dims=5)
    end
end
SurNoR_ind = 4   # index corresponding to SurNoR

# Extracting Log-likelihoods
Alg_Num = size(Agents_Tot_Alg)[4]
Part_Num = size(Agents_Tot_Alg)[5]
LogL_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))
LogL_Unif_Matrix = zeros((Sub_Num,5,2,Alg_Num,Part_Num))

for Sub = 1:Sub_Num
    for Epi = 1:5
        for Block = 1:2
            for Alg = 1:Alg_Num
                for Part_ind = 1:Part_Num
                    LogL_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].LogL
                    LogL_Unif_Matrix[Sub,Epi,Block,Alg,Part_ind] =
                            Agents_Tot_Alg[Sub,Block,Epi,Alg,Part_ind].LogL_Unif
                end
            end
        end
    end
end

LogL_Matrix_Epi1 = LogL_Matrix[:,1,1,:,:]
LogL_Unif_Matrix_Epi1 = LogL_Unif_Matrix[:,1,1,:,:]

LogL_Matrix_Epi1B2 = LogL_Matrix[:,1,2,:,:]
LogL_Unif_Matrix_Epi1B2 = LogL_Unif_Matrix[:,1,2,:,:]

LogL_Sum = sum(LogL_Matrix,dims=1)
LogL_Unif_Sum = sum(LogL_Unif_Matrix,dims=1)

LogL_Sum_Epi1 = sum(LogL_Matrix_Epi1,dims=1)
LogL_Unif_Sum_Epi1 = sum(LogL_Unif_Matrix_Epi1,dims=1)

LogL_Sum_Epi1B2 = sum(LogL_Matrix_Epi1B2,dims=1)
LogL_Unif_Sum_Epi1B2 = sum(LogL_Unif_Matrix_Epi1B2,dims=1)

dLogL_Sum_Block = sum(LogL_Sum - LogL_Unif_Sum,dims=2)
dLogL_Sum_Block_Epi1 = LogL_Sum_Epi1 - LogL_Unif_Sum_Epi1
dLogL_Sum_Block_Epi1B2 = LogL_Sum_Epi1B2 - LogL_Unif_Sum_Epi1B2

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 8A
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
x = 1:Alg_Num
Colors = ["#6E0D0D","#C785EA","#B159E0","#8723BE"]
Legend = ["MF+S+N","Hyb+S+U", "Hyb+N", "SurNoR"]

# All data
y = dLogL_Sum_Block[1,1,1,:,:] .+ dLogL_Sum_Block[1,1,2,:,:]
dy = std(y,dims=2)/sqrt(Part_Num)
y = mean(y,dims=2)
fig = figure(figsize=(6,6.5)); ax = gca()
for i = 1:Alg_Num
   ax.bar(x[i],y[i],color = Colors[i])
   ax.errorbar(x[i],y[i],yerr=dy[i],color="k",linewidth=1,drawstyle="steps",capsize=3)
end
title("Delta Log Evidence - w.r.t. RC - all epi. for both blocks")
ax.set_xticks(1:Alg_Num)
if n_set == 1
    ax.set_ylim([1400,1900])
elseif n_set == 2
    ax.set_ylim([2000,2300])
elseif n_set == 3
    ax.set_ylim([2300,2700])
end
ax.set_xticklabels(Legend,rotation=90)
ax.set_ylabel("delta Log Evidence");
fig.subplots_adjust(bottom = 0.2, left = 0.15)

# ------------------------------------------------------------------------------
# Save LogL as .MAT for Bayesian Model Selection (KE Stephan, et al. 2009)
# ------------------------------------------------------------------------------
LogL_MAT = mean(LogL_Matrix,dims=5)
LogL_MAT = cat(LogL_MAT[:,:,:,1:Alg_Num,1],
               LogL_Unif_Matrix[:,:,:,1,1], dims=4)

file = matopen(string("data/Log_evidences_simulated/nset_", string(n_set), "/LogL_Matrix.mat"), "w")
write(file, "LogL_Matrix", LogL_MAT)
close(file)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Fig 8B: Expected posterior
# WARNING: Before running this part, in principle, one needs to run the MATLAB
#          code "./data/Log_evidences_simulated/LogL_BMS.m" to use package SPM12 and
#          produce the results and save the file "./data/Log_evidences_simulated/nset_i/BMS_Data.mat"
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
BMS = matread(string("data/Log_evidences_simulated/nset_", string(n_set), "/BMS_Data.mat"))
Colors = ["#6E0D0D","#C785EA","#B159E0","#8723BE", "#474848"]

Legend = ["MF+S+N","Hyb+S+U", "Hyb+N", "SurNoR", "RC"]

fig = figure(figsize=(7,6)); ax = gca()
for i = 1:5
   ax.bar(i,BMS["BMS"]["exp_r"][i], color=Colors[i])
   if (i==SurNoR_ind)
        ax.text(i,BMS["BMS"]["exp_r"][i]+0.01,
                string("ϕ = ",string(round(BMS["BMS"]["pxp"][i],digits=5))),
                fontsize=12, horizontalalignment="center", rotation=0)
   end
end
ax.plot([0,6],[1,1] ./ 5, linestyle="dashed",linewidth=1,color="k")
title("Posterior Probabilities for Different Models - Uniform Prior")
ax.set_xticks(1:5)
ax.set_xticklabels(Legend,rotation=90)
ax.set_yticks(0:0.1:1)
ax.set_yticklabels(0:0.1:1)
ax.set_ylabel("Model Posterior Probability");
ax.set_ylim([0,1])
ax.set_xlim([0.25,5.75])
fig.subplots_adjust(bottom = 0.2, left = 0.15)


text:
module Nexus

import Base:
show, ==, empty!, getindex, setindex!, length, iterate, rand, eltype,
firstindex, lastindex, isempty, keys
using Dictionaries
using StatsBase
using DataStructures

include("graph.jl")
include("arc.jl")
include("node.jl")
include("path.jl")
include("show.jl")

export Graph, nn, na, fadj, badj,
add_arc!, rem_arc!, add_arcs_from!, rem_arcs_from!,
add_node!, rem_node!, add_nodes_from!, rem_nodes_from!, # graph.jl
Arc, source, target, key, has_arc, ArcIter, arcs, has_attrib, # arc. jl
Path, nodelist, keylist,
depth_first_search, dijkstra_shortest_path, # path.jl
nodes, in_neighbor, out_neighbor, has_node # node.jl

end


text:
module PowerFlowNetworks

using Graphs, MetaGraphs
using SQLite, DataFrames
using Dates
using Statistics: mean, median, var
using Serialization
using DataStructures
using Random
using DelimitedFiles
using UUIDs
using SparseArrays
using LinearAlgebra
using Combinatorics
import JSON
using MPI

# Solve module
using Logging
using CSV
using Printf
using JuMP
using Mosek
using MosekTools
import Base.Iterators: flatten

mutable struct PowerFlowNetwork
    name::String
    bus::DataFrame
    gen::DataFrame
    branch::DataFrame
    gencost::DataFrame
    baseMVA::Float64
end

include("core.jl")
include("utils/graphs/lbfs.jl")
include("utils/graphs/ischordal.jl")
include("utils/graphs/build_graph.jl")
include("utils/clique_cliquetree.jl")
include("io/read.jl")
include("io/write.jl")
include("read_features.jl")
include("graphs/graphs.jl")
include("graphs/operations.jl")
include("db/setup_db.jl")
include("db/infos.jl")
include("db/inserts.jl")
include("utils/graphs/operators.jl")
include("decompositions/chordal_extension/chordal_extension.jl")
include("decompositions/merge/merge.jl")
include("decompositions/combine/combine.jl")
include("db/operations.jl")
include("solve/solve.jl")
include("mpi.jl")

PowerFlowNetwork(path::AbstractString, format::AbstractString) = read_network(path, format)
SimpleGraph(network::PowerFlowNetwork) = to_simple_graph(network)

export nbus, nbranch, ngen, is_disjoint, has_bus, has_branch, has_gen,
       has_continuous_index, normalize_index, merge_duplicate_branch!
export ischordal
export read_cliquetree, read_clique, get_nv, get_ne, get_nb_lc
export nbranch_unique, ngen_unique
export PowerFlowNetwork
export write_pfn
export Graph
export add_edges_distance!, add_edges_random!, add_edges!
export setup_db
export has_opf_tables, count_missing_columns, table_count, nb_instances, nb_decompositions, get_table_ids
export load_instance_in_db!
export chordal_extension
export merge_dec
export combine_graph
export save_features_instances!, save_basic_features_instances!, save_single_features_instances!,
       serialize_instances!, generate_decompositions!, merge_decompositions!, combine_decompositions!,
       check_sanity, delete_duplicates!, export_matpowerm_instances!, load_matctr_instances!, solve_decompositions!
export solve_sdp
export execute_process_mpi

end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libidn2"
version = v"2.3.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/libidn/libidn2-$(version).tar.gz", "76940cd4e778e8093579a9d195b25fff5e936e9dc6242068528b437a76764f91")
]

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

 ./configure \
 --prefix=${prefix} \
 --build=${MACHTYPE} \
 --host=${target} \
 --enable-static=no \
 --enable-shared=yes \
 --disable-doc

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental = true)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531"))
    Dependency(PackageSpec(name="libunistring_jll", uuid="6db05002-db9d-53dd-a359-17d4854bdf22"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("Data.jl")

function load(instance::Symbol)
    instance = replace(string(instance), "_" => "-")
    file_name = joinpath(data_path, instance * ".dat")
    return load(file_name)
end

function load(file_name::AbstractString)
    if !isfile(file_name)
        println("Error loading file ", file_name)
        return nothing
    end

    values = parse.(Int64, split(read(file_name, String)))
    name = splitext(basename(file_name))[1]

    n = values[1]
    vertices = Vector{Vertex}()
    for v in 1:n
        push!(vertices, Vertex(v, Vector{Edge}(), Vector{Edge}()))
    end

    m = values[2]
    edges = Vector{Edge}()
    requireds = Vector{Edge}()

    i = 3
    for e = 1:m
        from = vertices[values[i] + 1]
        to = vertices[values[i + 1] + 1]
        cost = values[i + 2]
        demand = values[i + 3]
        i += 4

        edge = Edge(e, from, to, cost, demand)

        push!(edges, edge)
        push!(from.edges, edge)
        push!(to.edges, edge)

        if isRequired(edge)
            push!(requireds, edge)
            push!(from.requireds, edge)
            push!(to.requireds, edge)
        end
    end

    vehicles = values[i]
    capacity = values[i + 1]

    lb = values[i + 2]
    ub = values[i + 3]

    data = Data(name, capacity, vehicles, vertices, edges, requireds, lb, ub)

    return data
end


text:
using RayTracingWeekend
using BenchmarkTools
#using StaticArrays
#using Images

using Images, LinearAlgebra, Random, RandomNumbers.Xorshifts, StaticArrays

# const Vec3{T<:AbstractFloat} = SVector{3, T}

# # Adapted from @Christ_Foster's: https://discourse.julialang.org/t/ray-tracing-in-a-week-end-julia-vs-simd-optimized-c/72958/40
# @inline @inbounds function Base.getproperty(v::SVector{3}, name::Symbol)
#     name === :x || name === :r ? getfield(v, :data)[1] :
#     name === :y || name === :g ? getfield(v, :data)[2] :
#     name === :z || name === :b ? getfield(v, :data)[3] : getfield(v, name)
# end

# @inline @inbounds function Base.getproperty(v::SVector{2}, name::Symbol)
#     name === :x || name === :r ? getfield(v, :data)[1] :
#     name === :y || name === :g ? getfield(v, :data)[2] : getfield(v, name)
# end

# @inline squared_length(v::Vec3) = v ⋅ v
# @inline near_zero(v::Vec3) = squared_length(v) < 1e-5
# @inline rgb(v::Vec3) = RGB(v...)
# @inline rgb_gamma2(v::Vec3) = RGB(sqrt.(v)...)

# struct Ray{T}
# 	origin::Vec3{T} # Point
# 	dir::Vec3{T} # Vec3 # direction (unit vector)
# end

# @inline point(r::Ray{T}, t::T) where T <: AbstractFloat = r.origin .+ t .* r.dir # equivalent to C++'s ray.at()

# @inline function skycolor(ray::Ray{T}) where T
# 	white = SA[1.0, 1.0, 1.0]
# 	skyblue = SA[0.5, 0.7, 1.0]
# 	t = T(0.5)*(ray.dir.y + one(T))
#     (one(T)-t)*white + t*skyblue
# end

# claforte: if I comment out this TRNG re-definition, the multithreading crash happens!

# Instantiate 1 RNG (Random Number Generator) per thread, for performance
# Fix the random seeds, to make it easier to benchmark changes.
#const TRNG = [Xoroshiro128Plus(i) for i = 1:128]    # claforte CRASH: When I commented-out this line, the crash occured.

# reseed!() = for (i,rng) in enumerate(TRNG) Random.seed!(rng, i) end # reset the seed
# reseed!()

# @inline function trand() # thread-local rand()
#     @inbounds rng = TRNG[Threads.threadid()]
#     rand(rng)
# end

# @inline function trand(::Type{T}) where T # thread-local rand()
#     @inbounds rng = TRNG[Threads.threadid()]
#     rand(rng, T)
# end

# @inline function random_vec3_in_sphere(::Type{T}) where T # equiv to random_in_unit_sphere()
# 	while true
# 		p = random_vec3(T(-1), T(1))
# 		if p⋅p <= 1
# 			return p
# 		end
# 	end
# end

# @inline random_between(min::T=0, max::T=1) where T = trand(T)*(max-min) + min # equiv to random_double()
# @inline random_vec3(min::T, max::T) where T = @inbounds @SVector[random_between(min, max) for i ∈ 1:3]
# @inline random_vec2(min::T, max::T) where T = @inbounds @SVector[random_between(min, max) for i ∈ 1:2]

# "Random unit vector. Equivalent to C++'s `unit_vector(random_in_unit_sphere())`"
# @inline random_vec3_on_sphere(::Type{T}) where T = normalize(random_vec3_in_sphere(T))

# @inline function random_vec2_in_disk(::Type{T}) where T # equiv to random_in_unit_disk()
# 	while true
# 		p = random_vec2(T(-1), T(1))
# 		if p⋅p <= 1
# 			return p
# 		end
# 	end
# end

# "An object that can be hit by Ray"
# abstract type Hittable end

# """Materials tell us how rays interact with a surface"""
# abstract type Material{T <: AbstractFloat} end

# "Record a hit between a ray and an object's surface"
# struct HitRecord{T <: AbstractFloat}
# 	t::T # distance from the ray's origin to the intersection with a surface.

# 	p::Vec3{T} # point of the intersection between an object's surface and a ray
# 	n⃗::Vec3{T} # surface's outward normal vector, points towards outside of object?

# 	# If true, our ray hit from outside to the front of the surface.
# 	# If false, the ray hit from within.
# 	front_face::Bool

# 	mat::Material{T}

# 	@inline HitRecord(t::T,p,n⃗,front_face,mat) where T = new{T}(t,p,n⃗,front_face,mat)
# end

# struct Sphere{T <: AbstractFloat} <: Hittable
# 	center::Vec3{T}
# 	radius::T
# 	mat::Material{T}
# end

# """Equivalent to `hit_record.set_face_normal()`"""
# @inline @fastmath function ray_to_HitRecord(t::T, p, outward_n⃗, r_dir::Vec3{T}, mat::Material{T})::Union{HitRecord,Nothing} where T
# 	front_face = r_dir ⋅ outward_n⃗ < 0
# 	n⃗ = front_face ? outward_n⃗ : -outward_n⃗
# 	HitRecord(t,p,n⃗,front_face,mat)
# end

# struct Scatter{T<: AbstractFloat}
# 	r::Ray{T}
# 	attenuation::Vec3{T}

# 	# claforte: TODO: rename to "absorbed?", i.e. not reflected/refracted?
# 	reflected::Bool # whether the scattered ray was reflected, or fully absorbed
# 	@inline Scatter(r::Ray{T},a::Vec3{T},reflected=true) where T = new{T}(r,a,reflected)
# end

# struct Lambertian{T} <: Material{T}
# 	albedo::Vec3{T}
# end

# """Compute reflection vector for v (pointing to surface) and normal n⃗.

# 	See [diagram](https://raytracing.github.io/books/RayTracingInOneWeekend.html#metal/mirroredlightreflection)"""
# @inline @fastmath reflect(v::Vec3{T}, n⃗::Vec3{T}) where T = v - (2v⋅n⃗)*n⃗

# """Create a scattered ray emitted by `mat` from incident Ray `r`.

# 	Args:
# 		rec: the HitRecord of the surface from which to scatter the ray.

# 	Return `nothing`` if it's fully absorbed. """
# @inline @fastmath function scatter(mat::Lambertian{T}, r::Ray{T}, rec::HitRecord{T})::Scatter{T} where T
# 	scatter_dir = rec.n⃗ + random_vec3_on_sphere(T)
# 	if near_zero(scatter_dir) # Catch degenerate scatter direction
# 		scatter_dir = rec.n⃗
# 	else
# 		scatter_dir = normalize(scatter_dir)
# 	end
# 	scattered_r = Ray{T}(rec.p, scatter_dir)
# 	attenuation = mat.albedo
# 	return Scatter(scattered_r, attenuation)
# end

# @inline @fastmath function hit(s::Sphere{T}, r::Ray{T}, tmin::T, tmax::T)::Union{HitRecord,Nothing} where T
#     oc = r.origin - s.center
#     #a = r.dir ⋅ r.dir # unnecessary since `r.dir` is normalized
# 	a = 1
# 	half_b = oc ⋅ r.dir
#     c = oc⋅oc - s.radius^2
#     discriminant = half_b^2 - a*c
# 	if discriminant < 0 return nothing end # no hit!
# 	sqrtd = √discriminant

# 	# Find the nearest root that lies in the acceptable range
# 	root = (-half_b - sqrtd) / a
# 	if root < tmin || tmax < root
# 		root = (-half_b + sqrtd) / a
# 		if root < tmin || tmax < root
# 			return nothing # no hit!
# 		end
# 	end

# 	t = root
# 	p = point(r, t)
# 	n⃗ = (p - s.center) / s.radius
# 	return ray_to_HitRecord(t, p, n⃗, r.dir, s.mat)
# end

# const HittableList = Vector{Hittable}

# """Find closest hit between `Ray r` and a list of Hittable objects `h`, within distance `tmin` < `tmax`"""
# @inline function hit(hittables::HittableList, r::Ray{T}, tmin::T, tmax::T)::Union{HitRecord,Nothing} where T
#     closest = tmax # closest t so far
#     best_rec::Union{HitRecord,Nothing} = nothing # by default, no hit
#     @inbounds for i in eachindex(hittables) # @paulmelis reported gave him a 4X speedup?!
# 		h = hittables[i]
#         rec = hit(h, r, tmin, closest)
#         if rec !== nothing
#             best_rec = rec
#             closest = best_rec.t # i.e. ignore any further hit > this one's.
#         end
#     end
#     best_rec
# end

#@inline color_vec3_in_rgb(v::Vec3{T}) where T = 0.5normalize(v) + SA{T}[0.5,0.5,0.5]

# struct Metal{T} <: Material{T}
# 	albedo::Vec3{T}
# 	fuzz::T # how big the sphere used to generate fuzzy reflection rays. 0=none
# 	@inline Metal(a::Vec3{T}, f::T=0.0) where T = new{T}(a,f)
# end
# @inline @fastmath function scatter(mat::Metal{T}, r_in::Ray{T}, rec::HitRecord)::Scatter{T} where T
# 	reflected = normalize(RaytracingWeekend.reflect(r_in.dir, rec.n⃗) + mat.fuzz*random_vec3_on_sphere(T)) # claforte: fully-qualified, to reuse that module's definition
# 	Scatter(Ray(rec.p, reflected), mat.albedo)
# end

# "Scene with 2 Lambertian spheres"
# function scene_2_spheres(; elem_type::Type{T}) where T
# 	spheres = Sphere[]

# 	# small center sphere
# 	push!(spheres, Sphere((SA{T}[0,0,-1]), T(0.5), Lambertian(SA{T}[0.7,0.3,0.3])))

# 	# ground sphere
# 	push!(spheres, Sphere((SA{T}[0,-100.5,-1]), T(100), Lambertian(SA{T}[0.8,0.8,0.0])))
# 	HittableList(spheres)
# end

# """Scene with 2 Lambertian, 2 Metal spheres.

# 	See https://raytracing.github.io/images/img-1.11-metal-shiny.png"""
# function scene_4_spheres(; elem_type::Type{T}) where T
# 	scene = scene_2_spheres(; elem_type=elem_type)

# 	# left and right Metal spheres
# 	push!(scene, Sphere((SA{T}[-1,0,-1]), T(0.5), Metal((SA{T}[0.8,0.8,0.8]), T(0.3))))
# 	push!(scene, Sphere((SA{T}[ 1,0,-1]), T(0.5), Metal((SA{T}[0.8,0.6,0.2]), T(0.8))))
# 	return scene
# end

# struct Camera{T <: AbstractFloat}
# 	origin::Vec3{T}
# 	lower_left_corner::Vec3{T}
# 	horizontal::Vec3{T}
# 	vertical::Vec3{T}
# 	u::Vec3{T}
# 	v::Vec3{T}
# 	w::Vec3{T}
# 	lens_radius::T
# end

# """
# 	Args:
# 		vfov: vertical field-of-view in degrees
# 		aspect_ratio: horizontal/vertical ratio of pixels
#       aperture: if 0 - no depth-of-field
# """
# function default_camera(lookfrom::Vec3{T}=(SA{T}[0,0,0]),
# 						lookat::Vec3{T}=(SA{T}[0,0,-1]),
# 						vup::Vec3{T}=(SA{T}[0,1,0]),
# 						vfov::T=T(90), aspect_ratio::T=T(16/9),
# 						aperture::T=T(0), focus_dist::T=T(1)) where T
# 	viewport_height = T(2) * tand(vfov/T(2))
# 	viewport_width = aspect_ratio * viewport_height

# 	w = normalize(lookfrom - lookat)
# 	u = normalize(vup × w)
# 	v = w × u

# 	origin = lookfrom
# 	horizontal = focus_dist * viewport_width * u
# 	vertical = focus_dist * viewport_height * v
# 	lower_left_corner = origin - horizontal/T(2) - vertical/T(2) - focus_dist*w
# 	lens_radius = aperture/T(2)
# 	Camera{T}(origin, lower_left_corner, horizontal, vertical, u, v, w, lens_radius)
# end

# claforte CRASH: if I comment-out the following line, and run with
# `julia --project=. --threads=8 --inline=no repro_crash.jl`, I get a crash with an error message:
# `nested task error: BoundsError: attempt to access 1-element Vector{Xoroshiro128Plus} at index [2]`
# ... but if I run instead with:
# `julia --project=. --threads=8 repro_crash.jl`...
# ... I get the usual segmentation fault.
# My guess: anything that uses TRNG, directly or indirectly, needs to be redefined in this file instead of used from the module.
#

# @inline @fastmath function get_ray(c::Camera{T}, s::T, t::T) where T
# 	rd = SVector{2,T}(c.lens_radius * random_vec2_in_disk(T))
# 	offset = c.u * rd.x + c.v * rd.y #offset = c.u * rd.x + c.v * rd.y
#     Ray(c.origin + offset, normalize(c.lower_left_corner + s*c.horizontal +
# 									 t*c.vertical - c.origin - offset))
# end

# """Compute color for a ray, recursively

# 	Args:
# 		depth: how many more levels of recursive ray bounces can we still compute?"""
# @inline @fastmath function ray_color(r::Ray{T}, world::HittableList, depth=16) where T
#     if depth <= 0
# 		return SA{T}[0,0,0]
# 	end

# 	rec::Union{HitRecord,Nothing} = hit(world, r, T(1e-4), typemax(T))
#     if rec !== nothing
# 		# For debugging, represent vectors as RGB:
# 		# claforte TODO: adapt to latest code!
# 		# return color_vec3_in_rgb(rec.p) # show the normalized hit point
# 		# return color_vec3_in_rgb(rec.n⃗) # show the normal in RGB
# 		# return color_vec3_in_rgb(rec.p + rec.n⃗)
# 		# return color_vec3_in_rgb(random_vec3_in_sphere())
# 		#return color_vec3_in_rgb(rec.n⃗ + random_vec3_in_sphere())

#         s = scatter(rec.mat, r, rec)
# 		if s.reflected
# 			return s.attenuation .* ray_color(s.r, world, depth-1)
# 		else
# 			return SA{T}[0,0,0]
# 		end
#     else
#         skycolor(r)
#     end
# end

# """Render an image of `scene` using the specified camera, number of samples.

# 	Args:
# 		scene: a HittableList, e.g. a list of spheres
# 		n_samples: number of samples per pixel, eq. to C++ samples_per_pixel

# 	Equivalent to C++'s `main` function."""
# function render(scene::HittableList, cam::Camera{T}, image_width=400,
# 				n_samples=1) where T
# 	# Image
# 	aspect_ratio = T(16.0/9.0) # TODO: use cam.aspect_ratio for consistency
# 	image_height = convert(Int64, floor(image_width / aspect_ratio))

# 	# Render
# 	img = zeros(RGB{T}, image_height, image_width)
# 	f32_image_width = convert(Float32, image_width)
# 	f32_image_height = convert(Float32, image_height)

# 	# Reset the random seeds, so we always get the same images...
# 	# Makes comparing performance more accurate.
# 	reseed!()

# 	Threads.@threads for i in 1:image_height
# 		@inbounds for j in 1:image_width # iterate over each row (FASTER?!)
# 			accum_color = SA{T}[0,0,0]
# 			u = convert(T, j/image_width)
# 			v = convert(T, (image_height-i)/image_height) # i is Y-down, v is Y-up!

# 			for s in 1:n_samples
# 				if s == 1 # 1st sample is always centered
# 					δu = δv = T(0)
# 				else
# 					# Supersampling antialiasing.
# 					δu = trand(T) / f32_image_width
# 					δv = trand(T) / f32_image_height
# 				end
# 				ray = get_ray(cam, u+δu, v+δv)
# 				accum_color += ray_color(ray, scene)
# 			end
# 			img[i,j] = rgb_gamma2(accum_color / n_samples)
# 		end
# 	end
# 	img
# end

# """
# 	Args:
# 		refraction_ratio: incident refraction index divided by refraction index of
# 			hit surface. i.e. η/η′ in the figure above"""
# @inline @fastmath function refract(dir::Vec3{T}, n⃗::Vec3{T}, refraction_ratio::T) where T
# 	cosθ = min(-dir ⋅ n⃗, one(T))
# 	r_out_perp = refraction_ratio * (dir + cosθ*n⃗)
# 	r_out_parallel = -√(abs(one(T)-squared_length(r_out_perp))) * n⃗
# 	normalize(r_out_perp + r_out_parallel)
# end

# struct Dielectric{T} <: Material{T}
# 	ir::T # index of refraction, i.e. η.
# end

# @inline @fastmath function reflectance(cosθ, refraction_ratio)
# 	# Use Schlick's approximation for reflectance.
# 	# claforte: may be buggy? I'm getting black pixels in the Hollow Glass Sphere...
# 	r0 = (1-refraction_ratio) / (1+refraction_ratio)
# 	r0 = r0^2
# 	r0 + (1-r0)*((1-cosθ)^5)
# end

# @inline @fastmath function scatter(mat::Dielectric{T}, r_in::Ray{T}, rec::HitRecord{T}) where T
# 	attenuation = SA{T}[1,1,1]
# 	refraction_ratio = rec.front_face ? (one(T)/mat.ir) : mat.ir # i.e. ηᵢ/ηₜ
# 	cosθ = min(-r_in.dir⋅rec.n⃗, one(T))
# 	sinθ = √(one(T) - cosθ^2)
# 	cannot_refract = refraction_ratio * sinθ > one(T)
# 	if cannot_refract || reflectance(cosθ, refraction_ratio) > trand(T)
# 		dir = reflect(r_in.dir, rec.n⃗)
# 	else
# 		dir = refract(r_in.dir, rec.n⃗, refraction_ratio)
# 	end
# 	Scatter(Ray{T}(rec.p, dir), attenuation) # TODO: rename reflected -> !absorbed?
# end

# @inline function scene_diel_spheres(left_radius=0.5; elem_type::Type{T}) where T # dielectric spheres
# 	spheres = Sphere[]

# 	# small center sphere
# 	push!(spheres, Sphere((SA{T}[0,0,-1]), T(0.5), Lambertian(SA{T}[0.1,0.2,0.5])))

# 	# ground sphere (planet?)
# 	push!(spheres, Sphere((SA{T}[0,-100.5,-1]), T(100), Lambertian(SA{T}[0.8,0.8,0.0])))

# 	# # left and right spheres.
# 	# # Use a negative radius on the left sphere to create a "thin bubble"
# 	push!(spheres, Sphere((SA{T}[-1,0,-1]), T(left_radius), Dielectric(T(1.5))))
# 	push!(spheres, Sphere((SA{T}[1,0,-1]), T(0.5), Metal((SA{T}[0.8,0.6,0.2]), T(0))))
# 	HittableList(spheres)
# end

# function scene_blue_red_spheres(; elem_type::Type{T}) where T # dielectric spheres
# 	spheres = Sphere[]
# 	R = cos(pi/4)
# 	push!(spheres, Sphere((SA{T}[-R,0,-1]), R, Lambertian(SA{T}[0,0,1])))
# 	push!(spheres, Sphere((SA{T}[ R,0,-1]), R, Lambertian(SA{T}[1,0,0])))
# 	HittableList(spheres)
# end

# function scene_random_spheres(; elem_type::Type{T}) where T
# 	spheres = Sphere[]

# 	# ground
# 	push!(spheres, Sphere((SA{T}[0,-1000,-1]), T(1000),
# 						  Lambertian(SA{T}[0.5,0.5,0.5])))

# 	for a in -11:10, b in -11:10
# 		choose_mat = trand(T)
# 		center = SA[a + T(0.9)*trand(T), T(0.2), b + T(0.9)*trand(T)]

# 		# skip spheres too close?
# 		if norm(center - SA{T}[4,0.2,0]) < T(0.9) continue end

# 		if choose_mat < T(0.8)
# 			# diffuse
# 			albedo = @SVector[trand(T) for i ∈ 1:3] .* @SVector[trand(T) for i ∈ 1:3]
# 			push!(spheres, Sphere(center, T(0.2), Lambertian(albedo)))
# 		elseif choose_mat < T(0.95)
# 			# metal
# 			albedo = @SVector[random_between(T(0.5),T(1.0)) for i ∈ 1:3]
# 			fuzz = random_between(T(0.0), T(5.0))
# 			push!(spheres, Sphere(center, T(0.2), Metal(albedo, fuzz)))
# 		else
# 			# glass
# 			push!(spheres, Sphere(center, T(0.2), Dielectric(T(1.5))))
# 		end
# 	end

# 	push!(spheres, Sphere((SA{T}[0,1,0]), T(1), Dielectric(T(1.5))))
# 	push!(spheres, Sphere((SA{T}[-4,1,0]), T(1),
# 						  Lambertian(SA{T}[0.4,0.2,0.1])))
# 	push!(spheres, Sphere((SA{T}[4,1,0]), T(1),
# 						  Metal((SA{T}[0.7,0.6,0.5]), T(0))))
# 	HittableList(spheres)
# end

ELEM_TYPE = Float64
t_default_cam = default_camera(SA{ELEM_TYPE}[0,0,0])
render(scene_2_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

t_cam1 = default_camera([13,2,3], [0,0,0], [0,1,0], 20, 16/9, 0.1, 10.0; elem_type=ELEM_TYPE)

print("render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 200, 32):")
_scene_random_spheres = scene_random_spheres(; elem_type=ELEM_TYPE)
render(_scene_random_spheres, t_cam1, 200, 32)
@benchmark render($_scene_random_spheres, $t_cam1, 200, 32)


text:
mutable struct Simulation{T <: SimulationModel}
    status::BUILD_STATUS
    problem::Union{Nothing, SciMLBase.DEProblem}
    tspan::NTuple{2, Float64}
    sys::PSY.System
    perturbations::Vector{<:Perturbation}
    x0_init::Vector{Float64}
    initialized::Bool
    tstops::Vector{Float64}
    callbacks::DiffEqBase.CallbackSet
    simulation_folder::String
    inputs::Union{Nothing, SimulationInputs}
    results::Union{Nothing, SimulationResults}
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
    multimachine::Bool
end

get_system(sim::Simulation) = sim.sys
get_simulation_inputs(sim::Simulation) = sim.inputs
get_initial_conditions(sim::Simulation) = sim.x0_init
get_tspan(sim::Simulation) = sim.tspan

function Simulation(
    ::Type{T},
    sys::PSY.System;
    tspan,
    initial_conditions,
    initialize_simulation,
    perturbations,
    simulation_folder,
    console_level,
    file_level,
) where {T <: SimulationModel}
    PSY.set_units_base_system!(sys, "DEVICE_BASE")

    return Simulation{T}(
        BUILD_INCOMPLETE,
        nothing,
        tspan,
        sys,
        perturbations,
        initial_conditions,
        !initialize_simulation,
        Vector{Float64}(),
        DiffEqBase.CallbackSet(),
        simulation_folder,
        nothing,
        nothing,
        console_level,
        file_level,
        false,
    )
end

function Simulation!(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
) where {T <: SimulationModel}
    return Simulation!(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
end

function Simulation(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
) where {T <: SimulationModel}
    return Simulation(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
end

"""
Builds the simulation object and conducts the indexing process. The initial conditions are stored in the system.

# Accepted Key Words
- `initialize_simulation::Bool : Runs the initialization routine. If false, simulation runs based on the operation point stored in System`
- `initial_conditions::Vector{Float64} : Allows the user to pass a vector with the initial condition values desired in the simulation. If initialize_simulation = true, these values are used as a first guess and overwritten.
- `system_to_file::Bool`: Serializes the initialized system
- `console_level::Logging`: Sets the level of logging output to the console. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
- `file_level::Logging`: Sets the level of logging output to file. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
"""
function Simulation!(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
) where {T <: SimulationModel}
    check_folder(simulation_folder)
    # Instantiates the Simulation object
    sim = Simulation(
        T,
        system;
        tspan = tspan,
        initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
        initialize_simulation = get(kwargs, :initialize_simulation, true),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Debug),
    )

    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "initialized_system.json"))
    end
    return sim
end

"""
Initializes the simulations and builds the indexing. The input system is not modified during the initialization

# Accepted Key Words
- `initialize_simulation::Bool : Runs the initialization routine. If false, simulation runs based on the operation point stored in System`
- `initial_conditions::Vector{Float64} : Initial conditions provided by the user. Must be used with initialize_simulation = false to not be override.`
- `system_to_file::Bool`: Serializes the original input system
- `console_level::Logging`: Sets the level of logging output to the console. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debugg
- `file_level::Logging`: Sets the level of logging output to file. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
"""
function Simulation(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
) where {T <: SimulationModel}
    check_folder(simulation_folder)
    simulation_system = deepcopy(system)
    # Instantiates the Simulation object
    sim = Simulation(
        T,
        simulation_system;
        tspan = tspan,
        initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
        initialize_simulation = get(kwargs, :initialize_simulation, true),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Debug),
    )
    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "input_system.json"))
    end
    return sim
end

function reset!(sim::Simulation{T}) where {T <: SimulationModel}
    @info "Rebuilding the simulation after reset"
    sim.inputs = SimulationInputs(T(), get_system(sim), sim.inputs.tspan)
    build!(sim)
    @info "Simulation reset to status $(sim.status)"
    return
end

function configure_logging(sim::Simulation, file_mode; kwargs...)
    return IS.configure_logging(
        console = true,
        console_stream = stderr,
        console_level = get(kwargs, :console_level, sim.console_level),
        file = true,
        filename = joinpath(sim.simulation_folder, SIMULATION_LOG_FILENAME),
        file_level = get(kwargs, :file_level, sim.file_level),
        file_mode = file_mode,
        tracker = nothing,
        set_global = false,
    )
end

function _build_inputs!(
    sim::Simulation{T},
    frequency_reference,
) where {T <: SimulationModel}
    simulation_system = get_system(sim)
    sim.inputs = SimulationInputs(T, simulation_system, frequency_reference)
    @debug "Simulation Inputs Created"
    return
end

function _get_flat_start(inputs::SimulationInputs)
    bus_count = get_bus_count(inputs)
    var_count = get_variable_count(inputs)
    initial_conditions = zeros(var_count)
    initial_conditions[1:bus_count] .= 1.0
    return initial_conditions
end

function _initialize_state_space(sim::Simulation{T}) where {T <: SimulationModel}
    simulation_inputs = get_simulation_inputs(sim)
    x0_init = sim.x0_init
    if isempty(x0_init) && sim.initialized
        @warn "Initial Conditions set to flat start"
        sim.x0_init = _get_flat_start(simulation_inputs)
    elseif isempty(x0_init) && !sim.initialized
        sim.x0_init = _get_flat_start(simulation_inputs)
    elseif !isempty(x0_init) && sim.initialized
        if length(sim.x0_init) != get_variable_count(simulation_inputs)
            throw(
                IS.ConflictingInputsError(
                    "The size of the provided initial state space does not match the model's state space.",
                ),
            )
        end
    elseif !isempty(x0_init) && !sim.initialized
        @warn "initial_conditions were provided with initialize_simulation. User's initial_conditions will be overwritten."
        if length(sim.x0_init) != get_variable_count(simulation_inputs)
            sim.x0_init = _get_flat_start(simulation_inputs)
        end
    else
        @assert false
    end
    return
end

function _pre_initialize_simulation!(sim::Simulation)
    _initialize_state_space(sim)
    if sim.initialized != true
        @info("Pre-Initializing Simulation States")
        sim.initialized = precalculate_initial_conditions!(sim)
        if !sim.initialized
            error(
                "The simulation failed to find an adequate initial guess for the initialization. Check the intialization routine.",
            )
        end
    else
        @warn(
            "No Pre-initialization conducted. If this is unexpected, check the initialization keywords"
        )
        sim.status = SIMULATION_INITIALIZED
    end
    return
end

function _get_jacobian(sim::Simulation{T}) where {T <: SimulationModel}
    inputs = get_simulation_inputs(sim)
    x0_init = get_initial_conditions(sim)
    return JacobianFunctionWrapper(T(inputs, x0_init, JacobianCache), x0_init)
end

function _build_perturbations!(sim::Simulation)
    @info "Attaching Perturbations"
    if isempty(sim.perturbations)
        @debug "The simulation has no perturbations"
        return DiffEqBase.CallbackSet(), [0.0]
    end
    inputs = get_simulation_inputs(sim)
    perturbations = sim.perturbations
    perturbations_count = length(perturbations)
    callback_vector = Vector{DiffEqBase.DiscreteCallback}(undef, perturbations_count)
    tstops = Vector{Float64}(undef, perturbations_count)
    for (ix, pert) in enumerate(perturbations)
        @debug pert
        condition = (x, t, integrator) -> t in [pert.time]
        affect = get_affect(inputs, get_system(sim), pert)
        callback_vector[ix] = DiffEqBase.DiscreteCallback(condition, affect)
        tstops[ix] = pert.time
    end
    sim.tstops = tstops
    sim.callbacks = DiffEqBase.CallbackSet((), tuple(callback_vector...))
    return
end

function _get_diffeq_problem(
    sim::Simulation,
    model::SystemModel{ResidualModel},
    jacobian::JacobianFunctionWrapper,
)
    x0 = get_initial_conditions(sim)
    dx0 = zeros(length(x0))
    simulation_inputs = get_simulation_inputs(sim)
    sim.problem = SciMLBase.DAEProblem(
        SciMLBase.DAEFunction{true}(
            model;
            # Currently commented for Sundials compatibility
            #jac = jacobian,
            tgrad = (dT, u, p, t) -> dT .= false,
            #jac_prototype = jacobian.Jv,
        ),
        dx0,
        x0,
        get_tspan(sim),
        simulation_inputs,
        differential_vars = get_DAE_vector(simulation_inputs),
    )
    sim.status = BUILT
    return
end

function _get_diffeq_problem(
    sim::Simulation,
    model::SystemModel{MassMatrixModel},
    jacobian::JacobianFunctionWrapper,
)
    simulation_inputs = get_simulation_inputs(sim)
    sim.problem = SciMLBase.ODEProblem(
        SciMLBase.ODEFunction{true}(
            model,
            mass_matrix = get_mass_matrix(simulation_inputs),
            jac = jacobian,
            jac_prototype = jacobian.Jv,
            # Necessary to avoid unnecessary calculations in Rosenbrock methods
            tgrad = (dT, u, p, t) -> dT .= false,
        ),
        sim.x0_init,
        get_tspan(sim),
        simulation_inputs,
    )
    sim.status = BUILT
    return
end

function _build!(sim::Simulation{T}; kwargs...) where {T <: SimulationModel}
    check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
    # Branches are a super set of Lines. Passing both kwargs will
    # be redundant.
    TimerOutputs.reset_timer!(BUILD_TIMER)
    if get(kwargs, :disable_timer_outputs, false)
        TimerOutputs.disable_timer!(BUILD_PROBLEMS_TIMER)
    end

    TimerOutputs.@timeit BUILD_TIMER "Build Simulation" begin
        if get(kwargs, :all_branches_dynamic, false)
            TimerOutputs.@timeit BUILD_TIMER "AC Branch Transform to Dynamic" begin
                sys = get_system(sim)
                transform_branches_to_dynamic(sys, PSY.ACBranch)
            end
        elseif get(kwargs, :all_lines_dynamic, false)
            TimerOutputs.@timeit BUILD_TIMER "Line Transform to Dynamic" begin
                sys = get_system(sim)
                transform_branches_to_dynamic(sys, PSY.Line)
            end
        end
        TimerOutputs.@timeit BUILD_TIMER "Build Simulation Inputs" begin
            f_ref = get(kwargs, :frequency_reference, ReferenceBus)
            _build_inputs!(sim, f_ref)
            sim.multimachine = f_ref != ConstantFrequency
        end
        TimerOutputs.@timeit BUILD_TIMER "Pre-initialization" begin
            _pre_initialize_simulation!(sim)
        end
        if sim.status != BUILD_FAILED
            simulation_inputs = get_simulation_inputs(sim)
            try
                TimerOutputs.@timeit BUILD_TIMER "Calculate Jacobian" begin
                    jacobian = _get_jacobian(sim)
                end
                TimerOutputs.@timeit BUILD_TIMER "Make Model Function" begin
                    model = T(simulation_inputs, get_initial_conditions(sim), SimCache)
                end
                TimerOutputs.@timeit BUILD_TIMER "Initial Condition NLsolve refinement" begin
                    refine_initial_condition!(sim, model, jacobian)
                end
                TimerOutputs.@timeit BUILD_TIMER "Build Perturbations" begin
                    _build_perturbations!(sim)
                end
                TimerOutputs.@timeit BUILD_TIMER "Make DiffEq Problem" begin
                    _get_diffeq_problem(sim, model, jacobian)
                end
                @info "Simulations status = $(sim.status)"
            catch e
                bt = catch_backtrace()
                @error "$T failed to build" exception = e, bt
                sim.status = BUILD_FAILED
            end
        else
            @error "The simulation couldn't be initialized correctly. Simulations status = $(sim.status)"
        end
    end
    return
end

function build!(sim; kwargs...)
    logger = configure_logging(sim, "w")
    try
        Logging.with_logger(logger) do
            _build!(sim; kwargs...)
        end
    catch e
        @error "Build failed" exception = (e, catch_backtrace())
        return sim.status
    finally
        string_buffer = IOBuffer()
        TimerOutputs.print_timer(
            string_buffer,
            BUILD_TIMER,
            sortby = :firstexec,
            compact = true,
        )
        @info "\n$(String(take!(string_buffer)))\n"
        close(logger)
    end
    return
end

function simulation_pre_step!(sim::Simulation, reset_sim::Bool)
    reset_sim = sim.status == CONVERTED_FOR_SMALL_SIGNAL || reset_sim
    if sim.status == BUILD_FAILED
        error(
            "The Simulation status is $(sim.status). Can not continue, correct your inputs and build the simulation again.",
        )
    elseif sim.status != BUILT && !reset_sim
        error(
            "The Simulation status is $(sim.status). Use keyword argument reset_simulation = true",
        )
    end

    reset_sim && reset!(sim)
    return
end

function _prog_meter_enabled()
    return isa(stderr, Base.TTY) || (get(ENV, "CI", nothing) == "true")
end

function _execute!(sim::Simulation, solver; kwargs...)
    @debug "status before execute" sim.status
    simulation_pre_step!(sim, get(kwargs, :reset_simulation, false))
    sim.status = SIMULATION_STARTED
    time_log = Dict{Symbol, Any}()
    solution,
    time_log[:timed_solve_time],
    time_log[:solve_bytes_alloc],
    time_log[:sec_in_gc] = @timed SciMLBase.solve(
        sim.problem,
        solver;
        callback = sim.callbacks,
        tstops = sim.tstops,
        progress = get(kwargs, :enable_progress_bar, _prog_meter_enabled()),
        progress_steps = 1,
        kwargs...,
    )
    if solution.retcode == :Success
        sim.status = SIMULATION_FINALIZED
        sim.results = SimulationResults(
            get_simulation_inputs(sim),
            get_system(sim),
            time_log,
            solution,
        )
    else
        @error("The simulation failed with return code $(solution.retcode)")
        sim.status = SIMULATION_FAILED
    end
end

function execute!(sim::Simulation, solver; kwargs...)
    logger = configure_logging(sim, "a"; kwargs...)
    try
        Logging.with_logger(logger) do
            _execute!(sim, solver; kwargs...)
        end
    catch e
        @error "Execution failed" exception = (e, catch_backtrace())
        return sim.status = SIMULATION_FAILED
    finally
        close(logger)
        return sim.status
    end
end

function read_results(sim::Simulation)
    return sim.results
end

function get_dynamic_wrapper(sim::Simulation, name::AbstractString)
    return get_dynamic_wrapper(get_simulation_inputs(sim), name)
end

function get_setpoints(sim::Simulation)
    return get_setpoints(get_simulation_inputs(sim))
end


text:
module AnalysisTools

using LinearAlgebra
using Plots
using Plots.PlotMeasures
using DataFrames
using Glob
using JSON

using ACEtb.FitTools

export plot_bands, plot_bands!, tabulate_results, comp_plot

function plot_bands(fc::FitContext, phase::Symbol, fermi_level::Float64, bands::Matrix; plotobj=nothing, label=nothing, color=:black)
    @assert phase in [:FCC, :BCC]
    phase = String(phase)

    bandpath = JSON.parsefile(joinpath(fc.base_dir, "$(phase)_tb_bandpath.json"))
    kpoints = bandpath["kpoints"]
    special_points = bandpath["special_points"]

    k0 = 0.0
    x_k = []
    tick_p = [k0]
    tick_k = phase == "FCC" ? ["W"] : tick_k = ["G"]
    push!(x_k, k0)
    for xi=2:size(kpoints, 1)
        k0 += norm(kpoints[xi] .- kpoints[xi-1])
        push!(x_k, k0)
        indict = [k for (k, v) in special_points if v==kpoints[xi] ]
        if length(indict) > 0
            push!(tick_p, k0)
            push!(tick_k, indict[1])
        end
    end

    plotobj === nothing && (plotobj = plot())
    for i = 1:size(bands, 2)
        thislabel = i == 1 ? label : false
        plot!(plotobj, x_k, bands[:, i], color=color, lw=1, label=thislabel)
    end
    for p=1:length(tick_p)
        plot!(plotobj, [tick_p[p]], lw=0.5, c=:gray, seriestype=:vline, label=false)
    end
    plot!(plotobj, [fermi_level], lw=0.5, c=:black, ls=:dash, seriestype=:hline, label=false)
    plot!(plotobj; fg_legend=:transparent, bg_legend=:transparent,
          grid=false, framestyle = :box, legend=:bottomleft)
    xticks!(plotobj, ([t for t in tick_p],[k for k in tick_k]), fontsize=14)
    ylabel!(plotobj, "Energy (eV)")

    xlims!(plotobj, 0.0, maximum(x_k))
    ylims!(plotobj, -20, fermi_level + 10)

    plot(plotobj, dpi=120)
    return plotobj
end

function plot_bands(fc::FitContext, phase::Symbol, param_hash::String; label=nothing, kwargs...)
    fermi_level, bands = read_bands(fc, phase, param_hash)
    label === nothing && (label = param_hash)
    plot_bands(fc, phase, fermi_level, bands; kwargs...,
               label=label)
end

plot_bands!(plotobj::Plots.Plot, args...; kwargs...) = plot_bands(args...; plotobj=plotobj, kwargs...)

function tabulate_results(fc::FitContext)
    df = DataFrame(param_hash=String[],
                   order=Int[],
                   degree=Int[],
                   env_deg=Int[],
                   rcut=Float64[],
                   renv=Float64[],
                   pcut=Int[],
                   N_train=Int[],
                   N_bond=Int[],
                   N_basis=Int[],
                   error_train=Float64[],
                   error_test=Union{Missing,Float64}[],
                   band_error_FCC=Union{Missing,Float64}[],
                   band_error_BCC=Union{Missing,Float64}[],
                   pot_file=String[],
                   training_datasets=String[],
                   method=String[],
                   rtol=Union{Missing,Float64}[])
    for pot_file in glob(joinpath(splitpath(fc.output_dir)[end], "pot_*/pot.json"))
        param_hash = match(r"pot_(.*)/pot.json", pot_file).captures[1]
        pot_dict = JSON.parsefile(pot_file)
        band_error = get(pot_dict, "band_error", missing)
        band_error_FCC = band_error_BCC = missing
        if band_error !== missing
            band_error_FCC = band_error["FCC"]
            band_error_BCC = band_error["BCC"]
        end
        fit_params = pot_dict["model"]["fit_params"]
        cutoff_params = pot_dict["model"]["cutoff_params"]
        row = (param_hash,
               fit_params["order"],
               fit_params["degree"],
               fit_params["env_deg"],
               cutoff_params["rcut"],
               cutoff_params["renv"],
               cutoff_params["pcut"],
               length(pot_dict["training_datasets"]),
               pot_dict["size_A"][1],
               pot_dict["size_A"][2],
               pot_dict["error_train"],
               get(pot_dict, "test_error", missing),
               band_error_FCC,
               band_error_BCC,
               pot_file,
               join(pot_dict["training_datasets"], ","),
               get(fit_params, "method", "QR"),
               get(fit_params, "rtol", missing))
        push!(df, row)
    end

    # add a descriptive label for the types of training configs
    df.config = ifelse.(occursin.("005", df.training_datasets),
                        ifelse.(occursin.("010", df.training_datasets), "mixT", "lowT"),
                        "highT")

    return df
end

function ref_plots(fc::FitContext)
    fcc_fermi_level_FHI, fcc_bands_FHI = read_FHI_bands(fc, :FCC)
    bcc_fermi_level_FHI, bcc_bands_FHI = read_FHI_bands(fc, :BCC)
    fcc_plot = plot_bands(fc, :FCC, fcc_fermi_level_FHI, fcc_bands_FHI; label="FHIaims ref", color=:red)
    bcc_plot = plot_bands(fc, :BCC, bcc_fermi_level_FHI, bcc_bands_FHI; label="FHIaims ref", color=:red)
    display(plot(fcc_plot, bcc_plot, title=["FCC" "BCC"]))
    return fcc_plot, bcc_plot
end

function add_plots(fc::FitContext, potential, color, fcc_plot, bcc_plot; label=nothing)
    plot_bands!(fcc_plot, fc, :FCC, potential, color=color, label=label)
    plot_bands!(bcc_plot, fc, :BCC, potential, color=color, label=label)
    p = plot(fcc_plot, bcc_plot, title=["FCC" "BCC"])
    display(p)
    return p
end

function comp_plot(fc::FitContext, potential; color=:blue, label=nothing)
    fcc_plot, bcc_plot = ref_plots(fc)
    add_plots(fc, potential, color, fcc_plot, bcc_plot, label=label)
end

end


text:
abstract type GloVe{LANG} <: EmbeddingSystem{LANG} end

const glove_max_size = 4_000_000

function init(::Type{GloVe})
    vectors = [("glove.6B",
                "Trained on 6B tokens (Wikipedia 2014, Gigaword 5), 400K vocab, uncased. Includes 50d, 100d, 200d, & 300d vectors. 822 MB download.",
                "617afb2fe6cbd085c235baf7a465b96f4112bd7f7ccb2b2cbd649fed9cbcf2fb",
                ["50d", "100d", "200d", "300d"]),
               ("glove.42B.300d",
                "Trained on 42B tokens (Common Crawl), 1.9M vocab, uncased, Includes 300d vectors. 1.75 GB download.",
                "03d5d7fa28e58762ace4b85fb71fe86a345ef0b5ff39f5390c14869da0fc1970",
               []),
               ("glove.840B.300d",
                "Trained on 840B tokens (Common Crawl), 2.2M vocab, cased. Includes 300d vectors. 2.03 GB download.",
                "c06db255e65095393609f19a4cfca20bf3a71e20cc53e892aafa490347e3849f",
                []),
               ("glove.twitter.27B",
                "Trained on 2B tweets, 27B tokens, 1.2M vocab, uncased, 25d, 50d, 100d, & 200d vectors, 1.42 GB download.",
                "792af52f795d1a32c9842a3240f5f3fe5e941a8ff6df5eb0f9d668092ebc019c",
                ["25d", "50d", "100d", "200d"])]
    for (depname, description, sha, dims) in vectors
        register(DataDep(depname,
                         """
                         Pretrained GloVe word embeddings.
                         Website: https://nlp.stanford.edu/projects/glove/
                         Author: Jeffrey Pennington, Richard Socher, Christopher D. Manning
                         Year: 2014
                         Licence: Open Data Commons Public Domain Dedication and License (PDDL)

                         Paper:
                             Jeffrey Pennington, Richard Socher, Christopher D. Manning. GloVe: Global Vectors for Word Representation. In Proceedings of EMNLP, 2014.

                         $description
                         """,
                         "https://nlp.stanford.edu/data/$(depname).zip",
                         sha,
                         post_fetch_method = unpack))
        if length(dims) >= 1
            append!(language_files(GloVe{:en}), ["$(depname)/$(depname).$(dim).txt" for dim in dims])
        else
            push!(language_files(GloVe{:en}), "$depname/$depname.txt")
        end
    end
end

function _load_embeddings(::Type{<:GloVe}, embedding_file, max_vocab_size, keep_words)
    local LL, indexed_words, index
    if length(keep_words) > 0
        max_vocab_size = length(keep_words)
    end
    indexed_words = Vector{String}()
    LL = Vector{Vector{Float32}}()
    open(embedding_file) do f
        index = 1
        for line in eachline(f)
            xs = split(line, ' ')
            word = xs[1]
            if length(keep_words) == 0 || (word in keep_words)
                index > max_vocab_size && break
                push!(indexed_words, word)
                try
                    push!(LL, parse.(Float32, @view(xs[2:end])))
                catch err
                    err isa ArgumentError || rethrow()
                    @warn "Could not parse word vector" index word exception=err
                end
                index += 1
            end
        end
    end
    return reduce(hcat, LL), indexed_words
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.DAE (ModiaMath/DAE/_module.jl)
#

"""
    initInfo = InitInfo(...)

Data structure that is returned when calling initialize!(...)
"""
struct InitInfo
    y0::Vector{Float64}          # Initial y
    yp0::Vector{Float64}         # Initial yp0
    y_nominal::Vector{Float64}   # Nominal values of y (y_nominal[i]>0.0 required)
    y_errorControl::Vector{Bool} # = true: error control; = false, no error control
    maxTime::Float64             # Integrate at most until maxTime
    nextEventTime::Float64       # Next time event
    integrateToEvent::Bool       # = true, if exact integration to nextEventTime
    terminate::Bool              # = true, terminate simulation

    function InitInfo(y0::Vector{Float64},
                      yp0::Vector{Float64};
                      y_nominal::Vector{Float64}=ones(length(y0)),
                      y_errorControl::Vector{Bool}=fill(true, length(y0)),
                      maxTime::Float64=floatmax(Float64),
                      nextEventTime::Float64=floatmax(Float64),
                      integrateToEvent::Bool=true,
                      terminate::Bool=false)
        @assert(length(y0) == length(yp0))
        @assert(length(y0) == length(y_nominal))
        @assert(y_nominal[ argmin(y_nominal) ] > 0.0)

        new(y0, yp0, y_nominal, y_errorControl, maxTime, nextEventTime, integrateToEvent, terminate)
    end
end

"""
    eventInfo = EventInfo(...)

Data structure in which the return values of processEvent!(...) are stored
"""
mutable struct EventInfo
    restart::EventRestart
    maxTime::Float64         # Integrate at most until maxTime
    nextEventTime::Float64   # Next time event
    integrateToEvent::Bool   # = true, if exact integration to nextEventTime
    nextModel                # If restart==FullRestart, next model

    EventInfo() = new(Restart, floatmax(Float64), floatmax(Float64), true, nothing)
end

function reset!(einfo::EventInfo)
    einfo.restart = Restart
    einfo.maxTime = floatmax(Float64)
    einfo.nextEventTime = floatmax(Float64)
    einfo.integrateToEvent = true
    einfo.nextModel = nothing
end