text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module CantorMachines

export ModuleCantorMachines
export CantorMachine, CantorEnumeration, CantorPairing
export CantorBoustrophedonicMachine
export CantorBoustrophedonicEnumeration
export CantorBoustrophedonicPairing
export RosenbergStrongBoustrophedonicMachine
export RosenbergStrongBoustrophedonicEnumeration
export RosenbergStrongBoustrophedonicPairing
export V319514, L319514

"""

[Cantor's enumeration of N X N revisited](https://luschny.wordpress.com/2018/09/24/cantors-enumeration-of-n2-revisited/).

* Cantor-Machine, Cantor-Enumeration, Cantor-Pairing, Cantor-BoustrophedonicMachine, Cantor-BoustrophedonicEnumeration, Cantor-BoustrophedonicPairing, RosenbergStrong-BoustrophedonicMachine, RosenbergStrong-BoustrophedonicEnumeration, RosenbergStrong-BoustrophedonicPairing
"""
const ModuleCantorMachines = ""

# https://luschny.wordpress.com/2018/09/24/cantors-enumeration-of-n2-revisited/
# https://luschny.wordpress.com/2018/09/25/the-cantor-jump-machine/
# https://oeis.org/search?q=A319514&go=Search

"""

The Cantor enumeration implemented as a state machine to avoid the evaluation of the square root function.
"""
function CantorMachine(x, y, state)
    x == 0 && !state && return x, y + 1, !state
    y == 0 && state && return x + 1, y, !state
    state && return x + 1, y - 1, state
    return x - 1, y + 1, state
end

"""

The  Cantor enumeration of ℕ X ℕ where ℕ ``= {0, 1, 2, ...}``. If ``(x, y)`` and ``(x', y')`` are adjacent points on the trajectory of the map then max``(|x - x'|, |y - y'|)`` can become arbitrarily large. In this sense Cantor's enumeration is not continous.
"""
function CantorEnumeration(len)
    x, y, state = 0, 0, false
    for n ∈ 0:len
        #println("$n -> ($x, $y)")
        print("$x, $y, ")
        x, y, state = CantorMachine(x, y, state)
    end
end

"""

The inverse function of the Cantor enumeration (the pairing function), computes n for given ``(x, y)`` and returns ``(x + y)(x + y + 1)/2 + p`` where ``p = x`` if ``x - y`` is odd and ``y`` otherwise.
"""
function CantorPairing(x, y)
    p = isodd(x - y) ? x : y
    div((x + y) * (x + y + 1), 2) + p
end

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

"""

The boustrophedonic Cantor enumeration implemented as a state machine to avoid the evaluation of the square root function.
"""
function CantorBoustrophedonicMachine(x, y)
    x == 0 && return y >> 1 + 1, y - y >> 1
    x > y && return y, x
    return y + 1, x - 1
end

"""

The boustrophedonic Cantor enumeration of ℕ X ℕ where ``ℕ = {0, 1, 2, ...}``. If ``(x, y)`` and ``(x', y')`` are adjacent points on the trajectory of the map then max``(|x - x'|, |y - y'|)`` is always 1 whereas for the Cantor enumeration this quantity can become arbitrarily large. In this sense the boustrophedonic variant is continuous whereas Cantor's realization is not.
"""
function CantorBoustrophedonicEnumeration(len)
    x, y = 0, 0
    for n ∈ 0:len
        # println("$n -> ($x, $y)")
        print("$x, $y, ")
        x, y = CantorBoustrophedonicMachine(x, y)
    end
end

"""

The inverse function of the boustrophedonic Cantor enumeration (the pairing function), computes n for given ``(x, y)`` and returns ``(x + y)(x + y + 1)/2 + m`` where ``m = |x - y| - (x > y  ?  1 : 0)``.
"""
function CantorBoustrophedonicPairing(x, y)
    m = abs(x - y) - (x > y ? 1 : 0)
    div((x + y) * (x + y + 1), 2) + m
end

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

"""

The boustrophedonic Rosenberg-Strong enumeration as considered by Pigeon implemented as a state machine to avoid the evaluation of the square root function.
"""
function RosenbergStrongBoustrophedonicMachine(x, y, state)
    x == 0 && state == 0 && return x, y + 1, 1
    y == 0 && state == 2 && return x + 1, y, 3
    x == y && state == 1 && return x, y - 1, 2
    x == y && return x - 1, y, 0
    state == 0 && return x - 1, y, 0
    state == 1 && return x + 1, y, 1
    state == 2 && return x, y - 1, 2
    return x, y + 1, 3
end

"""

The boustrophedonic Rosenberg-Strong enumeration of ℕ X ℕ where ℕ ``= {0, 1, 2, ...}``. If ``(x, y)`` and ``(x', y')`` are adjacent points on the trajectory of the map then max``(|x - x'|, |y - y'|)`` is always 1 whereas the Rosenberg-Strong realization is not.
"""
function RosenbergStrongBoustrophedonicEnumeration(len)
    x, y, state = 0, 0, 0
    for n ∈ 0:len
        #println("$n -> ($x, $y)")
        print("$x, $y, ")
        x, y, state = RosenbergStrongBoustrophedonicMachine(x, y, state)
    end
end

"""

The inverse function of the boustrophedonic Rosenberg-Strong enumeration (the pairing function), computes ``n`` for given ``(x, y)``.
"""
function RosenbergStrongBoustrophedonicPairing(x::Int, y::Int)
    m = max(x, y)
    d = isodd(m) ? x - y : y - x
    m * (m + 1) + d
end

"""

Return the pair ``(x, y)`` for given n as given by the boustrophedonic Rosenberg-Strong enumeration.
"""
function V319514(n)
    k, r = Base.divrem(n, 2)
    m = x = isqrt(k)
    y = k - x^2
    x <= y && ((x, y) = (2x - y, x))
    isodd(m) ? (y, x)[r+1] : (x, y)[r+1]
end

"""

Return a list of pairs ``(x, y)`` given by the boustrophedonic Rosenberg-Strong enumeration.
"""
L319514(len) = [V319514(n) for n ∈ 0:len-1]

#START-TEST-########################################################

using Test

function test()

    @testset "CantorMachines" begin

        len = 64

        # CantorTest

        println()
        println("Running CantorMachine")

        x, y, state = 0, 0, false
        for n ∈ 0:len
            p = CantorPairing(x, y)
            #println("$n -> ($x, $y) -> $p")
            @test n == p
            x, y, state = CantorMachine(x, y, state)
        end

        # CantorBoustrophedonicTest

        #println()
        println("Running CantorBoustrophedonicMachine")

        x, y = 0, 0
        for n ∈ 0:len
            p = CantorBoustrophedonicPairing(x, y)
            #println("$n -> ($x, $y) -> $p")
            @test n == p
            x, y = CantorBoustrophedonicMachine(x, y)
        end

        # RosenbergStrongBoustrophedonicTest

        #println()
        println("Running RosenbergStrongBoustrophedonicMachine")

        x, y, state = 0, 0, 0
        for n ∈ 0:len
            p = RosenbergStrongBoustrophedonicPairing(x, y)
            #println("$n -> ($x, $y) -> $p")
            @test n == p
            x, y, state = RosenbergStrongBoustrophedonicMachine(x, y, state)
        end

        println()
    end
end

function demo()

    println()
    println("CantorEnumeration")
    CantorEnumeration(20)

    println()
    println()
    println("CantorBoustrophedonicEnumeration")
    CantorBoustrophedonicEnumeration(20)

    println()
    println()
    println("RosenbergStrongBoustrophedonicEnumeration")
    RosenbergStrongBoustrophedonicEnumeration(20)

    println()
    println()
    println("L319514")
    L319514(42) |> println
end

function perf() end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
using HCIDatasets
using HCIDatasets: BetaPictoris, HR8799, V471Tau
using Test

macro HCIDatasetTestset(dataset, keys, sizes, T=Float32)
    quote
        @testset $(string(dataset)) begin
            @test keys($dataset) == $keys
            @test $dataset[$keys] == $dataset[$keys...]
            for ((k, arr), key, sz) in zip(pairs($dataset), $keys, $sizes)
                @test k == key
                @test $dataset[key] == arr
                @test size(arr) == sz
                @test eltype(arr) == $T
            end
        end
    end
end

@HCIDatasetTestset BetaPictoris (:cube, :cube_empty, :pa, :psf) [(101, 101, 61), (101, 101, 61), (61,), (39, 39)]
@HCIDatasetTestset HR8799 (:cube, :pa) [(501, 501, 74), (74,)]
@HCIDatasetTestset V471Tau (:cube, :pa, :psf, :wl) [(132, 132, 24, 39), (24,), (25, 25, 39), (39,)]

ext: jl


text:
module VulkanGen

using StructArrays
using UnPack
using LightGraphs
using MLStyle
using DocStringExtensions

@template (FUNCTIONS, METHODS, MACROS) = """
                                         $(DOCSTRING)
                                         $(TYPEDSIGNATURES)
                                         """

@template TYPES = """
                  $(DOCSTRING)
                  $(TYPEDEF)
                  $(TYPEDSIGNATURES)
                  $(TYPEDFIELDS)
                  $(SIGNATURES)
                  """

include("types.jl")
include("spec.jl")
include("wrapper.jl")

end # module VulkanGen


text:
using BoardGames
using BoardGamesGames
using BoardGamesGames.TicTacToeGame
using BoardGamesGames.OthelloGame
using Test

b = TicTacToeBoard()
for i in 1:7
    global b = play(b, i)
end

function randomgame(board)
    while !isempty(getmoves(board))
        m = rand(getmoves(board))
        board = play(board, m)
    end
    board
end

othello_test_board = zeros(Int, (8,8))
othello_test_board[8,6] =  1
othello_test_board[7,7] = -1
othello_test_board[6,8] = 1
othello_test_board[8,7] = -1

othello_test_board = OthelloBoard(othello_test_board, -1)

new_othello = OthelloBoard()

function Base.:(==)(b1::OthelloBoard, b2::OthelloBoard)
    return b1.v == b2.v && b1.turn == b2.turn
end

@testset "BoardGamesGames.jl" begin
    # Write your tests here.
    @test TicTacToeBoard().turn == 1
    @test getmoves(TicTacToeBoard()) == 1:9
    @test getmoves(play(TicTacToeBoard(), 5)) == (1:4) ∪ (6:9)
    @test playerturn(play(TicTacToeBoard(), 5)) == 2
    @test b.winner == 1
    @test winner(b) == 1
    @test show(stdout, MIME"image/png"(), b) === nothing
    @test show(b) === nothing
    @test 0 <= winner(randomgame(TicTacToeBoard())) <= 2

    @test playerturn(OthelloBoard()) == 1
    @test playerturn(play(OthelloBoard(), (3,4))) == 2
    @test show(OthelloBoard()) === nothing
    @test 0 <= winner(randomgame(OthelloBoard())) <= 2
    @test getmoves(othello_test_board) == [(8,5)]
    @test winner(play(othello_test_board, (8,5))) == 2
    @test show(stdout, MIME"image/png"(), othello_test_board) === nothing
    @test show(othello_test_board) === nothing
    @test play!(new_othello, (3,4)) == play(OthelloBoard(), (3,4))

end


text:
using CoDa
using CSV
using RData
using Tables
using DataDeps
using Distances
using ScientificTypes
using TableTransforms
using LinearAlgebra
using Test, Random, Plots
using ReferenceTests, ImageIO

# workaround for GR warnings
ENV["GKSwstype"] = "100"

# accept downloads without interaction
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

# download and setup data dependencies
register(DataDep("juraset",
  "A geochemical dataset from the Swiss Jura",
  "https://github.com/cran/compositions/raw/master/data/juraset.rda",
  "0fcf23fbca1d6fcb58ae0de6f11365f39fa3df02828128708185ebf45b002382"))
rda  = joinpath(datadep"juraset", "juraset.rda")
jura = RData.load(rda)["juraset"]
CSV.write(joinpath(datadir,"jura.csv"), jura)

# list of tests
testfiles = [
  "compositions.jl",
  "codaarrays.jl",
  "scitypes.jl",
  "distances.jl",
  "transforms.jl",
  "covariances.jl",
  "matrices.jl"
]

@testset "CoDa.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
using Optim
using GR
using Plots

gr()

