text:
@with_kw struct SetupOptions
    allow_caching::Bool = true
    return_quotes::Bool = false
    background_only::Bool = false
    Ginv::Function = (G,f,params) -> inv(G(f,params))
    christoffel::Function = christoffel
    simplification::Int = 1 # level of simplification, bigger is more aggressive
    # the simplification function, try simplify(x) for small Nf
    # factor(expand(x)) is more efficient at higher Nf
    # see https://docs.sympy.org/latest/tutorial/simplification.html for more
    simp::Function = x-> cancel(x)
    attempt_HCA::Bool = true # attempt to apply the horizon-crossing approximation
end
@with_kw struct BackgroundOptions
    background_atol::Float64 = 1e-14
    background_rtol::Float64 = 1e-14
    max_efolds::Float64 = 1e4
    min_efolds::Float64 = 0.0 # force eom evolution even if ϵ > 1 for this long
    max_steps::Float64 = Int(1e8)
    background_solver = Tsit5()
    end_epsilon::Float64 = 1.0 # the value of ϵ that ends inflation
    save_everystep::Bool = false
    saveat::Float64 = 1e-3
    Nexit_to_end::Float64 = 55
end
@with_kw struct PerturbationOptions
    # perturbations
    force_transport::Bool = true # compute perturbations even if inflation doesn't end
    compute_isocurvature::Bool = true # when d > 1, compute scalar powerspectra from all non-adiabatic 2-point functions
    compute_tensors::Bool = true # compute tensor perturbations
    Nsubhorizon::Float64 = 8. # number of efolds before horizon exit to impose Bunch-Davies
    Nkrange::Float64 = 5.0 # number of e-folds spanned by k-values
    Nexit_to_end::Float64 = 55
    kpivot::Float64 = 0.002 # Mpc^(-1)
    desired_ksteps::Int = 6 # odd number recommended, to measure ns right at pivot scale
    # scalar eom
    scalars_max_steps::Int = 4e8
    scalars_atol::Float64 = 1e-13
    scalars_rtol::Float64 = 1e-13
    # picking the right solver can be very problem dependent. https://docs.juliadiffeq.org/stable/tutorials/advanced_ode_example/#Choosing-a-Good-Solver-1
    # if Tsit5 fails, try lowering the atol and rtol and switching to Rodas5
    # All worth a try: Rosenbrock23, Rodas5, lsoda, TRBDF2, KenCarp4
    scalars_solver = Tsit5()
    save_everystep::Bool = false
    saveat::Float64 = 1e-1
    # tensor eom
    tensors_max_steps::Int = scalars_max_steps
    tensors_atol::Float64 = scalars_atol
    tensors_rtol::Float64 = scalars_rtol
    tensors_solver = scalars_solver
end


text:
#! /usr/bin/julia

# Rosetta Code, Dutch national flag problem

const COLORS = ["red", "white", "blue"]

function dutchsort!(a::Array{ASCIIString,1}, lo=COLORS[1], hi=COLORS[end])
    i = 1
    j = 1
    n = length(a)
    while j <= n
        if a[j] == lo
            a[i], a[j] = a[j], a[i]
            i += 1
            j += 1
        elseif a[j] == hi
            a[j], a[n] = a[n], a[j]
            n -= 1
        else
            j += 1
        end
    end
    return a
end

function dutchsort(a::Array{ASCIIString,1}, lo=COLORS[1], hi=COLORS[end])
    dutchsort!(copy(a), lo, hi)
end

function formatdf(a::Array{ASCIIString,1})
    i = 0
    s = "    "
    for c in a
        s *= @sprintf "%6s" c
        i += 1
        i %= 8
        if i == 0
            s *= "\n    "
        end
    end
    return s
end

cnum = 20
d = [COLORS[rand(1:3)] for i in 1:cnum]
while d == dutchsort(d)
    d = [COLORS[rand(1:3)] for i in 1:cnum]
end

println("The original list is:")
println(formatdf(d))

print("Sorting with dutchsort, ")
@time e = dutchsort(d)
println(formatdf(e))

print("Sorting conventionally, ")
@time e = sort(d, by=x->findfirst(COLORS, x))
println(formatdf(e))

ext: jl

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
mutate!(ind::indiv, fmut!) = fmut!(ind.x)

function rand_flip!(bits)
    nb = length(bits)
    for i = 1:nb
        if rand() < 1/nb
            bits[i] = 1 - bits[i]
        end
    end
end
default_mutation!(b::BitArray) = rand_flip!(b)
default_mutation!(b::Vector{Bool}) = rand_flip!(b)

function rand_swap!(perm::Vector{Int})
    i = j = rand(1:length(perm))
    while j == i
        j = rand(1:length(perm))
    end
    perm[i], perm[j] = perm[j], perm[i]
end
default_mutation!(p::Vector{Int}) = rand_swap!(p)


text:
@testset "EGPnormal" begin
    @test_throws ArgumentError EGPnormal(-1, .2, 2)
    @test_throws ArgumentError EGPnormal(-1, .2, -2)

    σ = 1
    ξ = -0.1
    κ = 2
    model = EGPnormal(σ, ξ, κ)
    @test minimum(model) == 0
    @test maximum(model) == -σ/ξ

    σ = 1
    ξ = .1
    κ = 2
    model = EGPnormal(σ, ξ, κ)

    @test maximum(model) == Inf

    @test insupport(model, -1) == false
    @test insupport(model, 0) == true
    @test insupport(model, 1000) == true

    @test params(model) == (σ, ξ, κ)

    @test partype(model) == Float64

    @test_throws AssertionError quantile(model, 0)
    @test_throws AssertionError quantile(model, 1)
end


text:
function default_args(str)
    buf = codeunits(str)
    len = length(buf)
    return buf, 1, 1, 1, len
end

@testset "name" begin
    str = "name"
    @test GraphQLParser.parse_name(default_args(str)...)[1] == "name"

    str = "n"
    @test GraphQLParser.parse_name(default_args(str)...)[1] == "n"

    str = "1name"
    @test_throws ArgumentError GraphQLParser.parse_name(default_args(str)...)
end

@testset "Selection" begin
    str = "name"
    @test GraphQLParser.parse_selection(default_args(str)...)[1] == Field(nothing, "name", nothing, nothing, nothing)

    str = "alias: name"
    @test GraphQLParser.parse_selection(default_args(str)...)[1] == Field("alias", "name", nothing, nothing, nothing)
end

@testset "SelectionSet" begin
    str = "{name}"
    @test GraphQLParser.parse_selection_set(default_args(str)...)[1] == SelectionSet([Field(nothing, "name", nothing, nothing, nothing)])

    str = "{name1, name2}"
    @test GraphQLParser.parse_selection_set(default_args(str)...)[1] == SelectionSet([
        Field(nothing, "name1", nothing, nothing, nothing),
        Field(nothing, "name2", nothing, nothing, nothing),
    ])

    str = "{name1, alias2: name2{subfield\n}}"
    @test GraphQLParser.parse_selection_set(default_args(str)...)[1] == SelectionSet([
        Field(nothing, "name1", nothing, nothing, nothing),
        Field("alias2", "name2", nothing, nothing, SelectionSet([Field(nothing, "subfield", nothing, nothing, nothing)])),
    ])
end

@testset "number" begin
    str = "0"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 0

    str = "-0"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 0

    str = "11111"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 11111

    str = "0.0"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 0.0

    str = "0.01"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 0.01

    str = "-0.01"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === -0.01

    str = "-0.01e10"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === -0.01e10

    str = "1E+10"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 1e10

    str = "1E-10  asd"
    @test GraphQLParser.parse_number(default_args(str)...)[1] === 1e-10
end

@testset "Argument" begin
    str = "number: 1E-10"
    @test GraphQLParser.parse_argument(default_args(str)...)[1] == Argument("number", 1e-10)

    str = "bool: true"
    @test GraphQLParser.parse_argument(default_args(str)...)[1] == Argument("bool", true)

    str = "bool: false"
    @test GraphQLParser.parse_argument(default_args(str)...)[1] == Argument("bool", false)

    str = "null: null"
    @test GraphQLParser.parse_argument(default_args(str)...)[1] == Argument("null", nothing)

    str = "var: \$variable_name"
    @test GraphQLParser.parse_argument(default_args(str)...)[1] == Argument("var", Variable("variable_name"))
end

@testset "Arguments" begin
    str = "(number: 1E-10)"
    @test GraphQLParser.parse_arguments(default_args(str)...)[1] == [Argument("number", 1e-10)]

    str = "(number: 1E-10 bool: true)"
    @test GraphQLParser.parse_arguments(default_args(str)...)[1] == [Argument("number", 1e-10), Argument("bool", true)]

    @test_throws ArgumentError GraphQLParser.parse_arguments(default_args("()")...)
end

@testset "string" begin
    @testset "one line" begin
        str = raw"\"str\""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "str"

        str = raw"\"str\" \"str\""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "str"

        str = raw"\"str\\ \""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "str\\ "

        str = raw"\"\t\""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "\t"

        str = "\"\t\""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "\t"

        str = "\"\t\\\"\""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "\t\""

        str = raw"\"\u3949\""
        @test GraphQLParser.parse_string(default_args(str)...)[1] == "\u3949"
    end

    @testset "block" begin
        buf = [0x22, 0x22, 0x22, 0x61, 0x22, 0x22, 0x22] # """a"""
        len = length(buf)
        @test GraphQLParser.parse_string(buf, 1, 1, 1, len)[1] == "a"

        buf = [0x22, 0x22, 0x22, 0x61, 0x09, 0x22, 0x22, 0x22] # """a\t"""
        len = length(buf)
        @test GraphQLParser.parse_string(buf, 1, 1, 1, len)[1] == "a\t"

        buf = [0x22, 0x22, 0x22, 0x5c, 0x22, 0x22, 0x22, 0x22, 0x22, 0x22] # """\""""""
        len = length(buf)
        @test GraphQLParser.parse_string(buf, 1, 1, 1, len)[1] == "\"\"\""

        buf = [0x22, 0x22, 0x22, 0x5c, 0x22, 0x22, 0x22, 0x22, 0x22, 0x22] # """\""""""
        len = length(buf)
        @test GraphQLParser.parse_string(buf, 1, 1, 1, len)[1] == "\"\"\""
    end
end

@testset "Lists" begin
    str = "[1,2,3]"
    @test GraphQLParser.parse_list(default_args(str)...)[1] == Any[1,2,3]

    str = """["a","b","c"]"""
    @test GraphQLParser.parse_list(default_args(str)...)[1] == Any["a","b","c"]

    str = "[[1,2],[2,3]]"
    @test GraphQLParser.parse_list(default_args(str)...)[1] == Any[Any[1,2], Any[2,3]]

    str = "[]"
    @test GraphQLParser.parse_list(default_args(str)...)[1] == Any[]
end

@testset "Input Objects" begin
    str = "{}"
    @test GraphQLParser.parse_input_object(default_args(str)...)[1] isa GraphQLParser.InputObject

    str = "{field_name: 1}"
    @test GraphQLParser.parse_input_object(default_args(str)...)[1] isa GraphQLParser.InputObject

    str = "{field_name: {subfield: \"str\"}}"
    GraphQLParser.parse_input_object(default_args(str)...)[1]
end

@testset "type" begin
    str = "Int"
    @test GraphQLParser.parse_type(default_args(str)...)[1] == "Int"

    str = "Int!"
    @test GraphQLParser.parse_type(default_args(str)...)[1] == "Int!"

    str = "[[[Int!]]]!"
    @test GraphQLParser.parse_type(default_args(str)...)[1] == "[[[Int!]]]!"
end

@testset "VariableDefinition" begin
    str = "(\$name: Int)"
    @test GraphQLParser.parse_variable_definitions(default_args(str)...)[1] == [VariableDefinition("name", "Int", nothing, nothing)]

    str = "(\$name: Int=1)"
    @test GraphQLParser.parse_variable_definitions(default_args(str)...)[1] == [VariableDefinition("name", "Int", 1, nothing)]

    str = "(\$name: Int=1, \$name2: [MyType]!)"
    @test GraphQLParser.parse_variable_definitions(default_args(str)...)[1] == [VariableDefinition("name", "Int", 1, nothing), VariableDefinition("name2", "[MyType]!", nothing, nothing)]
end

@testset "directives" begin
    str = "@skip"
    @test GraphQLParser.parse_directive(default_args(str)...)[1] == Directive("skip", nothing)

    str = "@skip(condition: \$val)"
    @test GraphQLParser.parse_directive(default_args(str)...)[1] == Directive("skip", [Argument("condition", Variable("val"))])
end

@testset "fragments" begin
    @testset "definition" begin
        # Example 19 from 2021 spec
        str = """
        fragment friendFields on User {
            id
            name
            profilePic(size: 50)
        }"""
        intended_result = FragmentDefinition(
            "friendFields",
            "User",
            nothing,
            SelectionSet(
                [
                    Field(nothing, "id", nothing, nothing, nothing),
                    Field(nothing, "name", nothing, nothing, nothing),
                    Field(nothing, "profilePic", [Argument("size", 50)], nothing, nothing),
                ]
            )
        )
        @test GraphQLParser.parse_fragment_definition(default_args(str)...)[1] == intended_result

        str = """
        fragment friendFields on User @skip(cond: \$cond){
            id
            name
            profilePic(size: 50)
        }"""
        intended_result = FragmentDefinition(
            "friendFields",
            "User",
            [
                Directive("skip", [Argument("cond", Variable("cond"))])
            ],
            SelectionSet(
                [
                    Field(nothing, "id", nothing, nothing, nothing),
                    Field(nothing, "name", nothing, nothing, nothing),
                    Field(nothing, "profilePic", [Argument("size", 50)], nothing, nothing),
                ]
            )
        )
        @test GraphQLParser.parse_fragment_definition(default_args(str)...)[1] == intended_result
    end

    @testset "spread" begin
        # From example 19 from 2021 spec
        str = """...friendFields("""
        @test GraphQLParser.parse_fragment(default_args(str)...)[1] == FragmentSpread("friendFields", nothing)
    end

    @testset "inline" begin
        # From Example 23, 2021 spec
        str = """... on User {
            friends {
                count
            }
        }
        """
        intended_result = InlineFragment(
            "User",
            nothing,
            SelectionSet([
                Field(nothing, "friends", nothing, nothing, SelectionSet([Field(nothing, "count", nothing, nothing, nothing)]))
            ])
        )
        @test GraphQLParser.parse_fragment(default_args(str)...)[1] == intended_result

        str = """... {
            friends {
                count
            }
        }
        """
        intended_result = InlineFragment(
            nothing,
            nothing,
            SelectionSet([
                Field(nothing, "friends", nothing, nothing, SelectionSet([Field(nothing, "count", nothing, nothing, nothing)]))
            ])
        )
        @test GraphQLParser.parse_fragment(default_args(str)...)[1] == intended_result

        # From Example 24, 2021 spec
        str = """... @include(if: \$expandedInfo) {
            firstName
            lastName
            birthday
            }
        """
        intended_result = InlineFragment(
            nothing,
            [Directive("include", [Argument("if", Variable("expandedInfo"))])],
            SelectionSet([
                Field(nothing, "firstName", nothing, nothing, nothing),
                Field(nothing, "lastName", nothing, nothing, nothing),
                Field(nothing, "birthday", nothing, nothing, nothing),
            ])
        )
        @test GraphQLParser.parse_fragment(default_args(str)...)[1] == intended_result
    end
end

@testset "comments" begin
    str1 = """
    {
        myQuery{
            field # this is a useful field
            # not_this_field
        }
    }
    # finishing comments
    """

    str2 = """
    {
        myQuery{
            field
        }
    }
    """

    @test GraphQLParser.parse(str1) == GraphQLParser.parse(str2)
end

@testset "Full Examples" begin
    # Example No 5
    str = """
    mutation {
        likeStory(storyID: 12345) {
            story {
                likeCount
            }
        }
    }
    """
    @test GraphQLParser.parse(str) == Document([
        Operation(
            "mutation",
            nothing,
            nothing,
            nothing,
            SelectionSet([Field(
                nothing,
                "likeStory",
                [Argument("storyID", 12345)],
                nothing,
                SelectionSet([Field(
                    nothing,
                    "story",
                    nothing,
                    nothing,
                    SelectionSet([Field(
                        nothing,
                        "likeCount",
                        nothing,
                        nothing,
                        nothing
                    )])
                )])
            )])
        )
    ])

    # Example No 6
    str = """{
        field
    }"""
    @test GraphQLParser.parse(str) == Document([
        Operation(
            "query",
            nothing,
            nothing,
            nothing,
            SelectionSet([Field(nothing, "field", nothing, nothing, nothing)])
        )
    ])

    # Example 19
    str = """query withFragments {
        user(id: 4) {
            friends(first: 10) {
                ...friendFields
            }
            mutualFriends(first: 10) {
                ...friendFields
            }
        }
    }

    fragment friendFields on User {
        id
        name
        profilePic(size: 50)
    }
    """
    @test GraphQLParser.parse(str) == Document([
        Operation(
            "query",
            "withFragments",
            nothing,
            nothing,
            SelectionSet([
                Field(
                    nothing,
                    "user",
                    [Argument("id", 4)],
                    nothing,
                    SelectionSet([
                        Field(
                            nothing,
                            "friends",
                            [Argument("first", 10)],
                            nothing,
                            SelectionSet([FragmentSpread("friendFields", nothing)])
                        ),
                        Field(
                            nothing,
                            "mutualFriends",
                            [Argument("first", 10)],
                            nothing,
                            SelectionSet([FragmentSpread("friendFields", nothing)])
                        )
                    ])
                )
            ])
        ),
        FragmentDefinition(
            "friendFields",
            "User",
            nothing,
            SelectionSet([
                Field(nothing, "id", nothing, nothing, nothing),
                Field(nothing, "name", nothing, nothing, nothing),
                Field(nothing, "profilePic", [Argument("size", 50)], nothing, nothing),
            ])
        )
    ])

    str = """
    query withNestedFragments {
        user(id: 4) {
            friends(first: 10) {
                ...friendFields
            }
            mutualFriends(first: 10) {
                ...friendFields
            }
        }
    }

    fragment friendFields on User {
        id
        name
        ...standardProfilePic
    }

    fragment standardProfilePic on User {
        profilePic(size: 50)
    }
    """
    @test GraphQLParser.parse(str) == Document([
        Operation(
            "query",
            "withNestedFragments",
            nothing,
            nothing,
            SelectionSet([
                Field(
                    nothing,
                    "user",
                    [Argument("id", 4)],
                    nothing,
                    SelectionSet([
                        Field(
                            nothing,
                            "friends",
                            [Argument("first", 10)],
                            nothing,
                            SelectionSet([FragmentSpread("friendFields", nothing)])
                        ),
                        Field(
                            nothing,
                            "mutualFriends",
                            [Argument("first", 10)],
                            nothing,
                            SelectionSet([FragmentSpread("friendFields", nothing)])
                        )
                    ])
                )
            ])
        ),
        FragmentDefinition(
            "friendFields",
            "User",
            nothing,
            SelectionSet([
                Field(nothing, "id", nothing, nothing, nothing),
                Field(nothing, "name", nothing, nothing, nothing),
                FragmentSpread("standardProfilePic", nothing),
            ])
        ),
        FragmentDefinition(
            "standardProfilePic",
            "User",
            nothing,
            SelectionSet([
                Field(nothing, "profilePic", [Argument("size", 50)], nothing, nothing),
            ])
        )
    ])

    # Examples 25 and 26
    str_25 = """
    mutation {
        sendEmail(message: \"\"\"
            Hello,
                World!

            Yours,
                GraphQL.
        \"\"\")
    }"""
    @test GraphQLParser.parse(str_25) == Document([
        Operation(
            "mutation",
            nothing,
            nothing,
            nothing,
            SelectionSet([Field(
                    nothing,
                    "sendEmail",
                    [Argument("message", "Hello,\n    World!\n\nYours,\n    GraphQL.")],
                    nothing,
                    nothing
            )])
        )
    ])

    str_26 = """mutation {
        sendEmail(message: "Hello,\n    World!\n\nYours,\n    GraphQL.")
    }
    """
    @test GraphQLParser.parse(str_25) == GraphQLParser.parse(str_26)

    # Check for anonymous query second
    str = """
    fragment Fragment on User {name}

    {
        myQuery
    }
    """
    @test GraphQLParser.parse(str) == Document([
        FragmentDefinition(
            "Fragment",
            "User",
            nothing,
            SelectionSet([Field(nothing, "name", nothing, nothing, nothing)]),
        ),
        Operation(
            "query",
            nothing,
            nothing,
            nothing,
            SelectionSet([Field(nothing, "myQuery", nothing, nothing, nothing)]),
        )
    ])
end

@testset "Line and column" begin
    str = "name"
    name, pos, line, column = GraphQLParser.parse_name(default_args(str)...)
    @test pos == 5
    @test line == 1
    @test column == 5
    str = "{\nname\n}"
    name, pos, line, column = GraphQLParser.parse_selection_set(default_args(str)...)
    @test pos == 9
    @test line == 3
    @test column == 2
    str = "{a(b:\"\"\"\nn\n\nh\n\"\"\")}"
    name, pos, line, column = GraphQLParser.parse_selection_set(default_args(str)...)
    @test pos == 20
    @test line == 5
    @test column == 6

    # Example No 5
    str = """
    mutation {
        likeStory(storyID: 12345) {
            story {
                likeCount
            }
        }
    }
    """
    doc = GraphQLParser.parse(str)
    @test doc.definitions[1].loc == Loc(1,1)
    @test doc.definitions[1].selection_set.loc == Loc(1,10)
    @test doc.definitions[1].selection_set.selections[1].loc == Loc(2,5)
    @test doc.definitions[1].selection_set.selections[1].arguments[1].loc == Loc(2,15)
    @test doc.definitions[1].selection_set.selections[1].selection_set.loc == Loc(2,31)
    @test doc.definitions[1].selection_set.selections[1].selection_set.selections[1].loc == Loc(3,9)
    @test doc.definitions[1].selection_set.selections[1].selection_set.selections[1].selection_set.loc == Loc(3,15)
    @test doc.definitions[1].selection_set.selections[1].selection_set.selections[1].selection_set.selections[1].loc == Loc(4,13)
    str = """
    query myquery{
        alias: name @skip(if: \$var){
            ...FragmentDef
            ... on User{field}
        }
    }
    """
end


text:
function BinomialDistribution(arg0::jint, arg1::jdouble)
    return BinomialDistribution((jint, jdouble), arg0, arg1)
end

function cumulative_probability(obj::BinomialDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_number_of_trials(obj::BinomialDistribution)
    return jcall(obj, "getNumberOfTrials", jint, ())
end

function get_numerical_mean(obj::BinomialDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::BinomialDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_probability_of_success(obj::BinomialDistribution)
    return jcall(obj, "getProbabilityOfSuccess", jdouble, ())
end

function get_support_lower_bound(obj::BinomialDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::BinomialDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::BinomialDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_probability(obj::BinomialDistribution, arg0::jint)
    return jcall(obj, "logProbability", jdouble, (jint,), arg0)
end

function probability(obj::BinomialDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end


text:
module RipsererPrep

using DataFrames
using CSV

# Write your package code here.

export ripserer_convert

include("extra_file.jl")

end


text:
module DuneIstlSolvers

using CxxWrap
using duneistl_julia_jll
@wrapmodule duneistl_julia_jll.libduneistl

function __init__()
    @initcxx
end

include("documentation.jl")

export DuneIstlSolver

export print_matrix, print_rhs,
       construct_matrix,
       add_value_matrix, add_value_rhs, add_value_x,
       get_value_matrix, get_value_rhs, get_value_x,
       reset_matrix, reset_rhs, reset_x,
       set_solver_type, set_target_reduction, set_max_iter, set_gmres_restart,
       set_preconditioner_type, set_preconditioner_relax, set_ilu_n, set_ilu_resort,
       get_solver_type, get_target_reduction, get_max_iter, get_gmres_restart,
       get_preconditioner_type, get_preconditioner_relax, get_ilu_n, get_ilu_resort,
       solve, set_verbose, get_verbose,
       get_solution, get_iterations, get_converged, get_reduction, get_conv_rate, get_elapsed

end


text:
function toh06a()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
   p = ones(21)
   z = [0.95557280578614 + 0.29475517441090im
        0.95557280578614 - 0.29475517441090im
        0.82623877431600 + 0.56332005806362im
        0.82623877431600 - 0.56332005806362im
        0.62348980185873 + 0.78183148246803im
        0.62348980185873 - 0.78183148246803im
        0.36534102436640 + 0.93087374864420im
        0.36534102436640 - 0.93087374864420im
        0.07473009358642 + 0.99720379718118im
        0.07473009358642 - 0.99720379718118im
        -0.98883082622513 + 0.14904226617617im
        -0.98883082622513 - 0.14904226617617im
        -0.90096886790242 + 0.43388373911756im
        -0.90096886790242 - 0.43388373911756im
        -0.73305187182983 + 0.68017273777092im
        -0.73305187182983 - 0.68017273777092im
        -0.50000000000000 + 0.86602540378444im
        -0.50000000000000 - 0.86602540378444im
        -0.22252093395631 + 0.97492791218182im
        -0.22252093395631 - 0.97492791218182im];
    z = [z ones(20)]

    p, PolyZeros(z)
end


text:
using DataStructures

struct Node
    children::Vector{Node}
    ID::Int64
    Node(ID::Int64) = new(Vector{Node}(), ID)
end

function DFS_recursive(n::Node)
    # Here we are doing something...
    println(n.ID)

    for child in n.children
        DFS_recursive(child)
    end
end

function DFS_recursive_postorder(n::Node)

    for child in n.children
        DFS_recursive_postorder(child)
    end

    # Here we are doing something...
    println(n.ID)
end

# This assumes only 2 children
function DFS_recursive_inorder_btree(n::Node)

    if (length(n.children) > 2)
        println("Not a binary tree!")
        exit(1)
    end

    if (length(n.children) > 0)
        DFS_recursive_inorder_btree(n.children[0])
        println(n.ID)
        DFS_recursive_inorder_btree(n.children[1])
    else
        println(n.ID)
    end
end

function DFS_stack(n::Node)
    s = Stack(Node)
    push!(s, n)

    while(length(s) > 0)
        println(top(s).ID)
        temp = pop!(s)
        for child in temp.children
            push!(s, child)
        end
    end
end

function BFS_queue(n::Node)
    q = Queue(Node)
    enqueue!(q, n)

    while(length(q) > 0)
        println(front(q).ID)
        temp = dequeue!(q)
        for child in temp.children
            enqueue!(q, child)
        end
    end
end

# function to create a simple, balanced tree
function create_tree(num_row::Int64, num_child::Int64)
    ret = Node(num_row)
    if (num_row == 0)
        return ret
    end

    for i = 1:num_child
        child = create_tree(num_row - 1, num_child)
        push!(ret.children, child)
    end

    return ret
end

function main()

    println("Creating Tree")
    root = create_tree(2,3)

    println("Using recursive DFS:")
    DFS_recursive(root);

    println("Using stack-based DFS:")
    DFS_stack(root);

    println("Using queue-based BFS:")
    BFS_queue(root);
end

main()


text:
function use_sixel(img::AbstractArray)
    encoder_backend[1] == :Sixel || return false

    # Sixel requires at least 6 pixels in row direction and thus doesn't perform very well for vectors.
    # ImageInTerminal encoder is good enough for vector case.
    ndims(img) == 1 && return false

    # Small images really do not need sixel encoding.
    # `100` is a randomly chosen value; it's not the best because
    # 100×100 image will be very small in terminal after sixel encoding.
    any(size(img) .<= 12) && return false
    all(size(img) .<= 100) && return false
    return true
end

"""
    imshow([stream], img, [depth::TermColorDepth], [maxsize])

Displays the given image `img` using unicode characters and
terminal colors (defaults to 256 colors).
`img` has to be an array of `Colorant`.

If working in the REPL, the function tries to choose the encoding
based on the current display size. The image will also be
downsampled to fit into the display (using `restrict`).
"""
function imshow(
        io::IO,
        img::AbstractArray{<:Colorant},
        colordepth::TermColorDepth,
        maxsize::Tuple = displaysize(io))
    use_sixel(img) && return sixel_encode(io, img)

    # otherwise, use our own implementation
    print_matrix(io, x) = imshow(io, x, colordepth, maxsize)
    Base.show_nd(io, img, print_matrix, true)
end

function imshow(
        io::IO,
        img::AbstractMatrix{<:Colorant},
        colordepth::TermColorDepth,
        maxsize::Tuple = displaysize(io))
    use_sixel(img) && return sixel_encode(io, img)

    # otherwise, use our own implementation
    io_h, io_w = maxsize
    img_h, img_w = map(length, axes(img))
    str = if img_h <= io_h-4 && 2img_w <= io_w
        first(encodeimg(BigBlocks(),   colordepth, img, io_h-4, io_w))
    else
        first(encodeimg(SmallBlocks(), colordepth, img, io_h-4, io_w))
    end
    for (idx, line) in enumerate(str)
        print(io, line)
        idx < length(str) && println(io)
    end
end

# colorant vector
function imshow(
        io::IO,
        img::AbstractVector{<:Colorant},
        colordepth::TermColorDepth,
        maxsize::Tuple = displaysize(io))
    @assert !use_sixel(img) "Sixel should be disabled for Vector colorant"

    io_h, io_w = maxsize
    img_w = length(img)
    str = if 3img_w <= io_w
        first(encodeimg(BigBlocks(),   colordepth, img, io_w))
    else
        first(encodeimg(SmallBlocks(), colordepth, img, io_w))
    end
    for (idx, line) in enumerate(str)
        print(io, line)
        idx < length(str) && println(io)
    end
end

imshow(io::IO, img, args...) = imshow(io, img, colormode[1], args...)
imshow(img, args...) = imshow(stdout, img, colormode[1], args...)
imshow(io::IO, img, colordepth::TermColorDepth, args...) = throw(ArgumentError("imshow only supports colorant arrays with 1 or 2 dimensions"))

"""
    imshow256([stream], img, [maxsize])

Displays the given image `img` using unicode characters and
the widely supported 256 terminal colors.
`img` has to be an array of `Colorant`.

If working in the REPL, the function tries to choose the encoding
based on the current display size. The image will also be
downsampled to fit into the display (using `restrict`).
"""
imshow256(io::IO, img, args...) = imshow(io, img, TermColor256(), args...)
imshow256(img, args...) = imshow256(stdout, img, args...)

"""
    imshow24bit([stream], img, [maxsize])

Displays the given image `img` using unicode characters and
the 24 terminal colors that some modern terminals support.
`img` has to be an array of `Colorant`.

If working in the REPL, the function tries to choose the encoding
based on the current display size. The image will also be
downsampled to fit into the display (using `restrict`).
"""
imshow24bit(io::IO, img, args...) = imshow(io, img, TermColor24bit(), args...)
imshow24bit(img, args...) = imshow24bit(stdout, img, args...)


text:
export build_SDP_Instance_from_SDPDual

function build_SDP_Instance_from_SDPDual(sdpdual::SDPDual)
    sdp_pb = SDP_Problem()

    ## First. Deal fully with moment matrices
    blocktovars = SortedDict{String, SortedSet}()
    n_sdp = 0
    for ((ctrname, cliquename), matrix) in sdpdual.constraints
        if ctrname == get_momentcstrname()
            n_sdp += 1
            blockname = cliquename
            block = SDP_Block(n_sdp, blockname)

            sdp_pb.name_to_sdpblock[blockname] = block
            sdp_pb.id_to_sdpblock[n_sdp] = block

            @assert matrix.matrixkind == :SDP
            blocktovars[blockname] = SortedSet{String}()
            for (α, β) in keys(matrix)
                var1 = format_string(α)
                var2 = format_string(β)

                push!(blocktovars[blockname], var1)
                push!(blocktovars[blockname], var2)
            end
        end
    end
    for (blockname, vars) in blocktovars
        block = sdp_pb.name_to_sdpblock[blockname]
        n=0
        for var in vars
            n += 1
            block.var_to_id[var] = n
        end
    end

    # warn("Done with moment constraints.\n")

    ## Second. Filling obj and ctrs coefficients
    ## Dealing with objective
    ctr_name = ("1", "1", "objective")
    for (moment, fαβ) in sdpdual.objective
        block_name, α, β = split_moment(moment)

        α_str, β_str = format_string(α), format_string(β)
        var1, var2 = max(α_str, β_str), min(α_str, β_str)

        @assert haskey(sdp_pb.name_to_sdpblock, block_name)
        @assert haskey(sdp_pb.name_to_sdpblock[block_name].var_to_id, var1)
        @assert haskey(sdp_pb.name_to_sdpblock[block_name].var_to_id, var2)

        if product(α, β) == Exponent()
            !haskey(sdp_pb.cst_ctr, ctr_name) && (sdp_pb.cst_ctr[ctr_name] = 0)

            sdp_pb.cst_ctr[ctr_name] += fαβ
        else
            sdp_pb.matrices[(ctr_name, block_name, var1, var2)] = fαβ * (var1!=var2 ? 0.5 : 1)
        end
    end

    ## Dealing with constraints
    for ((ctrname, cliquename), matrix) in sdpdual.constraints

        if ctrname != get_momentcstrname()

            # Name of auxiliary SDP constrained matrix
            S_name = get_auxSDPmatrix_name(ctrname, cliquename)

            # for all matrix coefficients
            # Deal with S[γ, δ]
            for ((γ, δ), momentpoly) in matrix
                ctr_name = get_dual_ctr_name(ctrname, cliquename, γ, δ)
                γ_str = format_string(γ)
                δ_str = format_string(δ)

                for (moment, fαβ) in momentpoly
                    block_name, α, β = split_moment(moment)

                    α_str, β_str = format_string(α), format_string(β)
                    var1, var2 = max(α_str, β_str), min(α_str, β_str)

                    @assert haskey(sdp_pb.name_to_sdpblock, block_name)
                    @assert haskey(sdp_pb.name_to_sdpblock[block_name].var_to_id, var1)
                    @assert haskey(sdp_pb.name_to_sdpblock[block_name].var_to_id, var2)

                    # @show (ctr_name, block_name, var1, var2), fαβ
                    if product(α, β) == Exponent()
                        !haskey(sdp_pb.cst_ctr, ctr_name) && (sdp_pb.cst_ctr[ctr_name] = 0)

                        sdp_pb.cst_ctr[ctr_name] += fαβ
                    else
                        sdp_pb.matrices[(ctr_name, block_name, var1, var2)] = fαβ * (var1!=var2 ? 0.5 : 1)
                    end
                end

                # If the matrix is SDP, add link to auxiliary SDP var
                if matrix.matrixkind == :SDP
                    block_name = S_name
                    var1 = max(γ_str, δ_str)
                    var2 = min(γ_str, δ_str)

                    sdp_pb.matrices[(ctr_name, block_name, var1, var2)] = -1 * (var1!=var2 ? 0.5 : 1)
                else
                    # Nothing to do, constraint will be scalar, ==0 by default.
                end
            end
        end
    end

    ## Third. Add auxiliary SDP matrices standing for SDP localizing constraints
    # println("--------------------------------------------------------")
    momentvars = SortedSet(collect(keys(sdp_pb.name_to_sdpblock)))

    for block_name in momentvars
        ctr_name = (block_name*"_1ctr", "1", "1")
        sdp_pb.matrices[(ctr_name, block_name, "1", "1")] = 1
        sdp_pb.cst_ctr[ctr_name] = -1
    end

    # @show momentvars
    # println()
    empty!(blocktovars)

    n_sdp = length(sdp_pb.name_to_sdpblock)
    for ((ctr_name, blockname, var1, var2), f_αβ) in sdp_pb.matrices
        if blockname ∉ momentvars
            if !haskey(sdp_pb.name_to_sdpblock, blockname)
                n_sdp += 1
                block = SDP_Block(n_sdp, blockname)

                sdp_pb.name_to_sdpblock[blockname] = block
                sdp_pb.id_to_sdpblock[n_sdp] = block
            else
                block = sdp_pb.name_to_sdpblock[blockname]
            end

            !haskey(blocktovars, blockname) && (blocktovars[blockname] = SortedSet{String}())
            push!(blocktovars[blockname], var1)
            push!(blocktovars[blockname], var2)
        end
    end
    for (blockname, vars) in blocktovars
        block = sdp_pb.name_to_sdpblock[blockname]
        n=0
        for var in vars
            n += 1
            block.var_to_id[var] = n
        end
    end

    ## Deal with equality constaints
    if length(sdpdual.moments_overlap) != 0
        for (expo, cliques) in sdpdual.moments_overlap
            if expo != Exponent()
                @assert length(cliques) > 1

                clique_ref = first(cliques)
                moment_ref = Moment(expo, clique_ref)

                block_name, α, β = split_moment(moment_ref)
                α_str, β_str = format_string(α), format_string(β)
                var1, var2 = max(α_str, β_str), min(α_str, β_str)

                for clique in setdiff(cliques, Set([clique_ref]))
                    ctr_name = get_coupling_ctr_name(clique_ref, clique, var1, var2)

                    sdp_pb.matrices[(ctr_name, clique_ref, var1, var2)] = 1
                    sdp_pb.matrices[(ctr_name, clique, var1, var2)] = -1
                end
            end
        end
    end

    # Build structural information and string to id maps
    # set_blockvartypes!(sdp_pb)

    # set_blocks!(sdp_pb, sdpdual)

    # MathProgComplex.set_blocks!(sdp_pb)

    MathProgComplex.set_constraints!(sdp_pb)
    MathProgComplex.set_scalvars!(sdp_pb)

    return sdp_pb
end

function get_dual_ctr_name(ctrname, cliquename, γ, δ)
    return replace(ctrname, " ", "")*"_"*replace(cliquename, " ", ""), format_string(γ), format_string(δ)
end

function get_auxSDPmatrix_name(ctrname, cliquename)
    return "S_"*replace(ctrname, " ", "")*"_"*replace(cliquename, " ", "")
end

function get_coupling_ctr_name(clique_ref, clique, var1, var2)
    return var1, var2, clique_ref*"_"*clique
end

"""
    α, β = split_expo(expo::Exponent)

Split the exponent into two exponents of conjugated and explicit variables in the complex case.
Real case is not supported yet.
"""
function split_moment(moment::Moment)
    α, β = Exponent(), Exponent()

    if moment.conj_part != Exponent()
        # Complex moment, separation is trivial
        α = moment.conj_part
        β = moment.expl_part

    else
        # Real moment, more complicated...
        explsum, conjsum = MathProgComplex.get_sumdegs(moment.expl_part)
        @assert conjsum == 0

        α_deg = β_deg = 0
        α_curadd = β_curadd = 0
        for (var, degree) in moment.expl_part
            @assert degree.conjvar == 0

            if α_deg < β_deg
                α_curadd = ceil(degree.explvar / 2)
                β_curadd = floor(degree.explvar / 2)
            else
                α_curadd = floor(degree.explvar / 2)
                β_curadd = ceil(degree.explvar / 2)
            end

            product!(α, Exponent(SortedDict(var=>Degree(α_curadd, 0))))
            α_deg += α_curadd
            product!(β, Exponent(SortedDict(var=>Degree(β_curadd, 0))))
            β_deg += β_curadd
        end

        @assert product(α, β) == product(moment.expl_part, moment.conj_part)
        αexplsum, αconjsum = get_sumdegs(α)
        βexplsum, βconjsum = get_sumdegs(β)

        @assert min(αexplsum, βexplsum) == floor(explsum/2)
        @assert max(αexplsum, βexplsum) == ceil(explsum/2)
    end

    return moment.clique, α, β
end

"""
set_blockvartypes!(sdp_pb::SDP_Problem)

# Set attributes `name_to_sdpblock` and `id_to_sdpblock` given a primal sdp.
"""
function set_blockvartypes!(sdp_pb::SDP_Problem)
    n_sdp = 0

    blocknames = SortedSet([k[2] for k in keys(sdp_pb.matrices)])

    for blockname in blocknames
        n_sdp += 1

        block = SDP_Block(n_sdp, blockname)
        sdp_pb.name_to_sdpblock[blockname] = block
        sdp_pb.id_to_sdpblock[n_sdp] = block
    end
end

function set_blocks!(sdp_pb::SDP_Problem, sdpdual::SDPDual)
    # Collecting moments from each moment constraint
    for ((ctrname, cliquename), matrix) in sdpdual.constraints

        # moment constraint
        if ctrname == get_momentcstrname()
            blockname = cliquename
            block = sdp_pb.name_to_sdpblock[blockname]

            for ((α, β), momentpoly) in matrix
                @assert length(momentpoly) == 1

                moment = first(momentpoly)[1]  #y_αβ
                # @show moment
                α_str = format_string(moment.conj_part)
                β_str = format_string(moment.expl_part)

                if !haskey(block.var_to_id, α_str)
                    block.var_to_id[α_str] = length(block.var_to_id) + 1
                end
                if !haskey(block.var_to_id, β_str)
                    block.var_to_id[β_str] = length(block.var_to_id) + 1
                end
            end

        end
    end
end

###########################################################################
function get_formatedexpo(expo::Exponent)
    if first(get_sumdegs(expo)) == 2
        if length(expo) == 1
            var1 = string(first(keys(expo)))
            var2 = string(first(keys(expo)))
        else
            var1, var2 = string.(collect(keys(expo)))
        end
        linetype = "QUAD"
    elseif first(get_sumdegs(expo)) == 1
        var1 = string(first(keys(expo)))
        var2 = "NONE"
        linetype = "LIN"
    elseif first(get_sumdegs(expo)) == 0
        var1 = "NONE"
        var2 = "NONE"
        linetype = "CONST"
    end

    return var1, var2, linetype
end

function print_QCQPform(io::IO, momentrelax::SDPDual{T}) where T
    maxvarlen = 15
    maxcstrlen = 15

    println(io, "Moment Relaxation Problem, QCQP form:")
    for moment in sort(collect(keys(momentrelax.objective)))
        expo = product(moment.expl_part, moment.conj_part)

        var1, var2, linetype = get_formatedexpo(expo)

        cstrname = "OBJ"
        val1 = momentrelax.objective[moment]
        val2 = 0
        MathProgComplex.print_dat_line(io, linetype, cstrname, var1, var2, val1, val2, maxvarlen, maxcstrlen)
    end

    for (cstrname, blocname) in sort(collect(keys(momentrelax.constraints)))

        if !ismatch(r"_lo", cstrname)
            mmtmat = momentrelax.constraints[(cstrname, blocname)]

            for (key, momentpoly) in mmtmat
                for (moment, coeff) in momentpoly
                    expo = product(moment.expl_part, moment.conj_part)

                    var1, var2, linetype = get_formatedexpo(expo)

                    val1 = coeff
                    val2 = 0
                    MathProgComplex.print_dat_line(io, linetype, cstrname, var1, var2, val1, val2, maxvarlen, maxcstrlen)
                end
            end
        end

    end

    # println(io, "→ Moments clique overlap:")
    # if length(momentrelax.moments_overlap) > 0
    #     mmtlength = maximum(x->length(string(x)), keys(momentrelax.moments_overlap))
    #     for moment in sort(collect(keys(momentrelax.moments_overlap)))
    #         cliquenames = momentrelax.moments_overlap[moment]
    #         print(io, " → ")
    #         print_string(io, string(moment), mmtlength)
    #         for clique in sort(collect(cliquenames)) print(io, "$clique, ") end
    #         @printf(io, "\b\b \n")
    #     end
    # else
    #     print(io, "  None")
    # end
end


text:
function ismongotype(value)
    unit_type = Union{AbstractString,Nothing,Number,Bool}
    typeof(value)<:unit_type && return true
    if typeof(value)<:AbstractVector
        typeof(value)<:AbstractVector{T} where T<:unit_type && return true
        return all([ismongotype(v) for v in value])
    elseif typeof(value)<:AbstractDict
        return all(vcat([[typeof(v)<:AbstractString, ismongotype(k)] for (v, k) in value]...))
    end
    return false
end

function Base.setindex!(bson::Mongoc.BSON, uuid::UUID, index::AbstractString)
    bson[index] = string(uuid)
end

"""
    push_instance!(collection::Mongoc.Collection,
                   instance_name::String,
                   paths::Union{T, AbstractVector{T}}=Mongoc.BSON(),
                   features::AbstractDict=Mongoc.BSON()) where T<:AbstractDict

Insert an instance in `collection`.
"""
function push_instance!(collection::Mongoc.Collection,
                        instance_name::String,
                        paths::Union{T, AbstractVector{T}}=Mongoc.BSON(),
                        features::AbstractDict=Mongoc.BSON()) where T<:AbstractDict
    if !isinstance(collection, instance_name)
        document = Mongoc.BSON("_id" => instance_name,
                               "paths" => paths,
                               "decompositions" => [],
                               "features" => features,
                               "insertion_date" => Dates.now())
        push!(collection, document)
        return true
    else
        return false
    end
end

"""
    push_decomposition!(collection::Mongoc.Collection,
                        instance_name::String,
                        added_edges::Vector{Vector{Int}},
                        cliques::Vector{Vector{Int}},
                        cliquetree::Vector{Vector{Int}},
                        options_src::Dict,
                        options_dst::Dict,
                        path_MOSEK::Union{Nothing,String}=nothing,
                        features::AbstractDict=Mongoc.BSON())

Insert a decomposition in `collection`.
"""
function push_decomposition!(collection::Mongoc.Collection,
                             instance_name::String,
                             #cliques::AbstractVector,
                             #cliquetree::AbstractVector,
                             #options_src::Dict=Dict(),
                             #options_dst::Dict=Dict(),
                             #nb_added_edges::Int=-1,
                             #path_MOSEK::Union{Nothing,String}=nothing,
                             data::AbstractDict,
                             features::AbstractDict=Mongoc.BSON())
    rng = MersenneTwister(75321)
    oid = string(uuid1(rng))
    data["cliques"] = map(sort, data["cliques"])
    permutation = sortperm(data["cliques"])
    permutation_inv = [findall(x -> x == i, permutation)[1] for i in 1:length(data["cliques"])]
    data["cliques"] = data["cliques"][permutation]
    data["cliquetree"] = map(edge -> map(v -> permutation_inv[v], edge), data["cliquetree"])
    data["seed"] = string(Random.default_rng())
    document = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid),
                           "features" => features)
    data = filter(d -> ismongotype(d.second), data)
    merge!(document, data)
    push!(collection, document)
    return instance_name, oid
end

"""
    add_decomposition_in_instance!(collection::Mongoc.Collection, instance_name::String, added_edges::AbstractVector{N}) where N<:AbstractVector{Int}

Add `added_edges` to the decompositions field.
"""
function add_decomposition_in_instance!(collection::Mongoc.Collection,
                                        instance_name::String,
                                        oid::Union{UUID,AbstractString})
    selector = Mongoc.BSON("_id" => instance_name)
    update = Mongoc.BSON("\$addToSet" => Mongoc.BSON("decompositions" => oid))
    Mongoc.update_one(collection, selector, update)
end

function add_merge_in_decomposition!(collection::Mongoc.Collection,
                                     instance_name::String,
                                     oid::Union{UUID,AbstractString},
                                     merge_oid::Union{UUID,AbstractString},
                                     merge_type,
                                     merge_treshold_type,
                                     merge_treshold)
    select = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    merge_infos = Mongoc.BSON("oid" => merge_oid,
                              "merge_type" => merge_type,
                              "merge_treshold_type" => merge_treshold_type,
                              "merge_treshold" => merge_treshold)
    update = Mongoc.BSON("\$addToSet" => Mongoc.BSON("merge" => merge_infos))
    Mongoc.update_one(collection, select, update)
end

function setfields!(collection, selector, fields)
    update = Mongoc.BSON("\$set" => fields)
    Mongoc.update_one(collection, selector, update)
end

function setfields_decomposition!(collection::Mongoc.Collection,
                                 instance_name::String,
                                 oid::Union{UUID,AbstractString},
                                 fields::AbstractDict)
    selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    setfields!(collection, selector, fields)
end

"""
    setfeatures!(collection::Mongoc.Collection, selector::AbstractDict, features::AbstractDict)

Add/Update the field "features" in the document corresponding to `selector` in the collection `collection`.
"""
function setfeatures!(collection::Mongoc.Collection,
                      selector::AbstractDict,
                      features::AbstractDict)
    formated_dict_features = Dict("features.$(first(item))" => last(item) for item in features)
    update = Mongoc.BSON("\$set" => formated_dict_features)
    Mongoc.update_one(collection, selector, update)
end

"""
    setfeatures_instance!(collection::Mongoc.Collection, instance_name::String, features::AbstractDict)

Add/Update/Replace the field "features" in the document corresponding to "_id" == `instance_name`.
"""
function setfeatures_instance!(collection::Mongoc.Collection,
                               instance_name::String,
                               features::AbstractDict)
    selector = Mongoc.BSON("_id" => instance_name)
    setfeatures!(collection, selector, features)
end

"""
    setfeatures_decomposition!(collection::Mongoc.Collection, instance_name::String, added_edges::Vector{Vector{Int}}, features::AbstractDict)

Add/Update/Replace the field "features" in the document corresponding to "_id" == {"instance_name" : `instance_name`, "added_edges" : `added_edges`}.
"""
function setfeatures_decomposition!(collection::Mongoc.Collection,
                                    instance_name::String,
                                    oid::Union{UUID,AbstractString},
                                    features::AbstractDict)
    selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    setfeatures!(collection, selector, features)
end

function setMPIinfos(collection::Mongoc.Collection,
                     instance_name::String,
                     oid::Union{UUID,AbstractString},
                     mpi_infos::AbstractDict)
    formated_dict_mpi = Dict("mpi.$(first(item))" => last(item) for item in mpi_infos)
    selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    update = Mongoc.BSON("\$set" => formated_dict_mpi)
    Mongoc.update_one(collection, selector, update)
end

function delete_decomposition(collection::Mongoc.Collection,
                              instance_name::String,
                              oid::Union{UUID,AbstractString})
    selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name, "oid" => oid))
    Mongoc.delete_one(collection, selector)
end

function delete_similar_decomposition(collection::Mongoc.Collection; verbose=false)
    decompositions = collect(collection)
    for i in 1:length(collection) - 1
        for j in i + 1:length(collection)
            document1 = decompositions[i]
            document2 = decompositions[j]
            # The clique array is sorted when inserted so no need to sort them before the comparison
            if document1["_id"]["instance_name"] == document2["_id"]["instance_name"] && document1["cliques"] == document2["cliques"]
                # If one of the decomposition is the default decomposition
                if document1["nb_added_edges"] == 0
                    verbose && @info "Deleting $(document2["_id"]["instance_name"]), $(document2["_id"]["oid"])"
                    delete_decomposition(collection, document2["_id"]["instance_name"], document2["_id"]["oid"])
                    continue
                elseif document2["nb_added_edges"] == 0
                    verbose && @info "Deleting $(document1["_id"]["instance_name"]), $(document1["_id"]["oid"])"
                    delete_decomposition(collection, document1["_id"]["instance_name"], document1["_id"]["oid"])
                    continue
                end
                # If one of the decomposition is already solved, we keep it
                issolved1 = haskey(document1["features"], "solver") && haskey(document1["features"]["solver"], "solving_time")
                issolved2 = haskey(document2["features"], "solver") && haskey(document2["features"]["solver"], "solving_time")
                if issolved1 && !issolved2
                    verbose && @info "Deleting $(document2["_id"]["instance_name"]), $(document2["_id"]["oid"])"
                    delete_decomposition(collection, document2["_id"]["instance_name"], document2["_id"]["oid"])
                elseif !issolved1 && issolved2
                    verbose && @info "Deleting $(document1["_id"]["instance_name"]), $(document1["_id"]["oid"])"
                    delete_decomposition(collection, document1["_id"]["instance_name"], document1["_id"]["oid"])
                else
                    verbose && @info "Deleting $(document2["_id"]["instance_name"]), $(document2["_id"]["oid"])"
                    delete_decomposition(collection, document2["_id"]["instance_name"], document2["_id"]["oid"])
                end
            end
        end
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
function NewKrest(r)
    for side in 0:3
        GoToTheWallAndMark(r, HorizonSide(side))
    end
end

function GoToTheWallAndMark(r::Robot, side::HorizonSide)
    x = 0
    while !mainborder(r,side)
        putmarker!(r)
        trymove!(r, side)
        x = x + 1
    end
    putmarker!(r)
    while x != 0
        trymove!(r, HorizonSide((Int(side) + 2) % 4))
        x = x - 1
    end
end

function trymove!(r::Robot, side::HorizonSide)
    x = 0
    if isborder(r,side)
        while isborder(r, side)
            move!(r,HorizonSide((Int(side) + 1) % 4))
            x = x + 1
        end
        move!(r,side)
        while isborder(r,HorizonSide((Int(side) + 3) % 4))
            move!(r,side)
        end
        while x != 0
            move!(r,HorizonSide((Int(side) + 3) % 4))
            x = x - 1
        end
    else
        move!(r,side)
    end
end

function mainborder(r,side)
    x = 0
    b = true
    while isborder(r,side) && !isborder(r,HorizonSide((Int(side) + 3) % 4))
        move!(r,HorizonSide((Int(side) + 3) % 4))
        x = x + 1
    end
    if isborder(r,side)
        b = true
    else
        b = false
    end
    while x != 0
        move!(r,HorizonSide((Int(side) + 1) % 4))
        x = x - 1
    end
    return b
end


text:
"""
NGSIMTrajdata

The trajectory data stored in the original NGSIM dataset format.
The dataset is a white-space separated text file with columns:
    id                  - I - Vehicle identification number (ascending by time of entry into section)
    frame               - I - Frame Identification number (ascending by start time), units 1/10 of a second
    n_frames_in_dataset - I - Total number of frames in which the vehicle appears in this data set, units 1/10 of a second
    epoch               - I - Elapsed time since Jan 1, 1970, in milliseconds
    local_x             - F - Lateral (X) coordinate of the front center of the vehicle with respect to the left-most edge of the section in the direction of travel, in feet
    local_y             - F - Longitudinal (Y) coordinate of the front center of the vehicle with respect to the entry edge of the section in the direction of travel, in feet
    global_x            - F - X Coordinate of the front center of the vehicle based on CA State Plane III in NAD83
    global_y            - F - Y Coordinate of the front center of the vehicle based on CA State Plane III in NAD83
    length              - F - Length of the vehicle, in feet
    width               - F - Width of the vehicle, in feet
    class               - I - vehicle class, 1 - motorcycle, 2 - auto, 3 - truck
    speed               - F - Instantaneous velocity of vehicle, in ft/second
    acc                 - F - Instantaneous acceleration of vehicle, in ft/second^2
    lane                - I - Current lane position of vehicle
    carind_front        - I - Vehicle Id of the lead vehicle in the same lane. A value of '0' represents no preceding vehicle
    carind_rear         - I - Vehicle Id of the vehicle following the subject vehicle in the same lane. A value of '0' represents no following vehicle
    dist_headway        - F - Spacing provides the distance between the front-center of a vehicle to the front-center of the preceding vehicle, in feet
    time_headway        - F - Headway provides the time to travel from the front-center of a vehicle (at the speed of the vehicle) to the front-center of the preceding vehicle. A headway value of 9999.99 means that the vehicle is traveling at zero speed (congested conditions), in second
"""
mutable struct NGSIMTrajdata
    df         :: DataFrame
    car2start  :: Dict{Int, Int}         # maps carindex to starting index in the df
    frame2cars :: Dict{Int, Vector{Int}} # maps frame to list of carids in the scene

    function NGSIMTrajdata(input_path::String)

        @assert(isfile(input_path))

        df = readtable(input_path, separator=' ', header = false)
        col_names = [:id, :frame, :n_frames_in_dataset, :epoch, :local_x, :local_y, :global_x, :global_y, :length, :width, :class, :speed, :acc, :lane, :carind_front, :carind_rear, :dist_headway, :time_headway]
        for (i,name) in enumerate(col_names)
            rename!(df, Symbol(@sprintf("x%d", i)), name)
        end

        df[:global_heading] = fill(NaN, nrow(df))

        car2start = Dict{Int, Int}()
        frame2cars = Dict{Int, Vector{Int}}()

        for (dfind, carid) in enumerate(df[:id])
            if !haskey(car2start, carid)
                car2start[carid] = dfind
            end

            frame = convert(Int, df[dfind, :frame])
            if !haskey(frame2cars, frame)
                frame2cars[frame] = [carid]
            else
                frame2cars[frame] = push!(frame2cars[frame], carid)
            end
        end

        new(df, car2start, frame2cars)
    end
end

Records.nframes(trajdata::NGSIMTrajdata) = maximum(keys(trajdata.frame2cars))
carsinframe(trajdata::NGSIMTrajdata, frame::Int) = get(trajdata.frame2cars, frame, Int[]) # NOTE: memory allocation!
carid_set(trajdata::NGSIMTrajdata) = Set(keys(trajdata.car2start)) # NOTE: memory allocation!
nth_carid(trajdata::NGSIMTrajdata, frame::Int, n::Int) = trajdata.frame2cars[frame][n]
first_carid(trajdata::NGSIMTrajdata, frame::Int) = nth_carid(trajdata, frame, 1)
iscarinframe(trajdata::NGSIMTrajdata, carid::Int, frame::Int) = in(carid, carsinframe(trajdata, frame))

function car_df_index(trajdata::NGSIMTrajdata, carid::Int, frame::Int)
    #=
    given frame and carid, find index of car in trajdata
    Returns 0 if it does not exist
    =#

    df = trajdata.df

    lo = trajdata.car2start[carid]
    framestart = df[lo, :frame]

    retval = 0

    if framestart == frame
        retval = lo
    elseif frame ≥ framestart
        retval = frame - framestart + lo
        n_frames = df[lo, :n_frames_in_dataset]
        if retval > lo + n_frames
            retval = 0
        end
    end

    retval
end
function get_frame_range(trajdata::NGSIMTrajdata, carid::Int)
    lo = trajdata.car2start[carid]
    framestart = trajdata.df[lo, :frame]

    n_frames = trajdata.df[lo, :n_frames_in_dataset]
    frameend = framestart + n_frames - 1

    framestart:frameend
end

function pull_vehicle_headings!(trajdata::NGSIMTrajdata;
    v_cutoff::Float64 = 2.5, # speeds below this will use a linearly interpolated heading [fps]
    smoothing_width::Float64 = 0.5, # [s]
    )

    df = trajdata.df

    for carid in carid_set(trajdata)

        frames = collect(get_frame_range(trajdata, carid))
        states = get_state_list_global(trajdata, carid, frames)

        arr_x = map(p->p.x, states)-states[1].x
        arr_y = map(p->p.y, states)-states[1].y

        arr_v = map(i->hypot(states[i]-states[i-1]), 2:length(states)) / 0.1
        slow_indeces = find(v->v < v_cutoff, arr_v)

        slow_segments = Tuple{Int,Int}[]
        i = 0
        while i < length(slow_indeces)
            i += 1
            i_lo = i
            while i < length(slow_indeces) && slow_indeces[i+1] == slow_indeces[i]+1
                i += 1
            end
            i_hi = i
            push!(slow_segments, (slow_indeces[i_lo], slow_indeces[i_hi]))
        end
        slow_segments

        arr_heading = map(i->atan2(convert(VecE2, states[i]-states[i-1])), 2:length(states))
        unshift!(arr_heading, arr_heading[1])

        arr_x_smoothed = symmetric_exponential_moving_average(arr_x, smoothing_width)
        arr_y_smoothed = symmetric_exponential_moving_average(arr_y, smoothing_width)
        arr_dx_smoothed = arr_x_smoothed[2:end] - arr_x_smoothed[1:end-1]
        arr_dy_smoothed = arr_y_smoothed[2:end] - arr_y_smoothed[1:end-1]
        arr_heading2 = map(i->atan2(arr_y_smoothed[i]-arr_y_smoothed[i-1], arr_x_smoothed[i]-arr_x_smoothed[i-1]), 2:length(states))
        unshift!(arr_heading2, arr_heading2[1])

        arr_heading3 = deepcopy(arr_heading2)
        for (lo,hi) in slow_segments
            heading_lo = arr_heading2[max(lo-1, 1)]
            heading_hi = arr_heading2[min(hi+1, length(arr_heading2))]
            for f in lo : hi
                arr_heading3[f] = heading_lo + (heading_hi - heading_lo)*(f-lo)/(hi-lo+1)
            end
        end

        for (i,frame) in enumerate(frames)
            df_ind = car_df_index(trajdata, carid, frame)
            df[df_ind, :global_heading] = arr_heading3[i]
        end
    end

    @assert(findfirst(v->isnan(v), df[:global_heading]) == 0)

    trajdata
end

const NGSIM_TRAJDATA_PATHS = [
                        Pkg.dir("NGSIM", "data", "i101_trajectories-0750am-0805am.txt"),
                        Pkg.dir("NGSIM", "data", "i101_trajectories-0805am-0820am.txt"),
                        Pkg.dir("NGSIM", "data", "i101_trajectories-0820am-0835am.txt"),
                        Pkg.dir("NGSIM", "data", "i80_trajectories-0400-0415.txt"),
                        Pkg.dir("NGSIM", "data", "i80_trajectories-0500-0515.txt"),
                        Pkg.dir("NGSIM", "data", "i80_trajectories-0515-0530.txt"),
                       ]

load_ngsim_trajdata(i::Int) = NGSIMTrajdata(NGSIM_TRAJDATA_PATHS[i])


text:
# PyPlot: needs Python and matplotlib installed
# [matplotlib.pyplot docs](http://matplotlib.org/api/pyplot_api.html)
import PyPlot
# plot 5 random numbers in [0,1], PyPlot.plot creates a new figure
PyPlot.plot(rand(5))
# labeling the axes, creating a title:
PyPlot.xlabel("x-axis")
PyPlot.ylabel("y-axis")
PyPlot.title("Random")
# creating the legend:
l = ["random data"]
# legend takes an array of strings, but since strings can be indexed to get
# chars (e.g. l[1] == 'r') the brackets are necessary
PyPlot.legend(l)

# create a new figure:
f = PyPlot.figure(2)
# generating data:
y1 = randn(10)
y2 = randn(10)
x = 1:10 # PyPlot can handle ranges as well as arrays
# setup a grid for plotting:
PyPlot.subplot2grid((3, 1), (0, 0), rowspan=2)
# options for plotting are passed as keyword-arguments:
PyPlot.plot(x, y1, c="red", marker=".", linestyle="None")
# a second plot will plot in an existing figure/subplot
PyPlot.plot(x, y2)
# add text to the plot:
PyPlot.text(x[4]+.1, y1[4], "value #4", verticalalignment="center")
# PyPlot will adjust the axes automatically, but xlim and ylim can explicitly
# change the smallest and largest value displayed on either axis
PyPlot.xlim([0, 11])
PyPlot.legend(["dots", "line"], loc="upper left")

# change into the other subplot:
PyPlot.subplot2grid((3, 1), (2, 0), rowspan=2)
# plot some more data into the lower subplot
PyPlot.plot(x, y1-y2)
PyPlot.xlim([0, 11])
PyPlot.ylim([-3, 3])
# change the ticks used on the y-axis
PyPlot.yticks([-3, 0, 3])
# save the current figure
PyPlot.savefig("pyplot.png")

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
using RiskAdjustedLinearizations, ForwardDiff, Random, Test

Random.seed!(1793)

aR = rand(2)
bR = rand(3)
aD = rand(2)
bD = rand(3)

AR = rand(2, 2, 3)
BR = rand(3, 2, 3, 2)
AD = ForwardDiff.Dual.(rand(2, 2, 3))
BD = ForwardDiff.Dual.(rand(3, 2, 3, 2))

@testset "dualarray and dualvector" begin
    @test length(RiskAdjustedLinearizations.dualarray(AR, BR)) == length(AR)
    @test length(RiskAdjustedLinearizations.dualarray(BR, AR)) == length(BR)
    @test length(RiskAdjustedLinearizations.dualarray(AR, BD)) == length(similar(AR, eltype(BD)))
    @test length(RiskAdjustedLinearizations.dualarray(BD, AR)) == length(similar(BD, eltype(AR)))
    @test length(RiskAdjustedLinearizations.dualarray(AD, BD)) == length(AD)
    @test length(RiskAdjustedLinearizations.dualarray(BD, AD)) == length(BD)
    @test length(RiskAdjustedLinearizations.dualarray(AD, BR)) == length(AD)
    @test length(RiskAdjustedLinearizations.dualarray(BR, AD)) == length(similar(BR, eltype(AD)))

    @test length(RiskAdjustedLinearizations.dualvector(aR, bR)) == length(aR)
    @test length(RiskAdjustedLinearizations.dualvector(bR, aR)) == length(bR)
    @test length(RiskAdjustedLinearizations.dualvector(aR, bD)) == length(similar(aR, eltype(bD)))
    @test length(RiskAdjustedLinearizations.dualvector(bD, aR)) == length(similar(bD, eltype(aR)))
    @test length(RiskAdjustedLinearizations.dualvector(aD, bD)) == length(aD)
    @test length(RiskAdjustedLinearizations.dualvector(bD, aD)) == length(bD)
    @test length(RiskAdjustedLinearizations.dualvector(aD, bR)) == length(aD)
    @test length(RiskAdjustedLinearizations.dualvector(bR, aD)) == length(similar(bR, eltype(aD)))
end

nothing


text:
using Revise, ApproxOperator

elements, nodes = importmsh("./msh/bar.msh")
nₚ = length(nodes[:x])
nₑ = length(elements["Domain"])

set𝓖!(elements["Domain"],:SegGI2)
set𝓖!(elements["NBC"],:PoiGI1)
set𝓖!(elements["EBC"],:PoiGI1)
elements["EBC"] = Element{:Seg2}(elements["Domain"][1],elements["EBC"][1])

r = 3
prescribe!(elements["Domain"],:u,(x,y,z)->x^r)
prescribe!(elements["Domain"],:∂u∂x,(x,y,z)->r*x^abs(r-1))
prescribe!(elements["Domain"],:b,(x,y,z)->-r*(r-1)*x^abs(r-2))
prescribe!(elements["NBC"],:t,(x,y,z)->r*x^abs(r-1))
prescribe!(elements["EBC"],:g,(x,y,z)->x^r)

coefficient = (:k=>1.0,:α=>1e3)
ops = [Operator(:∫∇v∇uvbdΩ,coefficient...),
       Operator(:∫vtdΓ,coefficient...),
       Operator(:∫∇𝑛vgdΓ,coefficient...),
       Operator(:H₁)]

k = zeros(nₚ,nₚ)
f = zeros(nₚ)

ops[1](elements["Domain"],k,f)
ops[2](elements["NBC"],f)
ops[3](elements["EBC"],k,f)
d = k\f
push!(nodes,:d=>d)
h1, l2 = ops[4](elements["Domain"])

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
#= STRUCTS =#
abstract type AbstractRun end

mutable struct Run <: AbstractRun
    name::AbstractString
    populations::Array{Any}
    df::DataFrame
    parent::Union{Run, Nothing}
    children::Array{Run}
    function Run(populations::Array{Any}; records = CURRENT_RECORDS, name = missing)
        # Assign name or random hash
        name = (name !== missing) ? name : randstring(5)

        # Create results table
        if length(populations) !== 0
            df = vcat([DataFrame(population) for population in populations]...)
            df.REP = repeat(1:length(populations), inner = nrow(df) ÷ length(populations))
        else
            df = DataFrame()
        end

        # Write CSV
        if (name !== "")
            file = joinpath(records.dir, string(name) * ".csv")
            CSV.write(file, df)
        end

        # Return Run object
        return new(
            name,
            populations,
            df,
            nothing,
            Run[]
        )
    end
end

Base.show(io::IO, run::Run) = begin
    println(io, "Run: $(run.name)")
    println(io, "  parent: $(run.parent === nothing ? nothing : run.parent.name)")
    println(io, "  children: $(length(run.children))")
    println(io, "  runs: $(length(run.populations))")
end

mutable struct Records
    dir::AbstractString
    runs::AbstractString
    models::AbstractString

    size::Int
    sentinel::Run
    head::Run

    function Records(dir, runs, models)
        dummy = Run([]; name = "")
        @assert !isdir(dir) "Records directory already exists"
        @assert !isdir(joinpath(dir, runs)) "Runs directory already exists"
        @assert !isdir(joinpath(dir, models)) "Models directory already exists"
        mkdir(dir)
        mkdir(joinpath(dir, runs))
        mkdir(joinpath(dir, models))
        return new(dir, runs, models, 0, dummy, dummy)
    end
end

Base.show(io::IO, records::Records) = begin
    println(io, "Records")
    println(io, "  directory: $(records.dir)")
    print(io, "  size: $(records.size)")
end


text:
function gauleg(n::Int, xmin::FT, xmax::FT; norm=false) where FT
    ξ, w = gausslegendre(n)
    ξ = (xmax - xmin) / FT(2) * ξ .+ (xmin + xmax) / FT(2)
    norm ? w /= sum(w) : w *= (xmax - xmin) / FT(2)
    return FT.(ξ), FT.(w)
end


text:
# Based on https://math.stackexchange.com/questions/3176543/intersection-point-of-2-lines-defined-by-2-points-each/3176590
function intersection(e1::E, e2::E, tol = 1e-15)
    nx = e1.q.x - e1.p.x
    ny = e1.q.y - e1.p.y
    mx = e2.p.x - e2.q.x
    my = e2.p.y - e2.q.y
    px = e2.p.x - e1.p.x
    py = e2.p.y - e1.p.y

    D = nx*my - ny*mx
    Qx = my*px - mx*py
    Qy = nx*py - ny*px

    (((D > 0) & (Qx < 0)) | ((D > 0) & (Qy < 0)) | ((D < 0) & (Qx > 0)) | ((D < 0) & (Qy > 0))) && return false, e1.p
    ((abs(D) < abs(Qx)) | (abs(D) < abs(Qy))) && return false, e1.p
    abs(D) < tol  && return false, e1.p
    t = Qx/D
    vx = e1.p.x + t*nx
    vy = e1.p.y + t*ny

    return true, V(vx, vy)
end

function islesswithp(s1::E, s2::E, p::V, tol = 1e-15)
    side = sideof(s1, p, tol)
    if side == 1
        return true
    elseif side == -1
        return false
    end
    side = sideof(s1, s2.q, tol)

    # Warning! If two segments coincide, then this function return false!
    if side == 1
        return true
    else
        return false
    end
end

"""
    insert_node!(tree::RBTree, node::RBTreeNode, p::V, tol = 1e-15)

Inserts `node` at proper location by traversing through the `tree` in a binary-search-tree fashion.
"""
function insert_node!(tree::RBTree, node::RBTreeNode, p::V, tol = 1e-15)
    node_y = nothing
    node_x = tree.root

    while node_x !== tree.nil
        node_y = node_x
        if islesswithp(node_x.data, node.data, p, tol)
        # if node.data < node_x.data
            node_x = node_x.rightChild
        else
            node_x = node_x.leftChild
        end
    end

    node.parent = node_y
    if node_y == nothing
        tree.root = node
    elseif islesswithp(node_y.data, node.data, p, tol)
    # elseif node.data < node_y.data
        node_y.rightChild = node
    else
        node_y.leftChild = node
    end
end

"""
    insert!(tree, key, p, tol = 1e-15)

Inserts `key` in the `tree` if it is not present.
"""
function Base.insert!(tree::RBTree{K}, d::K, p::V, tol = 1e-15) where K
    # if the key exists in the tree, no need to insert
    # haskey(tree, d) && return tree

    # insert, if not present in the tree
    node = RBTreeNode{K}(d)
    node.leftChild = node.rightChild = tree.nil

    insert_node!(tree, node, p, tol)

    if node.parent == nothing
        node.color = false
    elseif node.parent.parent == nothing
        ;
    else
        fix_insert!(tree, node)
    end
    tree.count += 1
    return tree
end

"""
    delete!(tree::RBTree, key, p, tol)

Deletes `key` from `tree`, if present, else returns the unmodified tree.
"""
function Base.delete!(tree::RBTree{K}, d::K, p, tol = 1e-15) where K
    z = tree.nil
    node = tree.root

    while node !== tree.nil
        if node.data == d
            z = node
        end

        if islesswithp(d, node.data, p, tol)
        # if d < node.data
            node = node.leftChild
        else
            node = node.rightChild
        end
    end

    (z === tree.nil) && return tree

    y = z
    y_original_color = y.color
    x = RBTreeNode{K}()
    if z.leftChild === tree.nil
        x = z.rightChild
        rb_transplant(tree, z, z.rightChild)
    elseif z.rightChild === tree.nil
        x = z.leftChild
        rb_transplant(tree, z, z.leftChild)
    else
        y = minimum_node(tree, z.rightChild)
        y_original_color = y.color
        x = y.rightChild

        if y.parent == z
            x.parent = y
        else
            rb_transplant(tree, y, y.rightChild)
            y.rightChild = z.rightChild
            y.rightChild.parent = y
        end

        rb_transplant(tree, z, y)
        y.leftChild = z.leftChild
        y.leftChild.parent = y
        y.color = z.color
    end

    !y_original_color && delete_fix(tree, x)
    tree.count -= 1
    return tree
end

function Base.pop!(tree::RBTree{K}, p::V, tol = 1e-15) where K
    z = search_edge(tree, p, tol)

    (z === tree.nil) && return z

    y = z
    y_original_color = y.color
    x = RBTreeNode{K}()
    if z.leftChild === tree.nil
        x = z.rightChild
        rb_transplant(tree, z, z.rightChild)
    elseif z.rightChild === tree.nil
        x = z.leftChild
        rb_transplant(tree, z, z.leftChild)
    else
        y = minimum_node(tree, z.rightChild)
        y_original_color = y.color
        x = y.rightChild

        if y.parent == z
            x.parent = y
        else
            rb_transplant(tree, y, y.rightChild)
            y.rightChild = z.rightChild
            y.rightChild.parent = y
        end

        rb_transplant(tree, z, y)
        y.leftChild = z.leftChild
        y.leftChild.parent = y
        y.color = z.color
    end

    !y_original_color && delete_fix(tree, x)
    tree.count -= 1
    return z
end

function search_edge(tree::RBTree{K}, p::V, tol = 1e-15) where K
    node = tree.root
    while node !== tree.nil
        side = sideof(node.data, p, tol)
        side == 0 && return node
        if side == -1
            node = node.leftChild
        else
            node = node.rightChild
        end
    end
    return node
end

########################################
# AbstractTrees beautifiers
# WARNING!!! TYPE PIRACY!!! DO SOMETHING!!!
########################################

children(tree::RBTree) = children(tree.root)
function children(node::RBTreeNode)
    if node.leftChild === nothing
        if node.rightChild === nothing
            return ()
        else
            return (node.rightChild, )
        end
    else
        if node.rightChild === nothing
            return (node.leftChild, )
        else
            return (node.leftChild, node.rightChild)
        end
    end
end

## Things that make printing prettier
printnode(io::IO, tree::RBTree) = printnode(io, tree.root)
printnode(io::IO, node::RBTreeNode) = print(io, node.data)

function leftmost(tree)
    node = tree.root
    while node.leftChild !== tree.nil
        node = node.leftChild
    end
    return node
end


text:
################################################################################
#
#   METHODS FOR CONSTRUCTION OF BANDSTRUCTURES (ALONG PATHS)
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE BANDSTRUCTURE
#       - type definition
#       - printInfo function
#
#   2) CALCULATION OF BAND STRUCTURES OF UNTICELL OBJECTS
#
#   3) PLOTTING OF BAND STRUCTURES
#       - plotting of Bandstructure objects
#       - plotting of bandstructures of unitcells along paths
#
################################################################################

################################################################################
#
#   TYPE BANDSTRUCTURE
#       - type definition
#       - printInfo function
#
################################################################################
"""
    struct Bandstructure

The type that contains information on a band structure (energy values in momentum space). Fields are

    path   :: Path
    bands  :: Array{Array{Array{Float64, 1}, 1}, 1}

Note that the notation of `bands` is the following:
- `bands[i]`       gives all bands of path segment `i`
- `bands[i][j]`    gives all energy values for band `j` of path segment `i`
- `bands[i][j][k]` gives the energy value at kpoint index `k` of band `j` in path segment `i`

New `Bandstructure` objects can be created only by the default constructor or one of
the several functions to calculate band structures.

# Examples

```julia-repl
julia> bandstructure = Bandstructure(path, bands)
LatticePhysics.Bandstructure(...)
```
"""
struct Bandstructure

    # the path along which the band structure is calcualted
    path::Path

    # bands for each segment
    # bands[i] gives all bands of segment i
    # bands[i][j] gives all energy values for band j of segment i
    # bands[i][j][k] gives the energy value at kpoint index k of band j in segment i
    bands::Array{Array{Array{Float64, 1}, 1}, 1}

    # ONLY DEFAULT CONSTRUCTOR

end

# export the type
export Bandstructure

# INFORMATION FUNCTION
"""
    printInfo(bandstructure::Bandstructure)

Prints information about a `Bandstructure` in terms of
how many bands are in the given object and along which path the structure is calculated.

# Examples

```julia-repl
julia> printInfo(bandstructure)
...

```
"""
function printInfo(bandstructure::Bandstructure)
    # print the header
    println("Bandstructure information:")
    print("$(bandstructure.path.point_names[1])")
    for p in bandstructure.path.point_names[2:end]
        print("\t->-\t$(p)")
    end
    println("")
    # first the point
    print("|")
    # then all segments with closing points
    for s in 1:length(bandstructure.bands)
        # segment s starts
        print("\t")
        # calculate the number of bands
        print("$(length(bandstructure.bands[s])) bands")
        # print the next point
        print("\t|")
    end
    # print a new line
    println("")
end

################################################################################
#
#   BAND STRUCTURE CALCULATION
#
################################################################################
"""
    getBandstructure(
                unitcell::Unitcell,
                path::Path
             [; resolution::Int64=-1,
                enforce_hermitian::Bool=false ]
            )

    getBandstructure(
                matrixFunction::Function,
                path::Path
             [; resolution::Int64=-1 ]
            )

Calculates the band struture of a `Unitcell` object (or a matrix given by a function `matrixFunction`)
along some path given by a `Path` object and returns a `Bandstructure` object.

# Examples

```julia-repl
julia> bandstructure = getBandstructure(unitcell, path)
LatticePhysics.Bandstructure(...)

julia> bandstructure = getBandstructure(unitcell, path, resolution=1000)
LatticePhysics.Bandstructure(...)
```
"""
function getBandstructure(
                unitcell::Unitcell,
                path::Path;
                resolution::Int64=-1,
                enforce_hermitian::Bool=false
            )

    # maybe modify the path resolution
    if resolution > 0
        setTotalResolution!(path, resolution)
    end

    # build up the list of all bands of all segments
    segments_total = Array{Array{Float64,1},1}[]
    # iterate over all path segments and push empty lists into the segments list
    for i in 1:length(path.segment_resolution)
        # build an empty band structure for this segment
        segment = Array{Float64, 1}[]
        for b in 1:length(unitcell.basis)
            push!(segment, zeros(Float64, path.segment_resolution[i]))
        end
        # push the segment band structure into the complete segment list
        push!(segments_total, segment)
    end

    # iterate over all segments
    for s in 1:length(path.segment_resolution)
        # get the grid in between two points
        segment_resolution = path.segment_resolution[s]
        # get all multipliers of k vectors (i.e. all alpha in (1-alpha)*k_1 + alpha*k_2)
        multipliers = linspace(0, 1, segment_resolution)
        # get the local start and end point of the segment
        k1 = path.points[s]
        k2 = path.points[s+1]
        # calculate all energies
        for i in 1:segment_resolution
            # get the current k
            k = (k2 .* multipliers[i]) .+ (k1 .* (1-multipliers[i]))
            # get the interaction matrix for this k
            matrix = getInteractionMatrixKSpace(unitcell, k, enforce_hermitian=enforce_hermitian)
            # diagonalize the matrix
            eigenvalues = eigvals(matrix)
            # save all the eigenvalues to their lists
            for b in 1:length(unitcell.basis)
                segments_total[s][b][i] = eigenvalues[b]
            end
        end
    end

    # generate a new band structure object
    bandstructure = Bandstructure(path, segments_total)

    # return the band structure
    return bandstructure
end
function getBandstructure(
                matrixFunction::Function,
                path::Path;
                resolution::Int64=-1
            )

    # maybe modify the path resolution
    if resolution > 0
        setTotalResolution!(path, resolution)
    end

    # build up the list of all bands of all segments
    segments_total = Array{Array{Float64,1},1}[]
    # get the number of bands
    N_bands = size(matrixFunction(path.points[1]),1)
    # iterate over all path segments and push empty lists into the segments list
    for i in 1:length(path.segment_resolution)
        # build an empty band structure for this segment
        segment = Array{Float64, 1}[]
        for b in 1:N_bands
            push!(segment, zeros(Float64, path.segment_resolution[i]))
        end
        # push the segment band structure into the complete segment list
        push!(segments_total, segment)
    end

    # iterate over all segments
    for s in 1:length(path.segment_resolution)
        # get the grid in between two points
        segment_resolution = path.segment_resolution[s]
        # get all multipliers of k vectors (i.e. all alpha in (1-alpha)*k_1 + alpha*k_2)
        multipliers = linspace(0, 1, segment_resolution)
        # get the local start and end point of the segment
        k1 = path.points[s]
        k2 = path.points[s+1]
        # calculate all energies
        for i in 1:segment_resolution
            # get the current k
            k = (k2 .* multipliers[i]) .+ (k1 .* (1-multipliers[i]))
            # get the interaction matrix for this k
            matrix = matrixFunction(k)
            # diagonalize the matrix
            eigenvalues = eigvals(matrix)
            # save all the eigenvalues to their lists
            for b in 1:length(eigenvalues)
                segments_total[s][b][i] = eigenvalues[b]
            end
        end
    end

    # generate a new band structure object
    bandstructure = Bandstructure(path, segments_total)

    # return the band structure
    return bandstructure
end
export getBandstructure

################################################################################
#
#   BAND STRUCTURE PLOTTING
#
################################################################################

"""
    plotBandstructure(
            bandstructure::Bandstructure
         [; limits_energy="AUTO",
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE" ]
            )

    plotBandstructure(
            unitcell::Unitcell,
            path::Path
         [; resolution::Int64=-1,
            enforce_hermitian::Bool=false,
            ... ]
        )

    plotBandstructure(
            matrixFunction::Function,
            path::Path
         [; resolution::Int64=-1,
            ... ]
        )

Plots the band struture of a passed `Bandstructure` object along some its path and returns the plot as a `PyPlot.Figure` object.
Alternatively, one can pass a `Unitcell` (or a matrix function) and `Path` to calculate the band structure which is plotted.

Additional options include setting plotting related options of `PyPlot` as well as determining if the plot is saved or shown.

# Examples

```julia-repl
julia> plotBandstructure(unitcell, path)
PyPlot.Figure(...)

julia> plotBandstructure(unitcell, path, showPlot=false)
PyPlot.Figure(...)

julia> plotBandstructure(unitcell, save_filename="myplot.pdf")
PyPlot.Figure(...)

julia> plotBandstructure(bandstructure)
PyPlot.Figure(...)
```
"""
function plotBandstructure(
            bandstructure::Bandstructure;
            limits_energy="AUTO",
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # get the path from the bandstructure
    path = bandstructure.path

    # configure plot environment
    rc("font", family="serif")

    # create a new figure
    fig = figure(figsize=figsize)

    ###########################
    #   PLOT BANDS
    ###########################

    # plot the band structure
    for s in 1:length(bandstructure.bands)
        # plot the segment
        for b in bandstructure.bands[s]
            # plot the band
            plot(
                collect(1:path.segment_resolution[s]) .+ sum(path.segment_resolution[1:s-1]), # xvalues
                b,
                "-$(plot_color)"
            )
        end
    end

    ###########################
    #   SET ALL TICKS (POINTS)
    ###########################

    # get the current axis
    ax = gca()
    axx = ax[:get_xaxis]()
    # compile tick positions and labels
    point_pos = Int64[]
    push!(point_pos, 1)
    for l in 1:length(path.segment_resolution)
        push!(point_pos, sum(path.segment_resolution[1:l]))
    end
    point_labels = String[path.point_names[i] for i in 1:length(path.points)]
    # configure tick labels
    xticks(point_pos, point_labels)
    # configure ticks
    axx[:set_tick_params](which="both", direction="out")
    axx[:set_tick_params](which="top", color="none")
    axy = ax[:get_yaxis]()
    axy[:set_tick_params](which="both", direction="out")

    # plot vertical lines for each point
    for p in point_pos
        axvline(p,color=[0.6, 0.6, 0.6], linestyle="--")
    end

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # label the axis
    xlabel("momentum")
    ylabel("energy")

    # energy limits
    # check if specific boundaries are desired
    if !(limits_energy == "AUTO")
        ylim(limits_energy[1], limits_energy[2])
    end

    # momentum limits (x axis)
    xlim(0, maximum(point_pos)+1)

    # set the title
    if plot_title == "AUTO"
        # set the title to an automatically generated title
        title("energy spectrum along path $(getPathString(path))")
    elseif plot_title == ""
        # do nothing title related
    else
        # set the title to the given title
        title(plot_title)
    end

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # save the plot
    if save_filename != "NONE"
        # make sure the directory exists
        if contains(save_filename, "/")
    		# get the containing folder
    		folder = save_filename[1:findlast(save_filename, '/')]
    		# build the path to that folder
    		mkpath(folder)
    	end
        # save the plot
        savefig(save_filename)
    end

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return fig
end
function plotBandstructure(
            unitcell::Unitcell,
            path::Path;
            resolution::Int64=-1,
            enforce_hermitian::Bool=false,
            limits_energy="AUTO",
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )
    # calculate the bandstructure
    bandstructure = getBandstructure(unitcell, path, resolution=resolution, enforce_hermitian=enforce_hermitian)
    # call the respective function
    return plotBandstructure(
                bandstructure;
                limits_energy=limits_energy,
                plot_title=plot_title,
                plot_color=plot_color,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
end
function plotBandstructure(
            matrixFunction::Function,
            path::Path;
            resolution::Int64=-1,
            limits_energy="AUTO",
            plot_title::String="",
            plot_color="b",
            figsize::Tuple=(6,4),
            showPlot::Bool=true,
            save_filename::String="NONE"
        )
    # calculate the bandstructure
    bandstructure = getBandstructure(matrixFunction, path, resolution=resolution)
    # call the respective function
    return plotBandstructure(
                bandstructure;
                limits_energy=limits_energy,
                plot_title=plot_title,
                plot_color=plot_color,
                figsize=figsize,
                showPlot=showPlot,
                save_filename=save_filename
            )
end
export plotBandstructure


text:
#=
This file contains routines that are common to both BD and IBD
=#

function func_grad!(storage, x::AbstractVector, pa, attrib)
	global to

	# x to pa.optm.cal.s or pa.optm.cal.g
	x_to_model!(x, pa, attrib)

	@timeit to "F" F!(pa, x, attrib) # forward

	if(storage === nothing)
		# compute misfit and δdcal
		f = Misfits.error_squared_euclidean!(nothing, pa.optm.cal.d, pa.optm.obs.d,
				       nothing, norm_flag=false)
	else
		f = Misfits.error_squared_euclidean!(pa.optm.ddcal, pa.optm.cal.d,
				       pa.optm.obs.d, nothing, norm_flag=false)
		@timeit to "Fadj" Fadj!(pa, storage, x, pa.optm.ddcal, attrib)
	end
	return f

end

# exponential-weighted norm for the green functions
function func_grad_g_weights!(storage, x, pa)
	x_to_model!(x, pa)
	!(pa.attrib_inv == :g) && error("only for g inversion")
	if(!(storage === nothing)) #
		f = Misfits.error_weighted_norm!(pa.optm.dg,pa.optm.cal.g, pa.gx.weights) #
		for i in eachindex(storage)
			storage[i]=pa.optm.dg[i]
		end
	else
		f = Misfits.error_weighted_norm!(nothing,pa.optm.cal.g, pa.gx.weights)
	end
	return f
end

# exponential-weighted norm for the green functions
function func_grad_g_acorr_weights!(storage, x, pa)
	x_to_model!(x, pa)
	!(pa.attrib_inv == :g) && error("only for g inversion")

	if(!(storage === nothing)) #
		f = Misfits.error_acorr_weighted_norm!(pa.optm.dg,pa.optm.cal.g,
					 paconv=pa.g_acorr,dfds=pa.optm.dg_acorr) #
		for i in eachindex(storage)
			storage[i]=pa.optm.dg[i]
		end
	else
		f = Misfits.error_acorr_weighted_norm!(nothing,pa.optm.cal.g,
					 paconv=pa.g_acorr,dfds=pa.optm.dg_acorr)
	end

	return f
end
#