f(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

start = [-5.0,4.0]
iter = 80

GD_res = optimize(f, start, GradientDescent(),
               Optim.Options(g_tol = 1e-40,
                             iterations = iter,
                             store_trace = true,
                             extended_trace = true,
                             show_trace = false))

CGD_res = optimize(f, start, ConjugateGradient(),
               Optim.Options(g_tol = 1e-40,
                             iterations = iter,
                             store_trace = true,
                             extended_trace = true,
                             show_trace = false))

SA_res = optimize(f, start, SimulatedAnnealing(),
               Optim.Options(g_tol = 1e-40,
                             iterations = iter,
                             store_trace = true,
                             extended_trace = true,
                             show_trace = false))

l = iter + 1

GD_trace = Array{Float64}(undef, l, 2)
CGD_trace = Array{Float64}(undef, l, 2)
SA_trace = Array{Float64}(undef, l, 2)

for i in 1:l
    x1 = Optim.x_trace(GD_res)[i][1]
    x2 = Optim.x_trace(GD_res)[i][2]
    GD_trace[i,1] = x1
    GD_trace[i,2] = x2
end

for i in 1:l
    x1 = Optim.x_trace(CGD_res)[i][1]
    x2 = Optim.x_trace(CGD_res)[i][2]
    CGD_trace[i,1] = x1
    CGD_trace[i,2] = x2
end

for i in 1:l
    x1 = Optim.x_trace(SA_res)[i][1]
    x2 = Optim.x_trace(SA_res)[i][2]
    SA_trace[i,1] = x1
    SA_trace[i,2] = x2
end

contour(-5:0.005:5,-5:0.005:5,f)

plot!([GD_trace[:,1]],
    [GD_trace[:,2]],
    lab = "Gradient Descent",
    zcolor = reverse(GD_trace, dims = 1), m = (4, 0.9, :reds, Plots.stroke(0)),
    cbar = false,
    w = 1.5,
    linecolor = :red)

plot!([CGD_trace[:,1]],
    [CGD_trace[:,2]],
    lab = "Conjugated GD",
    zcolor = reverse(CGD_trace, dims = 1), m = (4, 0.9, :greens, Plots.stroke(0)),
    cbar = false,
    w = 1.5,
    linecolor = :green)

plot!([SA_trace[:,1]],
    [SA_trace[:,2]],
    lab = "Simulated Annealing",
    zcolor = reverse(SA_trace, dims = 1), m = (4, 0.9, :blues, Plots.stroke(0)),
    cbar = false,
    w = 1.5,
    linecolor = :blue)

max_stars_repo_name: lukaskln/Optim-Visualization

max_issues_repo_name: lukaskln/Optim-Visualization

max_forks_repo_name: lukaskln/Optim-Visualization


text:
function show(io::IO, evar::SingleLevel)
    print(
		io,
		"The Single-Level Variable \"$(evar.varID)\" has the following properties:\n",
		"    Variable ID    (varID) : ", evar.varID, '\n',
		"    Long Name      (lname) : ", evar.lname, '\n',
		"    Variable Name  (vname) : ", evar.vname, '\n',
		"    Variable Units (units) : ", evar.units, '\n',
	)
end

function show(io::IO, evar::PressureLevel)
    print(
		io,
		"The Pressure-Level Variable \"$(evar.varID)\" has the following properties:\n",
		"    Variable ID    (varID) : ", evar.varID, '\n',
		"    Long Name      (lname) : ", evar.lname, '\n',
		"    Variable Name  (vname) : ", evar.vname, '\n',
		"    Variable Units (units) : ", evar.units, '\n',
		"    Pressure Level  (hPa)  : ", evar.hPa,   '\n',
	)
end


text:
using RandomizedLinAlg
using Test, LinearAlgebra, Random

@testset "Randomized Linear Algebra" begin

    Random.seed!(1234321)

    m = n = 100
    B = randn(m, n)
    nB = opnorm(B)
    p = 1e-5 #Probability of failure

    for j = 1 : n
        @test rnorm(B, 2j, p) ≥ nB
        @test rnorms(B, j, p) ≥ nB
    end

    A = B * B'
    k = 1
    l, u = reigmin(A, k, p)
    @test l ≤ eigmin(A) ≤ u

    l, u = reigmax(A, k, p)
    @test l ≤ eigmax(A) ≤ u

    l, u = rcond(A, k, p)
    @test l ≤ cond(A) ≤ u

    @test_throws ArgumentError RandomizedLinAlg.randnn(Char, m)
    @test_throws ArgumentError RandomizedLinAlg.randnn(Char, m, n)
end


text:
using LinearAlgebra
import JSON3

function flatten_json(d, prefix_delim = ".")
    new_d = empty(d)
    for (key, value) in d
        if isa(value, JSON3.Object)
            flattened_value = flatten_json(value, prefix_delim)
            for (ikey, ivalue) in flattened_value
                new_d[Symbol(Symbol(key,".",ikey))] = ivalue
            end
        else
            new_d[key] = value
        end
    end
    return new_d
end

"""
p1 = (48.0929414, 11.5355774)
p2 = (48.0428369, 11.5109936)

# 51.5, 0, 38.8, -77.1
p1 = (51.5, 0)
p2 = (38.8, -77.1)
distance(p1, p2) # 5918.185064088764

p1 = (0,0)
p2 = (0,0)

distanceInKmBetweenEarthCoordinates(51.5, 0, 38.8, -77.1)

lat1, lon1, lat2, lon2 = 51.5, 0, 38.8, -77.1
"""
function distanceInKmBetweenEarthCoordinates(lat1, lon1, lat2, lon2)
    earthRadiusKm = 6371

    dLat = deg2rad(lat2-lat1) # -0.2216568150032799
    dLon = deg2rad(lon2-lon1) # -1.345648853287628

    rlat1 = deg2rad(lat1) # 0.8988445647770797
    rlat2 = deg2rad(lat2) # 0.6771877497737998

    a = sin(dLat/2) * sin(dLat/2) + # 0.012232728027071711
        sin(dLon/2) * sin(dLon/2) * # 0.3883749419945243
        cos(rlat1) * cos(rlat2) # 0.48514929005721846
        # 0.20065255541172858

    c = 2 * atan(sqrt(a), sqrt(1-a)) # 0.928925610436158

    return earthRadiusKm * c # 5918.185064088762
end

function distance(p1::Tuple, p2::Tuple; earthRadiusKm=6371)
    pdist = deg2rad.(p2 .- p1) # (-0.2216568150032799, -1.345648853287628)
    rlat = deg2rad.((p1[1], p2[1])) # (0.8988445647770797, 0.6771877497737998)

    # halfLapsAroundGlobe
    a = sum(sin.(pdist ./ 2) .^ 2) * dot(cos.(rlat)...)

    return (2 * atan(sqrt(a), sqrt(1 - a))) * earthRadiusKm
end


text:
using CasesLH, Test

mdl = CasesLH;

function modifier_test()
    @testset "Modifier" begin
        m1 = Modifier((:a, 1));
        s1 = make_string(m1);
        @test isequal(s1, "a1");

        m2 = Modifier((:a, 1.0));
        s2 = make_string(m2);
        @test isequal(s2, "a1p0");

        m4 = Modifier((:entryPrefScale, 2));
        @test isequal(make_string(m4), "entryPrefScale2");
        @test isequal(m4[2], 2);

        m3 = Modifier(:a);
        @test isequal(make_string(m3), "a");

        @test isequal(make_modifier((:z, 2)), Modifier((:z, 2)));
        @test isequal(make_modifier(m3), m3);
    end
end

@testset "Modifier" begin
    modifier_test();
end

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


text:
## Functions to handle input files.

using Chemfiles

include("clustering.jl")

"""
    parse_cif(file_path)

Parse a CIF file and return a dictionary where each identifier (without the
starting '_' character) is linked to its value.
Values are either a string or a vector of string (if defined in a loop).
"""
function parse_cif(file)
    all_data = Dict{String, Union{String, Vector{String}}}()
    inloop = false
    loopisspecified = false
    loopspec = String[]
    loop_n = 0

    l = read(file, String)
    i, j, x = nextword(l, 0)
    while i != 0

        if inloop
            if !loopisspecified
                if l[i] != '_' # This indicates the start of the values
                    loop_n = length(loopspec)
                    loopisspecified = true
                else # The identifier is part of the loop specification
                    push!(loopspec, l[i+1:j])
                    all_data[l[i+1:j]] = String[]
                    i, j, x = nextword(l, x); continue
                end
            end

            # From this point, the loop has been specified
            @assert loopisspecified
            if l[i] != '_' && l[i:j] != "loop_"
                for k in 1:loop_n
                    push!(all_data[loopspec[k]], l[i:j])
                    i, j, x = nextword(l, x)
                end
                continue
            end
            if l[i:j] == "loop_"
                loopisspecified = false
                loopspec = String[]
                i, j, x = nextword(l, x); continue
            end

            # This point can only be reached if we just quitted a loop
            inloop = false
        end

        @assert !inloop
        if l[i] == '_' # Simple identifier definition
            next_i, next_j, next_x = nextword(l, x)
            @assert next_i != 0
            all_data[l[i+1:j]] = l[next_i:next_j]
            x = next_x
        else
            if l[i:j] == "loop_"
                inloop = true
                loopisspecified = false
                loopspec = String[]
            elseif j-i > 4 && l[i:i+4] == "data_"
                @assert !haskey(all_data, "data")
                all_data["data"] = l[i+5:j]
            else
                k::Int = findprev(isequal('\n'), l, i)
                n::Int = count("\n", l[1:k])
                error("Unkown word \"$(l[i:j])\" at line $(n+1), position $(i-k):$(j-k)")
            end
        end

        i, j, x = nextword(l, x)
    end

    return all_data
end

function parsestrip(s)
    s = s[end] == ')' ? s[1:prevind(s, findlast('(', s))] : s
    return parse(BigFloat, s)
end

"""
    CIF(file_path::AbstractString)

Make a CIF object out of the parsed file.
"""
CIF(file::AbstractString) = CIF(parse_cif(file))
function CIF(parsed::Dict{String, Union{Vector{String},String}})
    natoms = length(parsed["atom_site_label"])
    equivalentpositions = pop!(parsed,
        haskey(parsed, "symmetry_equiv_pos_as_xyz") ?
            "symmetry_equiv_pos_as_xyz" : "space_group_symop_operation_xyz")
    refid = find_refid(equivalentpositions)
    cell = Cell(Symbol(haskey(parsed, "symmetry_cell_setting") ?
                            pop!(parsed, "symmetry_cell_setting") : ""),
                pop!(parsed, "symmetry_space_group_name_H-M"),
                haskey(parsed, "symmetry_Int_Tables_number") ?
                    parse(Int, pop!(parsed, "symmetry_Int_Tables_number")) : 0,
                parsestrip(pop!(parsed, "cell_length_a")),
                parsestrip(pop!(parsed, "cell_length_b")),
                parsestrip(pop!(parsed, "cell_length_c")),
                parsestrip(pop!(parsed, "cell_angle_alpha")),
                parsestrip(pop!(parsed, "cell_angle_beta")),
                parsestrip(pop!(parsed, "cell_angle_gamma")),
                parse.(EquivalentPosition, equivalentpositions, Ref(refid)))

    haskey(parsed, "symmetry_equiv_pos_site_id") && pop!(parsed, "symmetry_equiv_pos_site_id")
    removed_identity = false
    for i in eachindex(cell.equivalents)
        eq = cell.equivalents[i]
        if isone(eq.mat) && iszero(eq.ofs)
            deleteat!(cell.equivalents, i)
            removed_identity = true
            break
        end
    end
    @assert removed_identity

    labels =  pop!(parsed, "atom_site_label")
    _symbols = haskey(parsed, "atom_site_type_symbol") ?
                    pop!(parsed, "atom_site_type_symbol") : copy(labels)
    symbols = String[]
    @inbounds for x in _symbols
        i = findfirst(!isletter, x)
        push!(symbols, isnothing(i) ? x : x[1:i-1])
    end
    pos_x = pop!(parsed, "atom_site_fract_x")
    pos_y = pop!(parsed, "atom_site_fract_y")
    pos_z = pop!(parsed, "atom_site_fract_z")

    types = Symbol[]
    pos = Matrix{Float64}(undef, 3, natoms)
    correspondence = Dict{String, Int}()
    for i in 1:natoms
        @assert !haskey(correspondence, labels[i])
        correspondence[labels[i]] = i
        push!(types, Symbol(symbols[i]))
        pos[:,i] = parsestrip.([pos_x[i], pos_y[i], pos_z[i]])
        pos[:,i] .-= floor.(Int, pos[:,i])
    end

    invids = sortperm(types)
    types = types[invids]
    ids = invperm(invids)
    bonds = falses(natoms, natoms)
    if haskey(parsed, "geom_bond_atom_site_label_1") &&
       haskey(parsed, "geom_bond_atom_site_label_2")
        bond_a = pop!(parsed, "geom_bond_atom_site_label_1")
        bond_b = pop!(parsed, "geom_bond_atom_site_label_2")
        for i in 1:length(bond_a)
            x = correspondence[bond_a[i]]
            y = correspondence[bond_b[i]]
            bonds[x,y] = bonds[y,x] = 1
        end
    end

    return CIF(parsed, cell, ids, types, pos, bonds)
end

"""
    parse_arc(file)

Parse a .arc Systre archive such as the one used by the RCSR.
Return a pair `(flag, pairs)`.

`flag` is set if the archive corresponds to one generated by a compatible release
of CrystalNets. If unset, the genomes of the archive may not be the same as those
computed by CrystalNets for the same nets.
`pairs` is a `Dict{String,String}` whose entries have the form `genome => id` where `id`
is the name of the net and `genome` is the topological genome corresponding to this net
(given as a string of whitespace-separated values parseable by `PeriodicGraph`).
"""
function parse_arc(file)
    pairs = Dict{String,String}()
    curr_key = ""
    counter = 1
    firstline = readline(file)
    flag = startswith(firstline, "Made by CrystalNets.jl v")
    if flag
        flag = CRYSTAL_NETS_VERSION == VersionNumber(firstline[25:end])
    end
    for l in eachline(file)
        if length(l) > 3 && l[1:3] == "key"
            @assert isempty(curr_key)
            i = 4
            while isspace(l[i])
                i += 1
            end
            curr_key = l[i:end]
            @assert !haskey(pairs, curr_key)
        elseif length(l) > 2 && l[1:2] == "id"
            @assert !isempty(curr_key)
            i = 3
            while isspace(l[i])
                i += 1
            end
            pairs[curr_key] = l[i:end]
            curr_key = ""
        end
    end
    return flag, pairs
end

"""
    parse_arcs(file)

Parse a folder containing .arc Systre archives such as the one used by the RCSR.
Return a pair `(flag, pairs)` with the same convention than `parse_arc`
"""
function parse_arcs(path)
    combine(x, y) = x * ", " * y
    dict = Dict{String,String}()
    flag = true
    for f in readdir(path; sort=true)
        _flag, _dict = parse_arc(path*f)
        flag &= _flag
        mergewith!(combine, dict, _dict)
    end
    return flag, dict
end

function parse_atom_name(name::AbstractString)
    firstsep = findfirst(x -> ispunct(x) || isspace(x) || isnumeric(x), name)
    if firstsep isa Nothing
        return String(name)
    else
        firstsep::Int
        if firstsep != 1 && !any(isuppercase, @view name[nextind(name, firstsep):end])
            return String(name[1:prevind(firstsep)])
        else
            return String(name)
        end
    end
end

function parse_atom(name)
    atom = Chemfiles.Atom(name)
    set_type!(atom, parse_atom_name(name))
    return atom
end

function set_unique_bond_type!(frame::Frame, types, bond_length, bonded_atoms::Tuple{Symbol, Symbol}, onlykeep, tol)
    @ifwarn @info "To avoid guessing bonds, use a file format that contains the bonds."
    n = Int(size(frame))
    pos = Chemfiles.positions(frame)
    mat = Chemfiles.matrix(UnitCell(frame))'
    indices = [i for i in 1:n if types[i] ∈ onlykeep]
    #=@inbounds=# for _i in 1:length(indices)
        i = indices[_i]
        for _j in _i+1:length(indices)
            j = indices[_j]
            if minmax(types[i], types[j]) == bonded_atoms
                bonded = abs2(periodic_distance(pos[:,i], pos[:,j], mat) - bond_length) <= tol
                if bonded
                    Chemfiles.add_bond!(frame, i-1, j-1)
                end
            end
        end
    end
    nothing
end

function try_guess_bonds!(frame::Frame)
    # n = Int(size(frame))
    # types = Vector{Symbol}(undef, n)
    # for i in 1:n
    #     types[i] = Symbol(parse_atom_name(Chemfiles.type(Chemfiles.Atom(frame, i-1))))
    # end
    # unique_types = unique!(sort(types))
    # if unique_types == [:C] || unique_types == [:C, :H]
    #     if !(NOWARN::Bool)
    #         @warn "Guessing bonds. The structure seems to be made of only carbons (and possibly hydrogens): using an interatomic distance of 1.54±0.3 Å to assign edges between C atoms."
    #     end
    #     set_unique_bond_type!(frame, types, 1.54, (:C, :C), (:C,), 0.3)
    # elseif unique_types == [:O, :Si] || unique_types == [:Si]
    #     if !(NOWARN::Bool)
    #         @warn "Guessing bonds. The structure seems to be a zeolite: using an interatomic distance of 3.1±0.2 Å to assign edges between Si atoms."
    #     end
    #     set_unique_bond_type!(frame, types, 1.63, (:O, :Si), (:O, :Si), 0.1)
    # else
        @ifwarn begin
            @warn "Guessing bonds through Chemfiles. This may take a while for big structures and may be inexact."
            @info "To avoid guessing bonds, use a file format that contains the bonds."
        end
        guess_bonds!(frame)
    # end
end

function attribute_residues(topology, assert_use_existing_residues)
    m = Int(count_residues(topology))
    n = Int(size(topology))
    atoms_in_residues = m == 0 ? 0 : sum(length(atoms(Residue(topology, i))) for i in 0:(m-1))
    @assert atoms_in_residues <= n

    if atoms_in_residues < n
        if assert_use_existing_residues
            throw(ArgumentError("""
            Cannot use existing residues as vertices because not all atoms have an associated residue.
            To fix this, either assign a residue to each atom or provide another way to detect the vertices.
            """))
        end
        @ifwarn begin
            if atoms_in_residues > 0
                @warn "Some but not all atoms have an associated residue, so we cannot rely on existing residues"
            end
        end
        attributions = Int[]
    else
        attributions = zeros(Int, n)
        for i_residue in 1:m
            for atom in atoms(Residue(topology, i_residue-1))
                attributions[atom+1] = i_residue
            end
        end
    end
    return attributions
end

@static if !isdefined(Chemfiles, :atoms) # up to 0.9.3 included
    function atoms(residue::Residue)
        count = size(residue)
        result = Array{UInt64}(undef, count)
        Chemfiles.__check(Chemfiles.lib.chfl_residue_atoms(Chemfiles.__const_ptr(residue), pointer(result), count))
        return result
    end
end

function check_collision(pos, mat)
    _n, n = size(pos)
    @assert _n == 3
    invmat = inv(mat)
    toremove = Int[]
    for i in 1:n, j in (i+1):n
        if periodic_distance(invmat*pos[:,i], invmat*pos[:,j], mat) < 0.55
            push!(toremove, i)
            push!(toremove, j)
        end
    end
    tokeep = collect(1:n)
    if !isempty(toremove)
        @ifwarn @warn "This file contains multiple colliding atoms. All colliding atoms will be removed."
        unique!(sort!(toremove))
    end
    return toremove
end

hasHneighbor(types, graph, i) = any(x -> types[x.v] === :H, neighbors(graph,i))
function check_valence(types, graph)
    n = length(types)

    ## Small atoms valence check
    invalidatoms = Set{Symbol}()
    for i in 1:n
        d = degree(graph, i)
        t = types[i]
        if t === :O && d > 2 && (d > 3 || !hasHneighbor(types, graph, i))
            push!(invalidatoms, :O)
        elseif t === :C && d > 4
            push!(invalidatoms, :C)
        elseif t === :N && d > 3 && (d > 4 || !hasHneighbor(types, graph, i))
            push!(invalidatoms, :N)
        end
    end
    if !isempty(invalidatoms)
        s = String.(collect(invalidatoms))
        @ifwarn @warn "Found $(join(s, ',', " and ")) with too many bonds."
        return true
    end
    return false
end

"""
    @enum BondingMode

Selection mode for the detection of bonds. The choices are:
-   `InputBonds`: use the input bonds. Fail if those are not specified.
-   `ChemfilesBonds`: use chemfiles built-in bond detection mechanism.
-   `AutoBonds`: if the input specifies bonds, use them unless they look suspicious (too or too
    large according to a heuristic). Otherwise, fall back to `ChemfilesBonds`.
"""
@enum BondingMode begin
    InputBonds
    ChemfilesBonds
    AutoBonds
end

function sanity_checks!(pos, types, graph, mat, bondingmode)
    ## Bond length check
    removeedges = PeriodicEdge3D[]
    for e in edges(graph)
        s, d = e.src, e.dst.v
        bondlength = norm(pos[:,d] .+ (mat * e.dst.ofs) .- pos[:,s])
        if bondlength > 3
            @ifwarn @warn "Suspiciously large bond found."
            return true
        elseif bondlength < 0.85 && types[s] !== :H && types[d] !== :H
            push!(removeedges, e)
        end
    end
    if bondingmode == AutoBonds
        @ifwarn begin
            if !isempty(removeedges)
                @warn "Suspiciously small bond found. Such bonds are probably spurious and will be deleted."
                @info "To force retaining these bonds, use --bond-detect=input or --bond-detect=chemfiles"
            end
        end
        for e in removeedges
            rem_edge!(graph, e)
        end
    end
    return false
end

"""
       parse_chemfiles(path)

Parse a file given in any reckognised chemical format and extract the topological
information.
Such format can be .cif or any file format reckognised by Chemfiles.jl that
contains all the necessary topological information.
"""
function parse_chemfile(_path, exportto=tempdir(), bondingmode::BondingMode=AutoBonds, assert_use_existing_residues=false ; ignore_atoms=[])
    # Separate the cases unhandled by Chemfiles from the others
    path = expanduser(_path)
    cell = Cell()
    frame::Frame = if lowercase(last(splitext(path))) == ".cif"
        framecif = Frame()
        cif = expand_symmetry(CIF(path))
        a, b, c, α, β, γ = cell_parameters(cif.cell)
        set_cell!(framecif, UnitCell(a, b, c, α, β, γ))
        n = length(cif.ids)
        ignored = 0
        for i in 1:n
            typ = cif.types[cif.ids[i]]
            if typ ∈ ignore_atoms
                ignored += 1
                continue
            end
            pos = cif.cell.mat * cif.pos[:,i]
            atom = parse_atom(string(typ))
            add_atom!(framecif, atom, Vector{Float64}(pos))
        end
        n -= ignored
        if iszero(cif.bonds)
            try_guess_bonds!(framecif)
        else
            for i in 1:n, j in (i+1):n
                if cif.bonds[i,j]
                    add_bond!(framecif, i-1, j-1)
                end
            end
        end
        attributions = Int[]

        framecif
    else # The rest is handled by Chemfiles

        frameelse = read(Trajectory(path))
        for i in Int(size(frameelse))-1:-1:0
            if Symbol(type(Chemfiles.Atom(frameelse, i))) ∈ ignore_atoms
                Chemfiles.remove_atom!(frameelse, i)
            end
        end

        toremove = reverse(check_collision(positions(frameelse), matrix(UnitCell(frame))))
        for j in toremove
            Chemfiles.remove_atom!(frameelse, j-1)
        end

        topology = Topology(frameelse)
        n = Int(size(topology))

        if bonds_count(topology) == 0
            try_guess_bonds!(frameelse)
            topology = Topology(frameelse) # safer but useless since the underlying pointer is the same
        end

        attributions = attribute_residues(topology, assert_use_existing_residues)

        frameelse
    end

    cell = Cell(Cell(), SMatrix{3,3,BigFloat}(matrix(UnitCell(frame)))')
    if !all(isfinite, cell.mat) || iszero(det(cell.mat))
        @ifwarn @warn "Suspicious unit cell of matrix $(Float64.(cell.mat)). Is the input really periodic? Using a cubic unit cell instead"
        cell = Cell()
    end

    types = [Symbol(type(Chemfiles.Atom(frame, i))) for i in 0:(n-1)]
    poss = copy(positions(frame))

    adjacency = zeros(Bool, n, n)
    for (a,b) in eachcol(bonds(Topology(frame)))
        adjacency[a+1,b+1] = true
        adjacency[b+1,a+1] = true
    end
    mat = Float64.(cell.mat)
    graph = PeriodicGraph3D(n, edges_from_bonds(adjacency, mat, poss))

    if bondingmode != InputBonds
        bad_valence = check_valence(types, graph)
        recompute_bonds = bad_valence || sanity_checks!(poss, types, graph, mat, bondingmode)
        if recompute_bonds
            @ifwarn begin
                @warn "Disregarding all bonds from the input file."
                @info "To force retaining the initial bonds, use --bond-detect=input"
            end
            try_guess_bonds!(frame)
            topology = Topology(frame)
            adjacency = zeros(Bool, n, n)
            for (a,b) in eachcol(bonds(topology))
                adjacency[a+1,b+1] = true
                adjacency[b+1,a+1] = true
            end
            graph = PeriodicGraph3D(n, edges_from_bonds(adjacency, cell.mat, poss))
            @ifwarn if check_valence(types, graph)
                @warn "Remaining atoms with invalid valence. Proceeding anyway."
            end
            recompute_bonds = sanity_checks!(poss, types, graph, cell.mat, bondingmode)
            @ifwarn if recompute_bonds
                 @warn "Remaining bonds of suspicious lengths. Proceeding anyway."
            end
            attributions = attribute_residues(topology, assert_use_existing_residues)
        end
    end

    if isempty(attributions)
        crystalnothing = Crystal{Nothing}(cell, types, nothing, poss, graph)
        ifexport(crystalnothing, splitext(splitdir(path)[2])[1], exportto)
        return crystalnothing
    else
        crystalclusters = Crystal{Clusters}(cell, types, regroup_sbus(graph, attributions), poss, graph)
        ifexport(crystalclusters, splitext(splitdir(path)[2])[1], exportto)
        return crystalclusters
    end
end


text:
const C = Cairo

export applytransform!

function fill!(cc, f::Fill)
    fill!(cc, f.val)
end

function fill!(cc, g::Gradient)
    C.save(cc)
    pat = C.pattern_create_linear(g.from.xy..., g.to.xy...)
    for (stop, col) in zip(g.stops, g.colors)
        C.pattern_add_color_stop_rgba(pat, stop, rgba(col)...)
    end
    C.set_source(cc, pat)
    C.fill_preserve(cc)
    C.destroy(pat)
    C.restore(cc)
end

function fill!(cc, rg::RadialGradient)
    C.save(cc)
    pat = C.pattern_create_radial(rg.from.center.xy..., rg.from.radius, rg.to.center.xy..., rg.to.radius)
    for (stop, col) in zip(rg.stops, rg.colors)
        C.pattern_add_color_stop_rgba(pat, stop, rgba(col)...)
    end
    C.set_source(cc, pat)
    C.fill_preserve(cc)
    C.destroy(pat)
    C.restore(cc)
end

function fill!(cc, c::Colors.Colorant)
    C.save(cc)
    C.set_source_rgba(cc, rgba(c)...)
    C.fill_preserve(cc)
    C.restore(cc)
end

function fill!(cc, n::Nothing)
    # do nothing
end

function stroke!(cc, canvasmatrix, s::Stroke)
    stroke!(cc, canvasmatrix, s.val)
end

function stroke!(cc, canvasmatrix, c::Colors.Colorant)
    C.save(cc)
    C.set_source_rgba(cc, rgba(c)...)
    C.set_matrix(cc, canvasmatrix)
    C.stroke_transformed_preserve(cc)
    C.restore(cc)
end

function stroke!(cc, canvasmatrix, c::Nothing)
    # do nothing
end

function stroke!(cc, canvasmatrix, g::Gradient)
    C.save(cc)
    pat = C.pattern_create_linear(g.from.xy..., g.to.xy...)
    for (stop, col) in zip(g.stops, g.colors)
        C.pattern_add_color_stop_rgba(pat, stop, rgba(col)...)
    end
    C.set_source(cc, pat)
    C.set_matrix(cc, canvasmatrix)
    C.stroke_transformed_preserve(cc)
    C.destroy(pat)
    C.restore(cc)
end

function stroke!(cc, canvasmatrix, rg::RadialGradient)
    C.save(cc)
    pat = C.pattern_create_radial(rg.from.center.xy..., rg.from.radius, rg.to.center.xy..., rg.to.radius)
    for (stop, col) in zip(rg.stops, rg.colors)
        C.pattern_add_color_stop_rgba(pat, stop, rgba(col)...)
    end
    C.set_source(cc, pat)
    C.set_matrix(cc, canvasmatrix)
    C.stroke_transformed_preserve(cc)
    C.destroy(pat)
    C.restore(cc)
end

function clearpath!(cc)
    Cairo.new_path(cc)
end

function lineattrs!(cc, a::Attributes)
    style = a[Linestyle].val
    lw = a[Linewidth].val
    if style == :solid
        #
    elseif style == :dashed
        C.set_dash(cc, [4.0, 2] .* lw, 0)
    elseif style == :dotted
        C.set_dash(cc, [1.0, 1] .* lw, 0)
    else
        error("$style not implemented")
    end
    C.set_line_width(cc, lw)
end

function rgba(c::Colors.Colorant)
    rgba = Colors.RGBA(c)
    Float64.((rgba.r, rgba.g, rgba.b, rgba.alpha))
end

function draw_svg(canvas::Canvas)
    bufferdata = UInt8[]
    iobuffer = IOBuffer(bufferdata, read=true, write=true)
    _draw_svg(canvas, iobuffer)
    iobuffer
end

function _prep_context(canvas, surface)
    context = C.CairoContext(surface);

    C.rectangle(context, 0, 0, canvas.size...)
    C.set_source_rgba(context, rgba(canvas.color)...)
    C.fill(context)

    C.translate(context, (canvas.size ./ 2)...)
    # flip all y-coordinates
    C.scale(context, 1, -1)

    canvasmatrix = C.get_matrix(context)
    context, canvasmatrix
end

function _draw_svg(canvas, target)
    surface = C.CairoSVGSurface(target, canvas.size...);
    context, canvasmatrix = _prep_context(canvas, surface)
    draw!(context, canvasmatrix, canvas.toplayer)
    Cairo.finish(surface)
    Cairo.destroy(surface)
    nothing
end

function draw_svg(canvas::Canvas, filename::String)
    _draw_svg(canvas, filename)
end

function draw_pdf(canvas::Canvas, filename::String)

    c = C.CairoPDFSurface(filename, canvas.size...);
    cc = C.CairoContext(c);

    C.rectangle(cc, 0, 0, canvas.size...)
    C.set_source_rgba(cc, rgba(canvas.color)...)
    C.fill(cc)

    C.translate(cc, (canvas.size./2)...)

    canvasmatrix = C.get_matrix(cc)

    draw!(cc, canvasmatrix, canvas.toplayer)
    c
end

function draw_rgba(canvas::Canvas; px_per_pt=1/0.75)

    size_pixel = round.(Int, canvas.size .* px_per_pt)

    c = C.CairoARGBSurface(size_pixel...);
    cc = C.CairoContext(c);

    begin
        font_options_ptr = ccall((:cairo_font_options_create, C.libcairo), Ptr{Nothing}, ())

        CAIRO_HINT_STYLE_NONE = 1
        CAIRO_HINT_STYLE_SLIGHT = 2
        CAIRO_HINT_STYLE_FULL = 4
        ccall(
            (:cairo_font_options_set_hint_style, C.libcairo), Nothing,
            (Ptr{Nothing}, Int32),
            font_options_ptr, CAIRO_HINT_STYLE_NONE)

        CAIRO_HINT_METRICS_DEFAULT = 0
        CAIRO_HINT_METRICS_OFF = 1
        CAIRO_HINT_METRICS_ON = 2

        ccall(
            (:cairo_font_options_set_hint_metrics, C.libcairo), Nothing,
            (Ptr{Nothing}, Int32),
            font_options_ptr, CAIRO_HINT_METRICS_OFF)

        # ccall(
        #     (:cairo_font_options_set_antialias, C.libcairo), Nothing,
        #     (Ptr{Nothing}, Int32),
        #     font_options_ptr, 3)

        ccall(
            (:cairo_set_font_options, C.libcairo), Nothing,
            (Ptr{Nothing}, Ptr{Nothing}),
            cc.ptr, font_options_ptr)

        ccall(
            (:cairo_font_options_destroy, C.libcairo), Nothing,
            (Ptr{Nothing},),
            font_options_ptr)
    end

    C.rectangle(cc, 0, 0, size_pixel...)
    C.set_source_rgba(cc, rgba(canvas.color)...)
    C.fill(cc)

    C.translate(cc, (size_pixel./2)...)
    C.scale(cc, px_per_pt, px_per_pt)
    # flip all y-coordinates
    C.scale(cc, 1, -1)

    canvasmatrix = C.get_matrix(cc)

    draw!(cc, canvasmatrix, canvas.toplayer)
    c
end

function applytransform!(cc, t::Transform)
    C.translate(cc, t.translation...)
    C.scale(cc, t.scale, t.scale)
    C.rotate(cc, rad(t.rotation))
end

function Base.convert(::Type{Int32}, op::Operator)
    @match op.operator begin
        :over => C.OPERATOR_OVER
        :add => C.OPERATOR_ADD
        :mult => C.OPERATOR_MULTIPLY
        :screen => C.OPERATOR_SCREEN
        :sat => C.OPERATOR_SATURATE
        :darken => C.OPERATOR_DARKEN
        :overlay => C.OPERATOR_OVERLAY
    end
end

function draw!(cc::C.CairoContext, canvasmatrix, l::Layer)

    C.save(cc)
    C.push_group(cc)
    applytransform!(cc, gettransform!(l, cc))

    for content in l.content
        draw!(cc, canvasmatrix, content)
    end

    gr = C.pop_group(cc)
    C.set_source(cc, gr);
    C.destroy(gr)

    # I guess this should go here so it only takes effect on the layer after
    # that is drawn, so there will be better fringes than if two antialiased clips
    # are applied on perfectly overlapping objects
    setclippath!(cc, l)

    C.set_operator(cc, Base.convert(Int32, l.operator))

    C.paint_with_alpha(cc, l.opacity.opacity)
    C.restore(cc)
end

function getattributes(s::Shape{T}) where T
    needed = needed_attributes(T)
    attributes = Attributes(Dict(attr => getattribute(s, attr) for attr in needed))
end

function getattributes(s::Shapes{T}) where T
    needed = needed_attributes(T)
    attributes = Attributes(Dict(attr => getattribute(s, attr) for attr in needed))
end

function getattribute(l::LayerContent, attr)
    if haskey(l.attrs, attr)
        return l.attrs[attr]
    else
        if isnothing(l.parent)
            error("Attribute $(attr) could not be found in the layer chain.")
        else
            return getattribute(l.parent, attr)
        end
    end
end

function draw!(cc, canvasmatrix, s::Shape)
    geom = solve!(s, cc)
    attributes = getattributes(s)

    #fast exit for invisible shapes
    if !attributes[Visible].val
        return
    end

    C.save(cc)

    C.push_group(cc)
    draw!(cc, canvasmatrix, geom, attributes)
    gr = C.pop_group(cc)
    C.set_source(cc, gr)
    C.destroy(gr)

    setclippath!(cc, s)

    C.set_operator(cc, Base.convert(Int32, s.operator))

    C.paint_with_alpha(cc, s.opacity.opacity)

    C.restore(cc)
end

function getattributes(attr::Attributes, i::Int, n::Int, geom::GeometricObject)
    singleattrs = Attributes()
    for (typ, a) in attr.attrs
        if ismultiattr(a)
            if typeof(a.val) <: Pair{Symbol, <:Function}
                sym = a.val[1]
                func = a.val[2]
                value = @match sym begin
                    :i => func(i)
                    :frac => func((i - 1) / (n - 1))
                    :geom => func(geom)
                    _ => error("Not implemented")
                end
                insert!(singleattrs, typ(value))
            elseif typeof(a.val) <: AbstractArray
                insert!(singleattrs, typ(a.val[i]))
            else
                error("Not implemented")
            end
        else
            insert!(singleattrs, a)
        end
    end
    singleattrs
end

function draw!(cc, canvasmatrix, s::Shapes)
    geoms = solve!(s, cc)
    attributes = getattributes(s)

    #fast exit for invisible shapes
    if attributes[Visible].val == false
        return
    end

    C.save(cc)
    C.push_group(cc)

    n = length(geoms)
    for (i, g) in enumerate(geoms)
        singleattrs = getattributes(attributes, i, n, g)
        draw!(cc, canvasmatrix, g, singleattrs)
    end

    gr = C.pop_group(cc)
    C.set_source(cc, gr);
    C.destroy(gr)

    setclippath!(cc, s)

    C.set_operator(cc, Base.convert(Int32, s.operator))

    C.paint_with_alpha(cc, s.opacity.opacity)

    C.restore(cc)
end

needslineattrs(g) = false
needslineattrs(g::Union{Point, Line, Bezier, Rect, Circle, Polygon, Path}) = true
needsstroke(g) = false
needsstroke(g::Union{Point, Line, Bezier, Rect, Circle, Polygon, Path}) = true
needsfill(g) = false
needsfill(g::Union{Point, Bezier, Rect, Circle, Polygon, Path}) = true

function draw!(cc, canvasmatrix, g::GeometricObject, a::Attributes)
    if typeof(g) <: Point
        makepath!(cc, g, canvasmatrix, a[Marker], a[Markersize])
    else
        makepath!(cc, g)
    end

    needslineattrs(g) && lineattrs!(cc, a)
    needsfill(g) && fill!(cc, a[Fill])
    needsstroke(g) && stroke!(cc, canvasmatrix, a[Stroke])

    clearpath!(cc)
end

# function draw!(cc, canvasmatrix, p::Point, a::Attributes)
#     C.move_to(cc, (p + 0.5 * a[Markersize].size * X(1)).xy...)
#     C.arc(cc, p.x, p.y, 0.5 * a[Markersize].size, 0, 2pi)
#     fill!(cc, canvasmatrix, a[Fill])
#     stroke!(cc, canvasmatrix, a[Stroke])
#     clearpath!(cc)
# end

# function draw(ps::Points, a::Attributes)
#     PyPlot.scatter(
#         xs(ps.points), ys(ps.points),
#         s = a[Markersizes].sizes,
#         color = rgba.(a[Strokes].colors),
#         marker = a[Marker].marker,
#     )
# end

function makepath!(cc, p::Point, canvasmatrix, m::Marker, ms::Markersize)
    m = m.val
    s = ms.val

    # save the position of the point in device coordinates
    xdev, ydev = C.user_to_device!(cc, [p.xy...])
    C.save(cc)
    # set the canvasmatrix and convert the device space point to this space
    C.set_matrix(cc, canvasmatrix)
    xuser, yuser = C.device_to_user!(cc, [xdev, ydev])
    # translate onto the point
    C.translate(cc, xuser, yuser)

    # now the shapes can be defined with origin (0, 0) and size in canvas units
    shape = @match m begin
        :square => Rect(O, s, s, deg(0))
        :circle || :o => Circle(O, 0.5s)
        :cross || :+ => Polygon(ncross(O, 4, 0.5s, 1 / Base.MathConstants.golden))
        _ => error("Markertype $m is not implemented")
    end

    makepath!(cc, shape)
    C.restore(cc)
end

function makepath!(cc, l::Line)
    C.move_to(cc, l.from.xy...)
    C.line_to(cc, l.to.xy...)
end

function continuepath!(cc, l::Line)
    C.line_to(cc, l.to.xy...)
end

# function draw(ls::LineSegments, a::Attributes)
#     path = PyPlot.matplotlib.path.Path(ls, closed=false)
#     pathpatch = PyPlot.matplotlib.patches.PathPatch(
#         path,
#         edgecolor = rgba(a[Stroke].color),
#         linewidth = a[Linewidth].width,
#         lineattrs = a[Linestyle].style,
#         zorder=zorder(),
#     )
#     PyPlot.gca().add_patch(pathpatch)
# end

function makepath!(cc, p::Polygon)
    C.move_to(cc, p.points[1].xy...)

    for po in p.points[2:end]
        C.line_to(cc, po.xy...)
    end

    C.close_path(cc)
end

function makepath!(cc, a::Arc)
    pstart = fraction(a, 0)
    C.move_to(cc, pstart.xy...)

    if a.end_angle - a.start_angle >= rad(0)
        C.arc(cc, a.center.xy..., a.radius, rad(a.start_angle), rad(a.end_angle))
    else
        C.arc_negative(cc, a.center.xy..., a.radius, rad(a.start_angle), rad(a.end_angle))
    end
end

function continuepath!(cc, a::Arc)
    if a.end_angle - a.start_angle >= rad(0)
        C.arc(cc, a.center.xy..., a.radius, rad(a.start_angle), rad(a.end_angle))
    else
        C.arc_negative(cc, a.center.xy..., a.radius, rad(a.start_angle), rad(a.end_angle))
    end
end

alignments = Dict(
    :c => "center",
    :t => "top",
    :b => "bottom",
    :l => "left",
    :r => "right",
    :bl => "baseline",
    :cbl => "center_baseline"
)

function TextExtent(cc, t::Txt)
    C.save(cc)
    C.select_font_face(cc, t.font, Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_NORMAL)
    C.set_font_size(cc, t.size)
    e = C.text_extents(cc, t.text);
    C.restore(cc)
    TextExtent(e...)
    # typedef struct {
    #     double x_bearing;
    #     double y_bearing;
    #     double width;
    #     double height;
    #     double x_advance;
    #     double y_advance;
    # } cairo_text_extents_t;
end

function draw!(cc, canvasmatrix, t::Txt, a::Attributes)

    # if typeof(a[Fill]) <: Gradient
    #     g = a[Fill]
    #     pat = C.pattern_create_linear(g.from.xy..., g.to.xy...);
    #     for (stop, col) in zip(g.stops, g.colors)
    #         C.pattern_add_color_stop_rgba(pat, stop, rgba(col)...);
    #     end
    #     C.set_source(cc, pat);
    #     C.fill_preserve(cc);
    #     C.destroy(pat);
    #     end
    # end
    ex = TextExtent(cc, t)

    C.set_source_rgba(cc, rgba(a[Textfill].val)...)
    C.set_font_face(cc, "$(t.font) $(t.size)")

    C.set_text(cc, t.text, false)
    C.update_layout(cc)

    function get_layout_size(ctx::C.CairoContext)
        w = Vector{Int32}(undef, 2)
        ccall((:pango_layout_get_size, C.libpango), Nothing,
              (Ptr{Nothing},Ptr{Int32},Ptr{Int32}), ctx.layout, pointer(w,1), pointer(w,2))
        w ./ 1024
    end

    w, h = get_layout_size(cc)
    # C.set_font_size(cc, t.size)

    shiftx = @match t.halign begin
        :l => 0
        :c => -w/2
        :r => -w
        _ => error("Alignment $(t.halign) doesn't exist")
    end

    shifty = @match t.valign begin
        :t => 0
        :c => -h/2
        :b => -h
        _ => error("Alignment $(t.valign) doesn't exist")
        # :bl => 0
        # :b => -(ex.ybearing + ex.height)
    end

    shift = rotate(P(shiftx, shifty), -t.angle)
    shiftflipped = P(shift.x, -shift.y)

    pos = t.pos + shiftflipped

    C.save(cc) # this is needed if multiple texts are drawn in one go

    C.move_to(cc, pos.xy...)

    C.scale(cc, 1, -1)

    C.rotate(cc, rad(-t.angle))

    C.show_layout(cc)
    C.restore(cc)
end

function makepath!(cc, b::Bezier)
    C.move_to(cc, b.from.xy...)
    C.curve_to(cc, b.c1.xy..., b.c2.xy..., b.to.xy...)
end

function continuepath!(cc, b::Bezier)
    C.curve_to(cc, b.c1.xy..., b.c2.xy..., b.to.xy...)
end

start(a::Arc) = fraction(a, 0)
stop(a::Arc) = fraction(a, 1)
start(b::Bezier) = b.from
stop(b::Bezier) = b.to
start(l::Line) = l.from
stop(l::Line) = l.to

function makepath!(cc, m::Move)
    C.move_to(cc, m.p.xy...)
end

function makepath!(cc, l::Lineto)
    C.line_to(cc, l.p.xy...)
end

function makepath!(cc, c::CurveTo)
    C.curve_to(cc, c.c1.xy..., c.c2.xy..., c.p.xy...)
end

function makepath!(cc, c::Close)
    C.close_path(cc)
end

function makepath!(cc, p::Path)
    makepath!.(Ref(cc), p.commands)
    # makepath!(cc, p.segments[1])
    # for i in 2:length(p.segments)
    #     news = p.segments[i]
    #     olds = p.segments[i-1]
    #     if isapprox(start(news), stop(olds))
    #         continuepath!(cc, news)
    #     else
    #         makepath!(cc, news)
    #     end
    # end
    # if p.closed
    #     C.close_path(cc)
    # end
end

# function draw(bs::Paths, a::Attributes)
#     paths = PyPlot.matplotlib.path.Path.(bs.paths)
#     collection = PyPlot.matplotlib.collections.PathCollection(
#         paths,
#         edgecolors = rgbas(a[Strokes]),
#         facecolors = rgbas(a[Fills]),
#         linewidths = a[Linewidths].widths,
#         zorder=zorder(),
#         antialiaseds=true,
#         snap=false,
#     )
#     PyPlot.gca().add_collection(collection)
# end

# function PyPlot.plot(v::Vector{Point}; kwargs...)
#     xx = [p.x for p in v]
#     yy = [p.y for p in v]
#     PyPlot.plot(xx, yy; kwargs...)
# end

function makepath!(cc, c::Circle)
    # C.move_to(cc, c.center.xy...)
    C.new_sub_path(cc)
    C.arc(cc, c.center.xy..., c.radius, 0, 2pi)
end

function setclippath!(cc, l::LayerContent)

    sh = l.clip.shape
    if isnothing(sh)
        return
    end
    if typeof(sh) <: Shape
        s = transform_from_to(sh, l, cc) * solve!(sh, cc)
        makepath!(cc, s)
    else
        s_inner = transform_from_to(sh[1], l, cc) * solve!(sh[1], cc)
        s_outer = transform_from_to(sh[2], l, cc) * solve!(sh[2], cc)
        makepath!(cc, s_inner)
        makepath!(cc, s_outer)
    end
    C.clip(cc)
end

function makepath!(cc, r::Rect)
    C.move_to(cc, bottomleft(r).xy...)
    C.line_to(cc, bottomright(r).xy...)
    C.line_to(cc, topright(r).xy...)
    C.line_to(cc, topleft(r).xy...)
    C.close_path(cc)
end


text:
# ------------------------------------------------------------------------
# Protecting access via property interface

mutable struct FileContent
    path
    loaded
    contents
end

function FileContent(path)
    ss = lstat(path)
    return FileContent(path, false, zeros(UInt8, ss.size))
end

function load_contents!(fc::FileContent)
    open(fc.path) do io
        readbytes!(io, getfield(fc, :contents)) # not f.contents
        fc.loaded = true
    end
    nothing
end

# quick test

fc = FileContent("/etc/hosts")
fc.loaded
load_contents!(fc)
fc.loaded
fc.contents
#=
julia> fc = FileContent("/etc/hosts")
FileContent("/etc/hosts", false, UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00  …  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00])

julia> fc.loaded
false

julia> load_contents!(fc)

julia> fc.loaded
true

julia> fc.contents
359-element Array{UInt8,1}:
 0x23
 0x23
 0x0a
=#

# how does the compiler work with the dot notation?

# reference

Meta.lower(Main, :( fc.path ))
#=
julia> Meta.lower(Main, :( fc.path ))
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope'
1 ─ %1 = Base.getproperty(fc, :path)
└──      return %1
))))
=#

# assignment

Meta.lower(Main, :( fc.path = "/etc/hosts"))
#=
julia> Meta.lower(Main, :( fc.path = "/etc/hosts"))
:($(Expr(:thunk, CodeInfo(
    @ none within `top-level scope'
1 ─     Base.setproperty!(fc, :path, "/etc/hosts")
└──     return "/etc/hosts"
))))
=#

# take a look at Julia's source code :-)
#=
julia> @edit fc.path
=#

# got this:
#=
getproperty(Core.@nospecialize(x), f::Symbol) = getfield(x, f)
setproperty!(x, f::Symbol, v) = setfield!(x, f, convert(fieldtype(typeof(x), f), v))
=#

# extend Base.getproperty to redefine the meaning of dot notation
function Base.getproperty(fc::FileContent, s::Symbol)
    direct_passthrough_fields = (:path, )
    if s in direct_passthrough_fields
        return getfield(fc, s)
    end
    if s === :contents
        getfield(fc, :loaded) || load_contents!(fc)
        return getfield(fc, :contents)
    end
    error("Unsupported property: $s")
end

# lazy load
function load_contents!(fc::FileContent)
    open(getfield(fc, :path)) do io
        readbytes!(io, getfield(fc, :contents)) # not f.contents
        setfield!(fc, :loaded, true)
    end
    nothing
end

fc = FileContent("/etc/hosts")

# test

fc = FileContent("/etc/hosts")
fc.path
fc.contents
fc.loaded
#=
julia> fc = FileContent("/etc/hosts")
FileContent("/etc/hosts", false, UInt8[0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00  …  0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00])

julia> fc.path
"/etc/hosts"

julia> fc.contents
359-element Array{UInt8,1}:
 0x23
 0x23
 0x0a
    ⋮

julia> fc.loaded
ERROR: Unsupported property: loaded
=#

function Base.setproperty!(fc::FileContent, s::Symbol, value)
    if s === :path
        ss = lstat(value)
        setfield!(fc, :path, value)
        setfield!(fc, :loaded, false)
        setfield!(fc, :contents, zeros(UInt8, ss.size))
        println("Object re-initialized for $value (size $(ss.size))")
        return nothing
    end
    error("Property $s cannot be changed.")
end

fc.path
fc.path = "/etc/profile"
fc.contents = []
#=
julia> fc.path = "/etc/profile"
Object re-initialized for /etc/profile (size 189)
"/etc/profile"

julia> fc.contents = []
ERROR: Property contents cannot be changed.
=#

# before propertynames
#=
julia> fc.
contents loaded    path
=#

# expose property names to tools like the REPL
function Base.propertynames(fc::FileContent)
    return (:path, :contents)
end

#=
julia> fc.
contents path
=#


text:
using CSV
using Statistics
using StatsBase
using Random
using Distributions
using LinearAlgebra
using SparseArrays
using Bootstrap
using ArgParse
import Bootstrap:draw!
import Base:copy, length

using JLD

BLAS.set_num_threads(1)

function make_strats(num_strats)
    K = length(num_strats)
    strats = [1:n for n in num_strats]
    idx2strat = Dict{Int64, NTuple{K,Int64}}([n => s for (n, s) in enumerate(Iterators.product(strats...))]...)
    strat2idx = Dict{NTuple{K,Int64}, Int64}([s => n for (n, s) in enumerate(Iterators.product(strats...))]...)
    strats_tuples = [s  for s in Iterators.product(strats...)]
    alt_strats = []
    for k in 1:K
        palts = []
        for sidx in 1:length(strats_tuples)
            psalts = []
            s = idx2strat[sidx]
            for i in 1:num_strats[k]
                ns = ntuple(j->(j==k ? i : s[j]), length(s))
                push!(psalts, strat2idx[ns])
            end
            push!(palts, psalts)
        end
        push!(alt_strats, palts)
    end

    return strats_tuples, idx2strat, strat2idx, alt_strats
end

function compute_eta(num_strats)
    return 1.0 / sum(num_strats .- 1)
end

function makeCinfinite_dense!(C, strats, alts, players, payoffs, η, m)
    numS = length(strats)
    j=0
    f1 = 0.0
    f2 = 0.0
    tie = η / m
    for i in 1:numS
        for k in 1:players
            for j in alts[k][i]
                if j == i
                    continue
                else
                    f1 = payoffs(k, j)
                    f2 = payoffs(k, i)
                    if f1 > f2
                        C[i,j] = η
                    #elseif f2 > f1
                    #    C[i,j] = 0.0
                    elseif isapprox(f1,f2)
                        C[i,j] = tie
                    #else
                    #    throw(Error("unknown condition: f1: $(f1) f2: $(f2)"))
                    end
                end
            end
        end
    end
    for s in 1:numS
        C[s,s] = 1.0 - sum(view(C,s,:))
    end
    return C
end

function makeCinfinite_bounds(strats, alts, players, payoffs, η, m)
    Ic = Array{Int64,1}();
    Jc = Array{Int64,1}();
    Vl = Array{Float64,1}();
    Vu = Array{Float64,1}();
    num_uncertain = 0
    Eu = Array{Array{Int64,1},1}()
    numS = length(strats)
    tie = η/m
    for i in 1:numS
        uncertain = Array{Int64,1}()
        push!(uncertain, i)
        for k in 1:players
            for j in alts[k][i]
                if j == i
                    continue
                else
                    f1_low  = payoffs(k, j, :left)
                    f1_high = payoffs(k, j, :right)
                    f2_low  = payoffs(k, i, :left)
                    f2_high = payoffs(k, i, :right)
                    push!(Ic, i)
                    push!(Jc, j)
                    if (f1_low > f2_high)
                        push!(Vl, η)
                        push!(Vu, η)
                    elseif f2_low > f1_high
                        push!(Vl, 0)
                        push!(Vu, 0)
                    elseif isapprox(f1_low,f2_low) && isapprox(f1_high, f2_high)
                        push!(Vl, tie)
                        push!(Vu, tie)
                    else
                        push!(Vl, 0.0)
                        push!(Vu, η)
                        push!(uncertain, j)
                        num_uncertain += 1
                    end
                end
            end
        end
        push!(Eu, uncertain)
    end
    numS = length(strats)
    Clow = sparse(Ic, Jc, Vl, numS, numS)
    Chigh = sparse(Ic, Jc, Vu, numS, numS)
    dropzeros!(Clow)
    dropzeros!(Chigh)
    for s in 1:numS
        ch = 1.0 - sum(Clow[s,:])
        Clow[s,s] = 1.0 - sum(Chigh[s,:])
        Chigh[s,s] = ch
    end
    return Clow, Chigh, Eu, num_uncertain
end

function valueiteration_row4sp!(Pi, Pi_low, Pi_high, Ri, v, γ, numS, p0, rpert; Eu)
    Pi[Eu] .= Pi_low[Eu]
    if length(Eu) == 0
        return nothing
    end
    w = zeros(length(Eu))

    @. w =  (Ri + γ * v[Eu] + rpert[Eu])
    idxs = sortperm(w)
    p = p0

    for idx in idxs[end:-1:1]
        j = Eu[idx]
        pjh = Pi_high[j]
        pjl = Pi_low[j]
        dp = min(pjh-pjl, 1.0 - p)
        Pi[j] = pjl+dp
        p += dp
    end
end

vierror_bound(ϵ, γ) = (2*ϵ*γ)/(1-γ)

function valueiteration_contagg(P_low, P_high, Eu, Rs, numS, direction; γ=(1.0 - 1e-8))
    R = zeros(size(Rs))
    if direction == :max
        R .= Rs
    elseif direction == :min
        R .= -Rs
    else
        println("ERROR unrecognized symbol: ", direction)
        return NaN
    end
    v = zeros(numS)
    vold = zeros(numS)
    P = zeros((numS, numS))
    Ptmps = zeros((numS, numS))
    Ptmps .= P_low
    P .= P_low
    P0 = [sum(P_low[i, :]) for i in 1:numS]
    rper = rand(Float64, numS) .* (1e-10 / (1-γ))
    iteration = 0
    changed = true
    while changed
        t0 = time_ns()
        iteration += 1
        if iteration > 400
            ϵ = norm(v .- vold, Inf)
            bound = vierror_bound(ϵ, γ) # bound on max distance to optimal value function for any state
            bound = (1.0 - γ) * bound # bound on max error to bound on alpha rank (do this to converge regradless of gamma choice)
            println("hit iteration limit, epsilon = $(ϵ), bound = $(bound)")
            break
        end
        vold .= v
        changed = false
        for s in 1:numS
            edges = Eu[s]

            valueiteration_row4sp!(view(Ptmps, s, :), view(P_low, s, :), view(P_high, s, :), R[s], vold, γ, numS, P0[s], rper, Eu=edges)
            if any(abs.(view(Ptmps, s, edges) .- view(P, s, edges)) .>= 1e-8)
                changed=true
            end
            P[s, edges] .= view(Ptmps,s, edges)
        end
        t1 = time_ns()
        v .= inv(I - γ .* P) * R
        ϵ = norm(v .- vold, Inf)
        t2 = time_ns()
        tot = (t2-t1) / 1.0e9
        bound = vierror_bound(ϵ, γ) # bound on max distance to optimal value function for any state
        bound = (1.0 - γ) * bound # bound on max error to bound on alpha rank (do this to converge regradless of gamma choice)
        if bound < 1e-7
            changed = false
        end
        tf = time_ns()
        tot2 = (tf-t0) / 1.0e9
    end
    sol = (1.0 - γ) * mean(abs.(v))
    return sol
end

function (ecdf::ECDF)(x::Real, δ, tail=:both, method=:DKW)
    isnan(x) && return NaN
    p = ecdf(x)
    n = length(ecdf.sorted_values)
    ϵ = √(log(2.0/δ)/2n)
    if tail == :both
        return (max(p-ϵ, 0.0), min(p+ϵ, 1.0))
    elseif tail == :left
        return max(p-ϵ, 0.0)
    elseif tail == :right
        return min(p+ϵ, 1.0)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

function (ecdf::ECDF)(x::Real, δ, a::Real, b::Real, tail=:both, method=:DKW)
    #TODO make epsilon depend on tail, tail=both used delta/2
    isnan(x) && return NaN
    n = length(ecdf.sorted_values)
    ϵ = √(log(2.0/δ)/2n)
    pl = 0.0
    ph = 0.0
    if x < a
        pl = 0.0
        ph = 0.0
    elseif x ≥ b
        pl = 1.0
        ph = 1.0
    else
        p = ecdf(x)
        pl = max(p-ϵ, 0.0)
        ph = min(p+ϵ, 1.0)
    end

    if tail == :both
        return pl, ph
    elseif tail == :left
        return pl
    elseif tail == :right
        return ph
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

function andersons_mean(D, δ, a, b, tail=:both)
    # assume D sorted
    # alpha is cvar threshold
    # δ is CI to hold with probability 1-δ
    # a is the lower bound on the Data
    # b is the upper bound on the Data
    N = length(D)
    ϵ = √(log(2.0/δ)/(2*N))
    if tail == :both
        Z = vcat(D, b)
        tmp = clamp.(collect(Float64, 1:N) ./ N .- ϵ, 0.0, Inf)
        upper = b - sum(diff(Z) .* tmp)

        Z = vcat(a, D)
        tmp = clamp.(collect(Float64, 0:N-1) ./ N .+ ϵ, -Inf, 1.0)
        lower = D[end] - sum(diff(Z) .* tmp)
        return lower, upper
    elseif tail == :left
        Z = vcat(a, D)
        tmp = clamp.(collect(Float64, 0:N-1) ./ N .+ ϵ, -Inf, 1.0)
        lower = D[end] - sum(diff(Z) .* tmp)
        return lower
    elseif tail == :right
        Z = vcat(D, b)
        tmp = clamp.(collect(Float64, 1:N) ./ N .- ϵ, 0.0, Inf)
        upper = b - sum(diff(Z) .* tmp)
        return upper
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

function andersons_meangx(D, g, δ, a, b, tail=:both)
    # assume D sorted
    # alpha is cvar threshold
    # δ is CI to hold with probability 1-δ
    # a is the lower bound on the Data
    # b is the upper bound on the Data
    N = length(D)
    ϵ = √(log(2.0/δ)/(2*N))
    if tail == :both
        Z = vcat(D, b)
        gZ = g.(Z)
        tmp = clamp.(collect(Float64, 1:N) ./ N .- ϵ, 0.0, Inf)
        upper = gZ[end] - sum(diff(gZ) .* tmp)

        Z .= vcat(a, D)
        gZ .= g.(Z)
        tmp = clamp.(collect(Float64, 0:N-1) ./ N .+ ϵ, -Inf, 1.0)
        lower = gZ[end] - sum(diff(gZ) .* tmp)
        return lower, upper
    elseif tail == :left
        Z = vcat(a, D)
	gZ = g.(Z)
        tmp = clamp.(collect(Float64, 0:N-1) ./ N .+ ϵ, -Inf, 1.0)
        lower = gZ[end] - sum(diff(gZ) .* tmp)
        return lower
    elseif tail == :right
        Z = vcat(D, b)
        gZ = g.(Z)
        tmp = clamp.(collect(Float64, 1:N) ./ N .- ϵ, 0.0, Inf)
        upper = gZ[end] - sum(diff(gZ) .* tmp)
        return upper
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

function cdfnormalized_mean(Fx, Fy, δx, δy, a, b, tail=:both)
    # Fx is normalizing distribution
    # Fy is empirical cdf of input variable
    # δx is CI to hold with probability 1-δx for Fx
    # δy is CI to hold with probability 1-δx for Fy
    # a is the lower bound on the distribution of X and Y
    # b is the upper bound on the distribution of X and Y
    y = Fy.sorted_values
    if tail == :both
        gl(x) = Fx(x, δx, a, b, :left, :DKW)
        lower = andersons_meangx(y, gl, δy, a, b, :left)
        gu(x) = Fx(x, δx, a, b, :right, :DKW)
        upper = andersons_meangx(y, gu, δy, a, b, :right)
        return lower, upper
    elseif tail == :left
        g(x) = Fx(x, δx, a, b, :left, :DKW)
        lower = andersons_meangx(y, g, δy, a, b, :left)
        return lower
    elseif tail == :right
        g(x) = Fx(x, δx, a, b, :right, :DKW)
        upper = andersons_meangx(y, g, δy, a, b, :right)
        return upper
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

struct PerfDists
    dists::Dict{String,Dict{String,ECDF{Array{Float64,1}}}}
    algs::Array{String,1}
    envs::Array{String,1}
    bounds::Dict{String,Tuple{Float64,Float64}}
end

function copy(x::PerfDists)
    return deepcopy(x)
end

# function size(x::PerfDists)
#     return (length(x.algs), length(x.envs), length(x.dists[x.algs[1]][x.envs[1]].sorted_values))
# end

function length(x::PerfDists)
    tot = 0
    for alg in x.algs
        for env in x.envs
            tot += length(x.dists[alg][env].sorted_values)
        end
    end
    return tot
end

function draw!(x::PerfDists, o::PerfDists)
    algs = x.algs
    envs = x.envs
    for alg in algs
        for env in envs
            d = x.dists[alg][env].sorted_values
            idx = sample(1:length(d), length(d))
            o.dists[alg][env] = ecdf(d[idx])
        end
    end
end

function make_perfs_nobound(D::PerfDists)
    n,m = length(D.algs), length(D.envs)
    A = zeros(Float64, (n,m,n))
    N = length(D.dists[D.algs[1]][D.envs[1]].sorted_values)
    pts = zeros(N)
    for (i,alg) in enumerate(D.algs)
        for (j,env) in enumerate(D.envs)
            for (k,algnorm) in enumerate(D.algs)
                samedef = false
                if occursin("gw", env) || occursin("chain", env)
                    if (occursin("ac", alg) && occursin("ac", algnorm)) || (occursin("sarsa", alg) && occursin("sarsa", algnorm)) || (occursin("qlambda", alg) && occursin("qlambda", algnorm))
                        if (occursin("scaled", alg) && occursin("normal", algnorm)) || (occursin("normal", alg) && occursin("scaled", algnorm))
                            samedef = true
                        end
                    end
                end
                if i==k || samedef
                    A[i,j,k] = 0.5
                else
                    x = D.dists[alg][env].sorted_values
                    pts .= D.dists[algnorm][env].(x)
                    p = mean(pts)
                    A[i,j,k] = p
                end
            end
        end
    end
    return A
end

function make_perfs_bound(D::PerfDists, δ, bound=:DKW)
    n,m = length(D.algs), length(D.envs)
    δ′ = δ/(n*m)
    A = zeros(Float64, (n,m,n))
    Alow = zeros(Float64, (n,m,n))
    Ahigh = zeros(Float64, (n,m,n))
    for (i,alg) in enumerate(D.algs)
        for (j,env) in enumerate(D.envs)
            a,b = D.bounds[env]
            for (k,algnorm) in enumerate(D.algs)
                samedef = false
                if occursin("gw", env) || occursin("chain", env)
                    if (occursin("ac", alg) && occursin("ac", algnorm)) || (occursin("sarsa", alg) && occursin("sarsa", algnorm)) || (occursin("qlambda", alg) && occursin("qlambda", algnorm))
                        if (occursin("scaled", alg) && occursin("normal", algnorm)) || (occursin("normal", alg) && occursin("scaled", algnorm))
                            samedef = true
                        end
                    end
                end
                if i==k || samedef
                    A[i,j,k] = 0.5
                    Alow[i,j,k] = 0.5
                    Ahigh[i,j,k] = 0.5
                else
                    x = D.dists[alg][env].sorted_values
                    N = length(x)
                    pts = zeros(N)
                    pts .= D.dists[algnorm][env].(x)
                    p = mean(pts)
                    A[i,j,k] = p

                    if bound == :AndersonDKW
                        pts .= D.dists[algnorm][env].(x, δ′, a, b, :left, :DKW)
                        pl = andersons_mean(pts, δ′, 0, 1, :left)
                        Alow[i,j,k] = pl
                        pts .= D.dists[algnorm][env].(x, δ′, a, b, :right, :DKW)
                        pu = andersons_mean(pts, δ′, 0, 1, :right)
                        Ahigh[i,j,k] = pu
                    elseif bound == :DKW
                        pl, pu = cdfnormalized_mean(D.dists[algnorm][env], D.dists[alg][env], δ′, δ′, a, b, :both)
                        Alow[i,j,k] = pl
                        Ahigh[i,j,k] = pu
                    elseif bound == :TTest
                        δ2 = δ / (n*n*m)
                        tstar = quantile(TDist(N-1), 1 - δ2/2)
                        μ, σ = mean(pts), std(pts)
                        pl = μ - (σ / √N) * tstar
                        pu = μ + (σ / √N) * tstar
                        Alow[i,j,k] = pl
                        Ahigh[i,j,k] = pu
                    elseif bound == :TTest2
                        tstar = quantile(TDist(N-1), 1 -  δ′/2)
                        μ, σ = mean(pts), std(pts)
                        pl = μ - (σ / √N) * tstar
                        pu = μ + (σ / √N) * tstar
                        Alow[i,j,k] = pl
                        Ahigh[i,j,k] = pu
                    end
                end
            end
        end
    end
    return A, Alow, Ahigh
end

function load_results(paths, envs, algs)
    D = Dict{String,Dict{String,ECDF{Array{Float64,1}}}}()
    ab = Dict{String,Tuple{Float64,Float64}}()
    for env in envs
        ab[env] = (Inf, -Inf)
    end
    for alg in algs
        algD = Dict{String,ECDF{Array{Float64,1}}}()
        D[alg] = algD
    end
    for (path, env) in zip(paths, envs)
        df = CSV.read(path)
        for alg in algs
            D[alg][env] = ecdf(df[:, Symbol(alg)])
            a,b = D[alg][env].sorted_values[[1, end]]
            a = min(ab[env][1], a)
            b = max(ab[env][2], b)
            ab[env] = (a,b)
        end
    end
    perfs = PerfDists(D, algs, envs, ab)
    return perfs
end

function compute_aggregate(perfs::PerfDists, moves, strats, idx2strat, strat2idx, alts)
    A = make_perfs_nobound(perfs)
    numA = length(perfs.algs)
    numE = length(perfs.envs)
    aidx(i) = Int((i-1) ÷ numE)+1
    eidx(i) = Int((i-1) % numE)+1

    function payoffs(n,s::Int,tail=:none)
        if n == 1
            return A[s]
        else
            return -A[s]
        end
    end

    function payoffs(n,s::Tuple,tail=:none)
        m1, m2 = s
        j,k = eidx(m2), aidx(m2)

        if n == 1
            return A[m1,j,k]
        else
            return -A[m1,j,k]
        end
    end

    players = length(moves)

    η = compute_eta(moves);
    min_p = 0.0
    m = 50
    numS = length(strats)
    C = zeros((numS, numS))
    makeCinfinite_dense!(C, strats, alts, players, payoffs, η, m);
    R = zeros(numS)
    γ = (numS - 1.0) / Float64(numS)

    y = zeros(Float64, numA)
    invC = inv(I - γ .* C)
    for i in 1:numA
        for s in 1:numS
            R[s] = A[i,eidx(idx2strat[s][2]),aidx(idx2strat[s][2]),1]
        end
        v = invC * R
        y[i] = (1-γ) * mean(abs.(v))
    end
    return y
end

function compute_aggbound(perfs, moves, strats, idx2strat, strat2idx, alts, δ, bound=:DKW)
    A, Alow, Ahigh = make_perfs_bound(perfs, δ, bound)
    numA = length(perfs.algs)
    numE = length(perfs.envs)
    aidx(i) = Int((i-1) ÷ numE)+1
    eidx(i) = Int((i-1) % numE)+1

    function payoffs(n,s::Int,tail=:none)
        if tail==:none
            if n == 1
                return A[s]
            else
                return -A[s]
            end
        elseif tail==:left
            if n == 1
                return Alow[s]
            else
                return -Ahigh[s]
            end
        elseif tail==:right
            if n==1
                return Ahigh[s]
            else
                return -Alow[s]
            end
        elseif tail==:both
            if n ==1
                return Alow[s], Ahigh[s]
            else
                return -Ahigh[s], -Alow[s]
            end
        else
            throw(ArgumentError("tail=$(tail) is invalid"))
        end
    end

    players = length(moves)

    η = compute_eta(moves);
    min_p = 0.0
    m = 50
    numS = length(strats)
    C = zeros((numS, numS))
    fill!(C, 0.0)
    makeCinfinite_dense!(C, strats, alts, players, payoffs, η, m);

    R = zeros(numS)
    γ = (numS - 1.0) / Float64(numS)

    y = zeros(Float64, numA)
    ylow = zeros(Float64, numA)
    yhigh = zeros(Float64, numA)
    invC = inv(I - γ .* C)
    for i in 1:numA
        for s in 1:numS
            R[s] = A[i,eidx(idx2strat[s][2]),aidx(idx2strat[s][2])]
        end
        v = invC * R
        y[i] = (1-γ) * mean(abs.(v))
    end

    Clow, Chigh, Eu, num_uncertain = makeCinfinite_bounds(strats, alts, players, payoffs, η, m);

    for i in 1:numA
        for s in 1:numS
            R[s] = Alow[i,eidx(idx2strat[s][2]),aidx(idx2strat[s][2])]
        end
        ylow[i] = valueiteration_contagg(Clow, Chigh, Eu, R, numS, :min, γ=γ)
        for s in 1:numS
            R[s] = Ahigh[i,eidx(idx2strat[s][2]),aidx(idx2strat[s][2])]
        end
        yhigh[i] = valueiteration_contagg(Clow, Chigh, Eu, R, numS, :max, γ=γ)
    end
    return hcat(y, ylow, yhigh)
end

function bootstrap_aggbound(perfs, moves, strats, idx2strat, strat2idx, alts, δ, num_boot)
    bs = bootstrap((x)->(compute_aggregate(x, moves, strats, idx2strat, strat2idx, alts)),
                    perfs, BasicSampling(num_boot)
    )
    numA = length(perfs.algs)
    numE = length(perfs.envs)
    cfbs = confint(bs, PercentileConfInt(1 - (δ/(numA*numE))))
    Y = hcat([cfbs[i][1] for i in 1:numA], [cfbs[i][2] for i in 1:numA], [cfbs[i][3] for i in 1:numA])
    return Y
end

function load_env_result(base_dir, env, alg, colnames, perf_col)
    path = join([base_dir, env, string("allres_", alg, ".csv")], "/")
    if occursin("gw", env) || occursin("chain", env)
        header = colnames[4:end]
    else
        header = colnames
    end
    f = CSV.File(path, header=header)
    N = length(f)
    data = zeros(N)
    for (i, row) in enumerate(f)
        data[i] = row[perf_col]
    end
    return data
end

function combine_algs(base_dir, env, algs, acols, perf_col)
    allalg = [load_env_result(base_dir, env, alg, cols, perf_col) for (alg, cols) in zip(algs, acols)]
    return allalg
end

function split_data(data, nsplits, sizes, save_dir, env, algs)
    min_samples = sum(nsplits .* sizes)
    println(min_samples, " ", length.(data))
    @assert all(length.(data) .>= min_samples)
    foreach(shuffle!, data)
    iters = ones(Int, length(algs))

    for N in sizes
        for i in 1:nsplits
            df = DataFrame()
            for (j, alg) in enumerate(algs)
                idx = iters[j]
                df[!, Symbol(alg)] = data[j][idx:idx+(N-1)]
                iters[j] = idx + N
            end
            path = join([save_dir,string("allalgs_", env, "_$(lpad(N, 5, '0'))_$(lpad(i, 5, '0')).csv")], "/")
            CSV.write(path, df)
        end
    end
end

function load_and_split(base_dir, save_dir, env, algs, acols, perf_col, nsplits, sizes)
    data = combine_algs(base_dir, env, algs, acols, perf_col)
    split_data(data, nsplits, sizes, save_dir, env, algs)
end

function make_paths(base_dir, envs, N, split_idx)
    paths = []
    for env in envs
        path = join([base_dir, string("allalgs_", env, "_$(lpad(N, 5, '0'))_$(lpad(split_idx, 5, '0')).csv")], "/")
        push!(paths, path)
    end
    return paths
end

function get_aggbounds_samples(base_dir, envs, algs, δ, bound, sample_sizes)
    numA = length(algs)
    numE = length(envs)
    moves = [numA, numA*numE]
    strats, idx2strat, strat2idx, alts = make_strats(moves);
    results = []
    num_splits = 1000
    for ss in sample_sizes
        sres = zeros((num_splits, numA, 3))
        for i in 1:num_splits #Threads.@threads for i = 1:num_splits
            paths = make_paths(base_dir, envs, ss, i)
            perfs = load_results(paths, envs, algs)
            sres[i, :, :] .= compute_aggbound(perfs, moves, strats, idx2strat, strat2idx, alts, δ, bound)
        end
        push!(results, sres)
    end
    return results
end

function bootstrap_aggbounds_samples(base_dir, envs, algs, δ, num_boot, sample_sizes)
    numA = length(algs)
    numE = length(envs)
    moves = [numA, numA*numE]
    strats, idx2strat, strat2idx, alts = make_strats(moves);
    results = []
    num_splits = 1000
    for ss in sample_sizes
        sres = zeros((num_splits, numA, 3))
        Threads.@threads for i = 1:num_splits
            paths = make_paths(base_dir, envs, ss, i)
            perfs = load_results(paths, envs, algs)
            sres[i, :, :] .= bootstrap_aggbound(perfs, moves, strats, idx2strat, strat2idx, alts, δ, num_boot)
        end
        push!(results, sres)
    end
    return results
end

function num_bound_violations(res)
    N = size(res)[1] # num trials
    violations = 0
    for i in 1:N
        perf = res[i, :, 1]
        for j in 1:N
            if i==j
                continue
            end
            below = perf .< res[j,:,2]
            above = perf .> res[j,:,3]
            if any(above) || any(below)
                violations += 1
            end
        end
    end
    return violations / (N * (N-1))
end

function num_significant(res)
    N, numA = size(res)[1:2] # num trials, num algorithms
    num_sig = 0
    for i in 1:N
        for j in 1:(numA-1)
            for k in j:numA
                jlow = res[i, j, 2]
                jhigh = res[i, j, 3]
                klow = res[i, k, 2]
                khigh = res[i, k, 3]
                num_sig += (jlow > khigh) | (jhigh < klow)
            end
        end
    end
    numComp = (numA * (numA-1)) / 2
    return num_sig / (N * numComp)
end

function bound_res(res)
    violations = zeros(size(res)[1])
    sigs = zeros(size(res)[1])
    for i in 1:length(violations)
        violations[i] = num_bound_violations(res[i])
        sigs[i] = num_significant(res[i])
    end
    return violations, sigs
end

function main()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--conf", "-c"
            help = "confidence level"
            arg_type = Float64
            required = true
        "--dir", "-d"
            help = "directory"
            arg_type = String
            required = true
        "--sdir", "-s"
            help = "folder directory to store results"
            arg_type = String
            required = true
        "--boot", "-b"
            help = "number of bootstrap samples to use"
            arg_type = Int
            required = false
            default = 1000
    end

    parsed_args = parse_args(ARGS, s)
    base_dir = parsed_args["dir"]
    save_dir = parsed_args["sdir"]
    δ = parsed_args["conf"]
    num_boot = parsed_args["boot"]
    # δ = 0.05
    # num_boot = 10000
    # base_dir = "/home/sjordan/Documents/icmleval/complete/bounds/splits";
    # save_dir = "/home/sjordan/Documents/icmleval/complete/bounds/intervals2";
    println("data dir: ", base_dir)
    println("save dir: ", save_dir)
    println("confidence level: ", δ)
    println("num_bootstraps: ", num_boot)

    algs_bt = ["ac-parl2", "sarsa-scaled", "sarsa-parl2", "q-parl2"]
    envs_bt = ["mntcar_d", "chain10d", "chain10s", "chain50d", "chain50s", "gw10d", "gw10s", "gw5d", "gw5s", ]

    acp_cols = ["dorder", "iorder", "full", "gamma", "lam", "palpha", "logp", "life", "end"]
    qlp_cols = ["dorder", "iorder", "full", "gamma", "lam", "eps", "logp", "life", "end"]
    sls_cols = ["dorder", "iorder", "full", "gamma", "lam", "eps", "qalpha", "logp", "life", "end"]
    slp_cols = ["dorder", "iorder", "full", "gamma", "lam", "eps", "logp", "life", "end"]
    alg_cols_bt = [acp_cols, sls_cols, slp_cols, qlp_cols];
    sample_sizes = [10, 30, 100, 1000, 10000]

    DKW_res = get_aggbounds_samples(base_dir, envs_bt, algs_bt, δ, :DKW, sample_sizes);
    path = join([save_dir, "dkw3.jld"], "/")
    save(path, "DKW_res", DKW_res)
    # @save path DKW_res
    #@load path DKW_res
    dkw_bv, dkw_sig = bound_res(DKW_res)
    println(dkw_bv)
    println(dkw_sig)
    println("got DKW bounds")
    flush(stdout)

    TT_res = get_aggbounds_samples(base_dir, envs_bt, algs_bt, δ, :TTest, sample_sizes);
    path = join([save_dir, "tt3-AM.jld"], "/")
    save(path, "TT_res", TT_res)
    # @save path TT_res
    # @load path TT_res
    tt_bv, tt_sig = bound_res(TT_res)
    println(tt_bv)
    println(tt_sig)
    println("got TTest bounds")
    flush(stdout)

    TT_res2 = get_aggbounds_samples(base_dir, envs_bt, algs_bt, δ, :TTest2, sample_sizes);
    path = join([save_dir, "tt3-AAM.jld"], "/")
    save(path, "TT_res2", TT_res2)
    # @save path TT_res2
    # @load path TT_res2
    tt2_bv, tt2_sig = bound_res(TT_res2)
    println(tt2_bv)
    println(tt2_sig)
    println("got TTest2 bounds")
    flush(stdout)

    #@time BS_res = bootstrap_aggbounds_samples(base_dir, envs_bt, algs_bt, δ, num_boot, sample_sizes);
    #path = join([save_dir, "bs.jld"], "/")
    #save(path, "BS_res", BS_res)
    #@save path BS_res
    # @load path BS_res

    println("got Bootstrap bounds")
    flush(stdout)

    # dkw_bv, dkw_sig = bound_res(DKW_res)
    # tt_bv, tt_sig = bound_res(TT_res)
    # tt2_bv, tt2_sig = bound_res(TT_res2)
    #bs_bv, bs_sig = bound_res(BS_res);
    println("got all result measures")
    flush(stdout)

    bvdf = DataFrame(SampleSize = sample_sizes,
        vDKW=dkw_bv, vTTest=tt_bv, vTTest2=tt2_bv,
        #vBootstrap=bs_bv,
        sDKW=dkw_sig, sTTest=tt_sig, sTTest2=tt2_sig,
        # sBootstrap=bs_sig
                )
    println("Bound results")
    println(bvdf)

end

main()

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# Time series having length == number_of_hours * number_of_scenarios
function make_time_series(source::AbstractDict, lookup::AbstractDict, y::Int, sn_data::AbstractDict; number_of_hours::Int, number_of_scenarios::Int)
    target = Dict{String,Any}(
        "gen"        => Dict{String,Any}(),
        "load"       => Dict{String,Any}(),
        "storage"    => Dict{String,Any}(),
        "ne_storage" => Dict{String,Any}(),
    )

    for comp in source["scenarioDataInputFile"]["generators"]
        index = lookup["generators"][comp["id"]]
        if haskey(comp, "capacityFactor")
            p = ts_vector(comp, "capacityFactor", y; number_of_hours, number_of_scenarios) .* sn_data["gen"]["$index"]["pmax"]
            target["gen"]["$index"] = Dict{String,Any}("pmax" => p)
        end
    end

    for comp in source["scenarioDataInputFile"]["loads"]
        index = lookup["loads"][comp["id"]]
        pd    = ts_vector(comp, "demandReference", y; number_of_hours, number_of_scenarios)
        target["load"]["$index"] = Dict{String,Any}("pd" => pd)
    end

    for comp in source["scenarioDataInputFile"]["storage"]
        index = lookup["storage"][comp["id"]]
        target["storage"]["$index"] = Dict{String,Any}()
        if haskey(comp, "powerExternalProcess")
            p = ts_vector(comp, "powerExternalProcess", y; number_of_hours, number_of_scenarios)
            target["storage"]["$index"]["stationary_energy_inflow"] = max.(p,0.0)
            target["storage"]["$index"]["stationary_energy_outflow"] = -min.(p,0.0)
        else
            target["storage"]["$index"]["stationary_energy_inflow"] = zeros(number_of_hours*number_of_scenarios)
            target["storage"]["$index"]["stationary_energy_outflow"] = zeros(number_of_hours*number_of_scenarios)
        end
        if haskey(comp, "maxAbsActivePower")
            target["storage"]["$index"]["charge_rating"] = ts_vector(comp, "maxAbsActivePower", y; number_of_hours, number_of_scenarios) * sn_data["storage"]["$index"]["charge_rating"]
        end
        if haskey(comp, "maxInjActivePower")
            target["storage"]["$index"]["discharge_rating"] = ts_vector(comp, "maxInjActivePower", y; number_of_hours, number_of_scenarios) * sn_data["storage"]["$index"]["discharge_rating"]
        end
    end

    if haskey(source, "candidatesInputFile")
        for cand in source["candidatesInputFile"]["storage"]
            comp = cand["storageData"]
            index = lookup["cand_storage"][comp["id"]]
            target["ne_storage"]["$index"] = Dict{String,Any}()
            if haskey(comp, "powerExternalProcess")
                p = ts_vector(comp, "powerExternalProcess", y; number_of_hours, number_of_scenarios)
                target["ne_storage"]["$index"]["stationary_energy_inflow"] = max.(p,0.0)
                target["ne_storage"]["$index"]["stationary_energy_outflow"] = -min.(p,0.0)
            else
                target["ne_storage"]["$index"]["stationary_energy_inflow"] = zeros(number_of_hours*number_of_scenarios)
                target["ne_storage"]["$index"]["stationary_energy_outflow"] = zeros(number_of_hours*number_of_scenarios)
            end
            if haskey(comp, "maxAbsActivePower")
                target["ne_storage"]["$index"]["charge_rating"] = ts_vector(comp, "maxAbsActivePower", y; number_of_hours, number_of_scenarios) * sn_data["ne_storage"]["$index"]["charge_rating"]
            end
            if haskey(comp, "maxInjActivePower")
                target["ne_storage"]["$index"]["discharge_rating"] = ts_vector(comp, "maxInjActivePower", y; number_of_hours, number_of_scenarios) * sn_data["ne_storage"]["$index"]["discharge_rating"]
            end
        end
    end

    return target
end

# Time series data vector from JSON component dict
function ts_vector(comp::AbstractDict, key::String, y::Int; number_of_hours::Int, number_of_scenarios::Int)
    # `comp[key]` is a Vector{Vector{Vector{Any}}} containing data of each scenario, year and hour
    # Returned value is a Vector{Float64} containing data of year y and each scenario and hour
    return mapreduce(Vector{Float64}, vcat, comp[key][s][y][1:number_of_hours] for s in 1:number_of_scenarios)
end


text:
# univariate
function enclose_BranchandBound(f::Function, dom::Interval; order=10, tol=0.6)
    x0 = Interval(mid(dom))
    x = TaylorModel1(order, x0, dom)
    return branchandbound(f(x-x0).pol, dom, tol)
end

# multivariate
function enclose_BranchandBound(f::Function, dom::IntervalBox{D};
                                order=10, tol=0.6) where {D}
    x0 = mid(dom)
    set_variables(Float64, "x", order=2order, numvars=D)
    x = [TaylorModelN(i, order, IntervalBox(x0), dom) for i=1:D]
    return branchandbound(f(x...).pol, dom - x0, tol)
end

@inline _Rnext(R::Vector{<:Interval}) = Interval(minimum(inf.(R)), maximum(sup.(R)))

const Kmax = 1000

function branchandbound(p::Union{Taylor1, TaylorN},
                        dom::Interval_or_IntervalBox,
                        tol::Number)
    K = 1
    Rprev = evaluate(p, dom)
    D = bisect(dom)
    R = [evaluate(p, D[1]), evaluate(p, D[2])]
    Rnext = R[1] ∩ R[2]

    while  (Rprev.hi - Rnext.hi) >= tol * diam(Rnext) &&
           (Rprev.lo - Rnext.lo) >= tol * diam(Rnext)

        Rprev = _Rnext(R)
        R_x = sup.(R)
        R_n = inf.(R)
        max_range = maximum(R_x)
        max_index = findall(x->x == max_range, R_x)[1]
        min_range = minimum(R_n)
        min_index = findall(x->x == min_range, R_n)[1]
        l_D = length(D[1])
        K += 1
        if K == Kmax
            error("convergence not achieved in branch and bound after $Kmax steps")
        end

        D, R = divide_dom!(p, D, R, max_index)
        if max_index < min_index
            min_index += 1
            D, R = divide_dom!(p, D, R, min_index)
        end
        Rnext = _Rnext(R)
    end
    return _Rnext(R)
end

function divide_dom!(p::Union{TaylorN{T}, Taylor1{T}},
                     D::Union{Array{IntervalBox{N, W}, M}, Array{Interval{W}, M}},
                     R::Array{Interval{W}, M}, index::Number) where {N, T, M, W}
    BA = [ ((D[index][i]).hi - (D[index][i]).lo) for i = 1:length(D[1])]
    Beta1 = maximum(BA)
    β = findall(x->x == Beta1, BA)[1]
    D1, D2 = bisect(D[index], β)
    D[index] = D1
    DD = push!(D[1:index], D2)
    D = append!(DD, D[(index + 1):length(D)])
    R[index] = evaluate(p, D[index])
    RR = append!(R[1:index], evaluate(p, D[index + 1]))
    R = append!(RR, R[(index + 1):length(R)])
    return D, R
end

function enclose_binary(f, dom::Interval; kmax=3, tol=1e-3, algorithm=:IntervalArithmetic)
    y = enclose(f, dom, algorithm)
    yinf, ysup = inf(y), sup(y)
    kmax == 0 && return Interval(yinf, ysup)
    x = bisect(dom)
    fx1 = enclose(f, x[1], algorithm)
    fx2 = enclose(f, x[2], algorithm)
    ynew = hull(fx1, fx2)
    ynew_inf, ynew_sup = inf(ynew), sup(ynew)
    inf_close = abs(yinf - ynew_inf) <= tol
    sup_close = abs(ysup - ynew_sup) <= tol
    both_close = inf_close && sup_close
    inf_improves = ynew_inf > yinf
    sup_improves = ynew_sup < ysup
    both_improve = inf_improves && sup_improves
    if both_close || !both_improve
        return Interval(yinf, ysup)
    end
    yinf = max(yinf, ynew_inf)
    ysup = min(ysup, ynew_sup)
    if inf_improves
        yinf = ynew_inf
    end
    if sup_improves
        ysup = ynew_sup
    end
    e1 = enclose_binary(f, x[1], kmax=kmax-1, algorithm=algorithm)
    e2 = enclose_binary(f, x[2], kmax=kmax-1, algorithm=algorithm)
    return Interval(hull(e1, e2))
end


text:
function fractran(n::Integer, ratios::Vector{<:Rational}, steplim::Integer)
    rst = zeros(BigInt, steplim)
    for i in 1:steplim
        rst[i] = n
        if (pos = findfirst(x -> isinteger(n * x), ratios)) > 0
            n *= ratios[pos]
        else
            break
        end
    end
    return rst
end

using IterTools
macro ratio_str(s)
    a = split(s, r"[\s,/]+")
    return collect(parse(BigInt, n) // parse(BigInt, d) for (n, d) in partition(a, 2))
end

fracs = ratio"""17 / 91, 78 / 85, 19 / 51, 23 / 38, 29 / 33, 77 / 29, 95 / 23,
                77 / 19, 1 / 17, 11 / 13, 13 / 11, 15 / 14, 15 / 2, 55 / 1"""
println("The first 20 in the series are ", fractran(2, fracs, 20))

prmfound = 0
n = big(2)
while prmfound < 20
    if isinteger(log2(n))
        prmfound += 1
        println("Prime $prmfound found: $n is 2 ^ $(Int(log2(n)))")
    end
    n = fractran(n, fracs, 2)[2]
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode