text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using REPL.REPLCompletions
using Test
using Random
using REPL
    @testset "Check symbols previously not shown by REPL.doc_completions()" begin
    symbols = ["?","=","[]","[","]","{}","{","}",";","","'","&&","||","julia","Julia","new","@var_str"]
        for i in symbols
            @test REPL.doc_completions(i)[1]==i
        end
    end
let ex = quote
    module CompletionFoo
        using Random
        import Test

        mutable struct Test_y
            yy
        end
        mutable struct Test_x
            xx :: Test_y
        end
        type_test = Test_x(Test_y(1))
        (::Test_y)() = "", ""
        module CompletionFoo2

        end
        const bar = 1
        foo() = bar
        macro foobar()
            :()
        end
        macro barfoo(ex)
            ex
        end
        macro error_expanding()
            error("cannot expand @error_expanding")
            :()
        end
        macro error_lowering_conditional(a)
            if isa(a, Number)
                return a
            end
            throw(AssertionError("Not a Number"))
            :()
        end
        macro error_throwing()
            return quote
                error("@error_throwing throws an error")
            end
        end

        primitive type NonStruct 8 end
        Base.propertynames(::NonStruct) = (:a, :b, :c)
        x = reinterpret(NonStruct, 0x00)

        # Support non-Dict AbstractDicts, #19441
        mutable struct CustomDict{K, V} <: AbstractDict{K, V}
            mydict::Dict{K, V}
        end

        Base.keys(d::CustomDict) = collect(keys(d.mydict))
        Base.length(d::CustomDict) = length(d.mydict)

        test(x::T, y::T) where {T<:Real} = pass
        test(x::Real, y::Real) = pass
        test(x::AbstractArray{T}, y) where {T<:Real} = pass
        test(args...) = pass

        test1(x::Type{Float64}) = pass

        test2(x::AbstractString) = pass
        test2(x::Char) = pass
        test2(x::Cmd) = pass

        test3(x::AbstractArray{Int}, y::Int) = pass
        test3(x::AbstractArray{Float64}, y::Float64) = pass

        test4(x::AbstractString, y::AbstractString) = pass
        test4(x::AbstractString, y::Regex) = pass

        test5(x::Array{Bool,1}) = pass
        test5(x::BitArray{1}) = pass
        test5(x::Float64) = pass
        const a=x->x
        test6()=[a, a]
        test7() = rand(Bool) ? 1 : 1.0
        test8() = Any[1][1]
        test9(x::Char) = pass
        test9(x::Char, i::Int) = pass

        test10(a, x::Int...) = pass
        test10(a::Integer, b::Integer, c) = pass
        test10(a, y::Bool...) = pass
        test10(a, d::Integer, z::Signed...) = pass
        test10(s::String...) = pass

        test11(a::Integer, b, c) = pass
        test11(u, v::Integer, w) = pass
        test11(x::Int, y::Int, z) = pass
        test11(_, _, s::String) = pass

        kwtest(; x=1, y=2, w...) = pass
        kwtest2(a; x=1, y=2, w...) = pass

        array = [1, 1]
        varfloat = 0.1

        const tuple = (1, 2)

        test_y_array=[CompletionFoo.Test_y(rand()) for i in 1:10]
        test_dict = Dict("abc"=>1, "abcd"=>10, :bar=>2, :bar2=>9, Base=>3,
                         occursin=>4, `ls`=>5, 66=>7, 67=>8, ("q",3)=>11,
                         "α"=>12, :α=>13)
        test_customdict = CustomDict(test_dict)

        macro teststr_str(s) end
        macro tϵsτstρ_str(s) end
        macro testcmd_cmd(s) end
        macro tϵsτcmδ_cmd(s) end

        end
        test_repl_comp_dict = CompletionFoo.test_dict
        test_repl_comp_customdict = CompletionFoo.test_customdict
        test_dict_ℂ = Dict(1=>2)
    end
    ex.head = :toplevel
    Core.eval(Main, ex)
end

function map_completion_text(completions)
    c, r, res = completions
    return map(completion_text, c), r, res
end

test_complete(s) = map_completion_text(@inferred(completions(s, lastindex(s))))
test_scomplete(s) =  map_completion_text(@inferred(shell_completions(s, lastindex(s))))
test_bslashcomplete(s) =  map_completion_text(@inferred(bslash_completions(s, lastindex(s)))[2])
test_complete_context(s, m) =  map_completion_text(@inferred(completions(s,lastindex(s), m)))
test_complete_foo(s) = test_complete_context(s, Main.CompletionFoo)
test_complete_noshift(s) = map_completion_text(@inferred(completions(s, lastindex(s), Main, false)))

module M32377 end
test_complete_32377(s) = map_completion_text(completions(s,lastindex(s), M32377))

let s = ""
    c, r = test_complete(s)
    @test "CompletionFoo" in c
    @test isempty(r)
    @test s[r] == ""
end

let s = "using REP"
    c, r = test_complete_32377(s)
    @test count(isequal("REPL"), c) == 1
    # issue #30234
    @test !Base.isbindingresolved(M32377, :tanh)
    # check what happens if REPL is already imported
    M32377.eval(:(using REPL))
    c, r = test_complete_32377(s)
    @test count(isequal("REPL"), c) == 1
end

let s = "Comp"
    c, r = test_complete(s)
    @test "CompletionFoo" in c
    @test r == 1:4
    @test s[r] == "Comp"
end

let s = "Main.Comp"
    c, r = test_complete(s)
    @test "CompletionFoo" in c
    @test r == 6:9
    @test s[r] == "Comp"
end

let s = "Main.CompletionFoo."
    c, r = test_complete(s)
    @test "bar" in c
    @test r === 20:19
    @test s[r] == ""
end

let s = "Main.CompletionFoo.f"
    c, r = test_complete(s)
    @test "foo" in c
    @test r == 20:20
    @test s[r] == "f"
    @test !("foobar" in c)
end

# test method completions when `!` operator precedes
let
    s = "!is"
    c, r = test_complete(s)
    @test "isa" in c
    @test s[r] == "is"
    @test !("!" in c)

    s = "!!is"
    c, r = test_complete(s)
    @test "isa" in c
    @test s[r] == "is"
    @test !("!" in c)
end

# issue #6424
let s = "Main.CompletionFoo.@f"
    c, r = test_complete(s)
    @test "@foobar" in c
    @test r == 20:21
    @test s[r] == "@f"
    @test !("foo" in c)
end

let s = "Main.CompletionFoo.type_test.x"
    c, r = test_complete(s)
    @test "xx" in c
    @test r == 30:30
    @test s[r] == "x"
end

let s = "Main.CompletionFoo.bar.no_val_available"
    c, r = test_complete(s)
    @test length(c)==0
end

#cannot do dot completion on infix operator
let s = "+."
    c, r = test_complete(s)
    @test length(c)==0
end

# To complete on a variable of a type, the type T of the variable
# must be a concrete type, hence Base.isstructtype(T) returns true,
# for the completion to succeed. That why `xx :: Test_y` of `Test_x`.
let s = "Main.CompletionFoo.type_test.xx.y"
    c, r = test_complete(s)
    @test "yy" in c
    @test r == 33:33
    @test s[r] == "y"
end

# issue #6333
let s = "Base.return_types(getin"
    c, r = test_complete(s)
    @test "getindex" in c
    @test r == 19:23
    @test s[r] == "getin"
end

# issue #23193: after `using`, identifiers can be prefixed by module names
let s = "using Test, Random"
    c, r = test_complete(s)
    @test !("RandomDevice" in c)
end

# issue #23226: identifiers must be separated by a comma (not a newline)
let s = "using Base\nusi"
    c, r = test_complete(s)
    @test "using" in c
end

# issue 23292
let
    @test_nowarn test_complete("test7().")
    c, r = test_complete("test8().")
    @test isempty(c)
end

# inexistent completion inside a string
let s = "Base.print(\"lol"
    c, r, res = test_complete(s)
    @test res == false
end

# inexistent completion inside a cmd
let s = "run(`lol"
    c, r, res = test_complete(s)
    @test res == false
end

# test latex symbol completions
let s = "\\alpha"
    c, r = test_bslashcomplete(s)
    @test c[1] == "α"
    @test r == 1:length(s)
    @test length(c) == 1
end

# test latex symbol completions after unicode #9209
let s = "α\\alpha"
    c, r = test_bslashcomplete(s)
    @test c[1] == "α"
    @test r == 3:sizeof(s)
    @test length(c) == 1
end

# test emoji symbol completions
let s = "\\:koala:"
    c, r = test_bslashcomplete(s)
    @test c[1] == "🐨"
    @test r == 1:sizeof(s)
    @test length(c) == 1
end

let s = "\\:ko"
    c, r = test_bslashcomplete(s)
    @test "\\:koala:" in c
end

# test emoji symbol completions after unicode #9209
let s = "α\\:koala:"
    c, r = test_bslashcomplete(s)
    @test c[1] == "🐨"
    @test r == 3:sizeof(s)
    @test length(c) == 1
end

# test latex symbol completions in strings should not work when there
# is a backslash in front of `\alpha` because it interferes with path completion on windows
let s = "cd(\"path_to_an_empty_folder_should_not_complete_latex\\\\\\alpha"
    c, r, res = test_complete(s)
    @test length(c) == 0
end

# test latex symbol completions in strings
let s = "\"C:\\\\ \\alpha"
    c, r, res = test_complete(s)
    @test c[1] == "α"
    @test r == 7:12
    @test length(c) == 1
end

# test latex symbol completion in getindex expressions (#24705)
let s = "tuple[\\alpha"
    c, r, res = test_complete_foo(s)
    @test c[1] == "α"
    @test r == 7:12
    @test length(c) == 1
end

let s = "\\a"
    c, r, res = test_complete(s)
    "\\alpha" in c
    @test r == 1:2
    @test s[r] == "\\a"
end

# `cd("C:\U should not make the repl crash due to escaping see comment #9137
let s = "cd(\"C:\\U"
    c, r, res = test_complete(s)
end

# Test method completions
let s = "max("
    c, r, res = test_complete(s)
    @test !res
    @test let found = false
        for m in methods(max)
            if !found
                found = (c[1] == string(m))
            end
        end
        found
    end
    @test r == 1:3
    @test s[r] == "max"
end

# test method completions when `!` operator precedes
let
    s = "!("
    c, r, res = test_complete(s)
    @test !res
    @test all(m -> string(m) in c, methods(!))
    @test s[r] == s[1:end-1]

    s = "!isnothing("
    c, r, res = test_complete(s)
    @test !res
    @test all(m -> string(m) in c, methods(isnothing))
    @test s[r] == s[1:end-1]

    s = "!!isnothing("
    c, r, res = test_complete(s)
    @test !res
    @test all(m -> string(m) in c, methods(isnothing))
    @test s[r] == s[1:end-1]
end

# Test completion of methods with input concrete args and args where typeinference determine their type
let s = "CompletionFoo.test(1, 1, "
    c, r, res = test_complete(s)
    @test !res
    @test c[1] == string(first(methods(Main.CompletionFoo.test, Tuple{Int, Int})))
    @test c[2] == string(first(methods(Main.CompletionFoo.test, Tuple{}))) # corresponding to the vararg
    @test length(c) == 2
    # In particular, this checks that test(x::Real, y::Real) is not a valid completion
    # since it is strictly less specific than test(x::T, y::T) where T
    @test r == 1:18
    @test s[r] == "CompletionFoo.test"
end

let s = "CompletionFoo.test(CompletionFoo.array,"
    c, r, res = test_complete(s)
    @test !res
    @test c[1] == string(first(methods(Main.CompletionFoo.test, Tuple{Array{Int, 1}, Any})))
    @test length(c) == 2
    @test r == 1:18
    @test s[r] == "CompletionFoo.test"
end

let s = "CompletionFoo.test(1,1,1,"
    c, r, res = test_complete(s)
    @test !res
    @test c[1] == string(first(methods(Main.CompletionFoo.test, Tuple{Any, Any, Any})))
    @test length(c) == 1
    @test r == 1:18
    @test s[r] == "CompletionFoo.test"
end

let s = "CompletionFoo.test1(Int,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 0
    @test r == 1:19
    @test s[r] == "CompletionFoo.test1"
end

let s = "CompletionFoo.test1(Float64,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test r == 1:19
    @test s[r] == "CompletionFoo.test1"
end

let s = "prevind(\"θ\",1,"
    c, r, res = test_complete(s)
    @test c[1] == string(first(methods(prevind, Tuple{String, Int})))
    @test r == 1:7
    @test s[r] == "prevind"
end

for (T, arg) in [(String,"\")\""),(Char, "')'")]
    s = "(1, CompletionFoo.test2($arg,"
    c, r, res = test_complete(s)
    @test length(c) == 1
    @test c[1] == string(first(methods(Main.CompletionFoo.test2, Tuple{T,})))
    @test r == 5:23
    @test s[r] == "CompletionFoo.test2"
end

let s = "(1, CompletionFoo.test2(`')'`,"
    c, r, res = test_complete(s)
    @test length(c) == 1
    @test c[1] == string(first(methods(Main.CompletionFoo.test2, Tuple{Cmd})))
end

let s = "CompletionFoo.test3([1, 2] .+ CompletionFoo.varfloat,"
    c, r, res = test_complete(s)
    @test !res
    @test_broken only(c) == string(first(methods(Main.CompletionFoo.test3, Tuple{Array{Float64, 1}, Float64})))
end

let s = "CompletionFoo.test3([1.,2.], 1.,"
    c, r, res = test_complete(s)
    @test !res
    @test c[1] == string(first(methods(Main.CompletionFoo.test3, Tuple{Array{Float64, 1}, Float64})))
    @test r == 1:19
    @test length(c) == 1
    @test s[r] == "CompletionFoo.test3"
end

let s = "CompletionFoo.test4(\"e\",r\" \","
    c, r, res = test_complete(s)
    @test !res
    @test c[1] == string(first(methods(Main.CompletionFoo.test4, Tuple{String, Regex})))
    @test r == 1:19
    @test length(c) == 1
    @test s[r] == "CompletionFoo.test4"
end

# (As discussed in #19829, the Base.REPLCompletions.get_type function isn't
#  powerful enough to analyze anonymous functions.)
let s = "CompletionFoo.test5(broadcast((x,y)->x==y, push!(Base.split(\"\",' '),\"\",\"\"), \"\"),"
    c, r, res = test_complete(s)
    @test !res
    @test_broken only(c) == string(first(methods(Main.CompletionFoo.test5, Tuple{BitArray{1}})))
end

# test partial expression expansion
let s = "CompletionFoo.test5(Bool[x==1 for x=1:4],"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test c[1] == string(first(methods(Main.CompletionFoo.test5, Tuple{Array{Bool,1}})))
end

let s = "CompletionFoo.test4(CompletionFoo.test_y_array[1]()[1], CompletionFoo.test_y_array[1]()[2], "
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test c[1] == string(first(methods(Main.CompletionFoo.test4, Tuple{String, String})))
end

# Test that string escaption is handled correct
let s = """CompletionFoo.test4("\\"","""
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 2
end

########## Test where the current inference logic fails ########
# Fails due to inference fails to determine a concrete type for arg 1
# But it returns AbstractArray{T,N} and hence is able to remove test5(x::Float64) from the suggestions
let s = "CompletionFoo.test5(AbstractArray[[]][1],"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 2
end

# equivalent to above but due to the time macro the completion fails to find the concrete type
let s = "CompletionFoo.test3(@time([1, 2] + CompletionFoo.varfloat),"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 2
end
#################################################################

# method completions with kwargs
let s = "CompletionFoo.kwtest( "
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test occursin("x, y, w...", c[1])
end

for s in ("CompletionFoo.kwtest(;",
          "CompletionFoo.kwtest(; x=1, ",
          "CompletionFoo.kwtest(; kw=1, ",
          )
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test occursin("x, y, w...", c[1])
end

for s in ("CompletionFoo.kwtest2(1; x=1,",
          "CompletionFoo.kwtest2(1; kw=1, ",
          )
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test occursin("a; x, y, w...", c[1])
end

#################################################################

# method completion with `?` (arbitrary method with given argument types)
let s = "CompletionFoo.?([1,2,3], 2.0)"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 1
    @test occursin("test(x::AbstractArray{T}, y) where T<:Real", c[1])
    # In particular, this checks that test(args...) is not a valid completion
    # since it is strictly less specific than test(x::AbstractArray{T}, y)
end

let s = "CompletionFoo.?([1,2,3], 2.0"
    c, r, res = test_complete(s)
    @test !res
    @test  any(str->occursin("test(x::AbstractArray{T}, y) where T<:Real", str), c)
    @test  any(str->occursin("test(args...)", str), c)
    @test !any(str->occursin("test3(x::AbstractArray{Int", str), c)
    @test !any(str->occursin("test4", str), c)
end

let s = "CompletionFoo.?('c')"
    c, r, res = test_complete(s)
    @test !res
    @test  any(str->occursin("test9(x::Char)", str), c)
    @test  any(str->occursin("test10(a, ", str), c)
    @test !any(str->occursin("test9(x::Char, i::Int", str), c)
end

let s = "CompletionFoo.?('c'"
    c, r, res = test_complete(s)
    @test !res
    @test  any(str->occursin("test9(x::Char)", str), c)
    @test  any(str->occursin("test10(a, ", str), c)
    @test  any(str->occursin("test9(x::Char, i::Int", str), c)
end

let s = "CompletionFoo.?(false, \"a\", 3, "
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 2
    @test occursin("test(args...)", c[1])
    @test occursin("test11(a::Integer, b, c)", c[2])
end

let s = "CompletionFoo.?(false, \"a\", 3, "
    c, r, res = test_complete_noshift(s)
    @test !res
    @test length(c) == 1
    @test occursin("test11(a::Integer, b, c)", c[1])
end

let s = "CompletionFoo.?(\"a\", 3, "
    c, r, res = test_complete(s)
    @test !res
    @test  any(str->occursin("test10(a, x::$Int...)", str), c)
    @test !any(str->occursin("test10(a, y::Bool...)", str), c)
    @test !any(str->occursin("test10(s::String...)", str), c)
end

let s = "CompletionFoo.?()"
    c, r, res = test_complete(s)
    @test !res
    @test any(str->occursin("foo()", str), c)
    @test any(str->occursin("kwtest(;", str), c)
    @test any(str->occursin("test(args...)", str), c)
end

let s = "CompletionFoo.?()"
    c, r, res = test_complete_noshift(s)
    @test !res
    @test length(c) == 1
    @test occursin("test10(s::String...)", c[1])
end

#################################################################

# Test method completion with varargs
let s = "CompletionFoo.test10(z, Integer[]...,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 5
    @test all(startswith("test10("), c)
    @test allunique(c)
end

let s = "CompletionFoo.test10(3, Integer[]...,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 4
    @test all(startswith("test10("), c)
    @test allunique(c)
    @test !any(str->occursin("test10(s::String...)", str), c)
end

let s = "CompletionFoo.test10(3, 4,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 3
    @test any(str->occursin("test10(a, x::$Int...)", str), c)
    @test any(str->occursin("test10(a::Integer, b::Integer, c)", str), c)
    @test any(str->occursin("test10(a, d::Integer, z::Signed...)", str), c)
end

let s = "CompletionFoo.test10(3, 4, 5,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 3
    @test any(str->occursin("test10(a, x::$Int...)", str), c)
    @test any(str->occursin("test10(a::Integer, b::Integer, c)", str), c) # show it even though the call would result in an ambiguity error
    @test any(str->occursin("test10(a, d::Integer, z::Signed...)", str), c)
    # the last one is not eliminated by specificity since the complete call could be
    # test10(3, 4, 5, Int8(6)) for instance
end

let s = "CompletionFoo.test10(z, z, 0, "
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 3
    @test any(str->occursin("test10(a, x::$Int...)", str), c)
    @test any(str->occursin("test10(a::Integer, b::Integer, c)", str), c) # show it even though the call would result in an ambiguity error
    @test any(str->occursin("test10(a, d::Integer, z::Signed...)", str), c)
end

let s = "CompletionFoo.test10(\"a\", Union{Signed,Bool,String}[3][1], "
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 4
    @test all(startswith("test10("), c)
    @test allunique(c)
    @test !any(str->occursin("test10(a::Integer, b::Integer, c)", str), c)
end

# Test method completion with ambiguity
let s = "CompletionFoo.test11(Integer[false][1], Integer[14][1], "
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 4
    @test all(startswith("test11("), c)
    @test allunique(c)
end

let s = "CompletionFoo.test11(Integer[-7][1], Integer[0x6][1], 6,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 3
    @test any(str->occursin("test11(a::Integer, b, c)", str), c)
    @test any(str->occursin("test11(u, v::Integer, w)", str), c)
    @test any(str->occursin("test11(x::$Int, y::$Int, z)", str), c)
end

let s = "CompletionFoo.test11(3, 4,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 4
    @test any(str->occursin("test11(x::$Int, y::$Int, z)", str), c)
    @test any(str->occursin("test11(::Any, ::Any, s::String)", str), c)
end

let s = "CompletionFoo.test11(0x8, 5,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 3
    @test any(str->occursin("test11(a::Integer, b, c)", str), c)
    @test any(str->occursin("test11(u, v::Integer, w)", str), c)
    @test any(str->occursin("test11(::Any, ::Any, s::String)", str), c)
end

let s = "CompletionFoo.test11(0x8, 'c',"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 2
    @test any(str->occursin("test11(a::Integer, b, c)", str), c)
    @test any(str->occursin("test11(::Any, ::Any, s::String)", str), c)
end

let s = "CompletionFoo.test11('d', 3,"
    c, r, res = test_complete(s)
    @test !res
    @test length(c) == 2
    @test any(str->occursin("test11(u, v::Integer, w)", str), c)
    @test any(str->occursin("test11(::Any, ::Any, s::String)", str), c)
end

# Test of inference based getfield completion
let s = "(1+2im)."
    c,r = test_complete(s)
    @test length(c) == 2
    @test r == (lastindex(s) + 1):lastindex(s)
    @test c == ["im", "re"]
end

let s = "((1+2im))."
    c, r = test_complete(s)
    @test length(c) == 2
    @test r == (lastindex(s) + 1):lastindex(s)
    @test c == ["im", "re"]
end

let s = "CompletionFoo.test_y_array[1]."
    c, r = test_complete(s)
    @test length(c) == 1
    @test r == (lastindex(s) + 1):lastindex(s)
    @test c[1] == "yy"
end

let s = "CompletionFoo.Test_y(rand()).y"
    c, r = test_complete(s)
    @test length(c) == 1
    @test r == lastindex(s):lastindex(s)
    @test c[1] == "yy"
end

let s = "CompletionFoo.test6()[1](CompletionFoo.Test_y(rand())).y"
    c, r = test_complete(s)
    @test length(c) == 1
    @test r == lastindex(s):lastindex(s)
    @test c[1] == "yy"
end

# Test completion in multi-line comments
let s = "#=\n\\alpha"
    c, r, res = test_complete(s)
    @test c[1] == "α"
    @test r == 4:9
    @test length(c) == 1
end

# Test that completion do not work in multi-line comments
let s = "#=\nmax"
    c, r, res = test_complete(s)
    @test length(c) == 0
end

# Test completion of packages
path = joinpath(tempdir(),randstring())
pushfirst!(LOAD_PATH, path)
try
    mkpath(path)
    write(joinpath(path, "Project.toml"),
        """
        name = "MyProj"

        [deps]
        MyPack = "09ebe64f-f76c-4f21-bef2-bd9be6c77e76"
        """)
        c, r, res = test_complete("using MyP")
        @test "MyPack" in c
        @test "MyProj" in c
finally
    @test popfirst!(LOAD_PATH) == path
    rm(path, recursive=true)
end

path = joinpath(tempdir(),randstring())
push!(LOAD_PATH, path)
try
    # Should not throw an error even though the path do no exist
    test_complete("using ")
    Pack_folder = joinpath(path, "Test_pack")
    mkpath(Pack_folder)

    Pack_folder2 = joinpath(path, "Test_pack2", "src")
    mkpath(Pack_folder2)
    touch(joinpath(Pack_folder2, "Test_pack2.jl"))

    # Test it completes on folders
    local c, r, res # workaround for issue #24331
    c, r, res = test_complete("using Test_p")
    @test !("Test_pack" in c)
    @test "Test_pack2" in c
finally
    @test pop!(LOAD_PATH) == path
    rm(path, recursive=true)
end

# Test $ in shell-mode
let s = "cd \$(max"
    c, r, res = test_scomplete(s)
    @test "max" in c
    @test r == 6:8
    @test s[r] == "max"
end

# The return type is of importance, before #8995 it would return nothing
# which would raise an error in the repl code.
let c, r, res
    c, r, res = test_scomplete("\$a")
    @test c == String[]
    @test r === 0:-1
    @test res === false
end

if Sys.isunix()
let s, c, r
    #Assume that we can rely on the existence and accessibility of /tmp

    # Tests path in Julia code and closing " if it's a file
    # Issue #8047
    s = "@show \"/dev/nul"
    c,r = test_complete(s)
    @test "null\"" in c
    @test r == 13:15
    @test s[r] == "nul"

    # Tests path in Julia code and not closing " if it's a directory
    # Issue #8047
    s = "@show \"/tm"
    c,r = test_complete(s)
    @test "tmp/" in c
    @test r == 9:10
    @test s[r] == "tm"

    # Tests path in Julia code and not double-closing "
    # Issue #8047
    s = "@show \"/dev/nul\""
    c,r = completions(s, 15)
    c = map(completion_text, c)
    @test "null" in c
    @test r == 13:15
    @test s[r] == "nul"

    s = "/t"
    c,r = test_scomplete(s)
    @test "tmp/" in c
    @test r == 2:2
    @test s[r] == "t"

    s = "/tmp"
    c,r = test_scomplete(s)
    @test "tmp/" in c
    @test r == 2:4
    @test s[r] == "tmp"

    # This should match things that are inside the tmp directory
    s = tempdir()
    if !endswith(s, "/")
        s = string(s, "/")
    end
    if !isdir(joinpath(s, "tmp"))
        c,r = test_scomplete(s)
        @test !("tmp/" in c)
        @test r === length(s) + 1:0
        @test s[r] == ""
    end

    s = "cd \$(Iter"
    c,r = test_scomplete(s)
    @test "Iterators" in c
    @test r == 6:9
    @test s[r] == "Iter"

    # Pressing tab after having entered "/tmp " should not
    # attempt to complete "/tmp" but rather work on the current
    # working directory again.
    let s, c, r, file
        file = joinpath(path, "repl completions")
        s = "/tmp "
        c,r = test_scomplete(s)
        @test r === 6:5
    end

    # Test completing paths with an escaped trailing space
    let s, c, r, file
        file = joinpath(tempdir(), "repl completions")
        touch(file)
        s = string(tempdir(), "/repl\\ ")
        c,r = test_scomplete(s)
        @test ["repl\\ completions"] == c
        @test s[r] == "repl\\ "
        rm(file)
    end

    # Tests homedir expansion
    mktempdir() do tmphome
        withenv("HOME" => tmphome, "USERPROFILE" => tmphome) do
            path = homedir()
            dir = joinpath(path, "tmpfoobar")
            mkdir(dir)
            s = "\"" * path * "/tmpfoob"
            c,r = test_complete(s)
            @test "tmpfoobar/" in c
            l = 3 + length(path)
            @test r == l:l+6
            @test s[r] == "tmpfoob"
            s = "\"~"
            @test "tmpfoobar/" in c
            c,r = test_complete(s)
            rm(dir)
        end
    end

    # Tests detecting of files in the env path (in shell mode)
    let path, file
        path = tempdir()
        unreadable = joinpath(tempdir(), "replcompletion-unreadable")

        try
            file = joinpath(path, "tmp-executable")
            touch(file)
            chmod(file, 0o755)
            mkdir(unreadable)
            chmod(unreadable, 0o000)

            # PATH can also contain folders which we aren't actually allowed to read.
            withenv("PATH" => string(path, ":", unreadable)) do
                s = "tmp-execu"
                c,r = test_scomplete(s)
                @test "tmp-executable" in c
                @test r == 1:9
                @test s[r] == "tmp-execu"
            end
        finally
            rm(file)
            rm(unreadable)
        end
    end

    # Make sure completion results are unique in case things are in the env path twice.
    let
        file0 = joinpath(tempdir(), "repl-completion")
        dir = joinpath(tempdir(), "repl-completion-subdir")
        file1 = joinpath(dir, "repl-completion")

        try
            # Create /tmp/repl-completion and /tmp/repl-completion-subdir/repl-completion
            mkdir(dir)
            touch(file0)
            touch(file1)

            withenv("PATH" => string(tempdir(), ":", dir)) do
                s = string("repl-completio")
                c,r = test_scomplete(s)
                @test ["repl-completion"] == c
                @test s[r] == "repl-completio"
            end

        finally
            rm(file0)
            rm(file1)
            rm(dir)
        end
    end
end
end

#test that it does not crash on files for which `stat` errors
let current_dir, forbidden
    # Issue #36855
    if !Sys.iswindows() || Sys.windows_version() >= Sys.WINDOWS_VISTA_VER
        mktempdir() do path
            selfsymlink = joinpath(path, "selfsymlink")
            symlink(selfsymlink, selfsymlink)
            @test try
                stat(selfsymlink) # should crash with a IOError
                false
            catch e
                e isa Base.IOError && occursin("ELOOP", e.msg)
            end
            c, r = test_complete("\"$(joinpath(path, "selfsym"))")
            @test c == ["selfsymlink"]
        end
    end

    # Issue #32797
    forbidden = Sys.iswindows() ? "C:\\S" : "/root/x"
    test_complete(forbidden); @test true # simply check that it did not crash

     # Issue #19310
    if Sys.iswindows()
        current_dir = pwd()
        cd("C:\\")
        test_complete("C"); @test true
        test_complete("C:"); @test true
        test_complete("C:\\"); @test true
        if isdir("D:\\")
            cd("D:\\")
            test_complete("C"); @test true
            test_complete("C:"); @test true
            test_complete("C:\\"); @test true
        end
        cd(current_dir)
    end
end

#test that it can auto complete with spaces in file/path
mktempdir() do path
    space_folder = randstring() * " α"
    dir = joinpath(path, space_folder)
    dir_space = replace(space_folder, " " => "\\ ")

    mkdir(dir)
    cd(path) do
        open(joinpath(space_folder, "space .file"),"w") do f
            s = Sys.iswindows() ? "rm $dir_space\\\\space" : "cd $dir_space/space"
            c, r = test_scomplete(s)
            @test r == lastindex(s)-4:lastindex(s)
            @test "space\\ .file" in c

            s = Sys.iswindows() ? "cd(\"β $dir_space\\\\space" : "cd(\"β $dir_space/space"
            c, r = test_complete(s)
            @test r == lastindex(s)-4:lastindex(s)
            @test "space .file\"" in c
        end
        # Test for issue #10324
        s = "cd(\"$dir_space"
        c, r = test_complete(s)
        @test r == 5:15
        @test s[r] ==  dir_space

        #Test for #18479
        for c in "'`@\$;&"
            test_dir = "test$(c)test"
            mkdir(joinpath(path, test_dir))
            try
                if !(c in ['\'','$']) # As these characters hold special meaning
                    # in shell commands the shell path completion cannot complete
                    # paths with these characters
                    c, r, res = test_scomplete(test_dir)
                    @test c[1] == test_dir*(Sys.iswindows() ? "\\\\" : "/")
                    @test res
                end
                c, r, res = test_complete("\""*test_dir)
                @test c[1] == test_dir*(Sys.iswindows() ? "\\\\" : "/")
                @test res
            finally
                rm(joinpath(path, test_dir), recursive=true)
            end
        end
    end
    rm(dir, recursive=true)
end

# Test tilde path completion
let (c, r, res) = test_complete("\"~/julia")
    if !Sys.iswindows()
        @test res && c == String[homedir() * "/julia"]
    else
        @test !res
    end

    c, r, res = test_complete("\"foo~bar")
    @test !res
end

# Test the completion returns nothing when the folder do not exist
let (c, r) = test_complete("cd(\"folder_do_not_exist_77/file")
    @test length(c) == 0
end

if Sys.iswindows()
    tmp = tempname()
    touch(tmp)
    path = realpath(dirname(tmp))
    file = basename(tmp)
    temp_name = basename(path)
    cd(path) do
        s = "cd ..\\\\"
        c,r = test_scomplete(s)
        @test r == length(s)+1:length(s)
        @test temp_name * "\\\\" in c

        s = "ls $(file[1:2])"
        c,r = test_scomplete(s)
        @test r == length(s)-1:length(s)
        @test file in c

        s = "cd(\"..\\"
        c,r = test_complete(s)
        @test r == length(s)+1:length(s)
        @test temp_name * "\\\\" in c

        s = "cd(\"$(file[1:2])"
        c,r = test_complete(s)
        @test r == length(s) - 1:length(s)
        @test (length(c) > 1 && file in c) || (["$file\""] == c)
    end
    rm(tmp)
end

# auto completions of true and false... issue #14101
let s = "tru"
    c, r, res = test_complete(s)
    @test "true" in c
end

let s = "fals"
    c, r, res = test_complete(s)
    @test "false" in c
end

# Don't crash when attempting to complete a tuple, #15329
let s = "CompletionFoo.tuple."
    c, r, res = test_complete(s)
    @test isempty(c)
end

@testset "sub/superscripts" begin
    @test "⁽¹²³⁾ⁿ" in test_complete("\\^(123)n")[1]
    @test "ⁿ" in test_complete("\\^n")[1]
    @test "ᵞ" in test_complete("\\^gamma")[1]
    @test isempty(test_complete("\\^(123)nq")[1])
    @test "₍₁₂₃₎ₙ" in test_complete("\\_(123)n")[1]
    @test "ₙ" in test_complete("\\_n")[1]
    @test "ᵧ" in test_complete("\\_gamma")[1]
    @test isempty(test_complete("\\_(123)nq")[1])
end

# test Dicts
function test_dict_completion(dict_name)
    s = "$dict_name[\"ab"
    c, r = test_complete(s)
    @test c == Any["\"abc\"", "\"abcd\""]
    s = "$dict_name[\"abcd"
    c, r = test_complete(s)
    @test c == Any["\"abcd\"]"]
    s = "$dict_name[ \"abcd"  # leading whitespace
    c, r = test_complete(s)
    @test c == Any["\"abcd\"]"]
    s = "$dict_name[\"abcd]"  # trailing close bracket
    c, r = completions(s, lastindex(s) - 1)
    c = map(completion_text, c)
    @test c == Any["\"abcd\""]
    s = "$dict_name[:b"
    c, r = test_complete(s)
    @test c == Any[":bar", ":bar2"]
    s = "$dict_name[:bar2"
    c, r = test_complete(s)
    @test c == Any[":bar2]"]
    s = "$dict_name[Ba"
    c, r = test_complete(s)
    @test c == Any["Base]"]
    s = "$dict_name[occ"
    c, r = test_complete(s)
    @test c == Any["occursin]"]
    s = "$dict_name[`l"
    c, r = test_complete(s)
    @test c == Any["`ls`]"]
    s = "$dict_name[6"
    c, r = test_complete(s)
    @test c == Any["66", "67"]
    s = "$dict_name[66"
    c, r = test_complete(s)
    @test c == Any["66]"]
    s = "$dict_name[("
    c, r = test_complete(s)
    @test c == Any["(\"q\", 3)]"]
    s = "$dict_name[\"\\alp"
    c, r = test_complete(s)
    @test c == String["\\alpha"]
    s = "$dict_name[\"\\alpha"
    c, r = test_complete(s)
    @test c == String["α"]
    s = "$dict_name[\"α"
    c, r = test_complete(s)
    @test c == Any["\"α\"]"]
    s = "$dict_name[:\\alp"
    c, r = test_complete(s)
    @test c == String["\\alpha"]
    s = "$dict_name[:\\alpha"
    c, r = test_complete(s)
    @test c == String["α"]
    s = "$dict_name[:α"
    c, r = test_complete(s)
    @test c == Any[":α]"]
    s = "$dict_name["
    c, r = test_complete(s)
    @test c == sort!(repr.(keys(Main.CompletionFoo.test_dict)))
end
test_dict_completion("CompletionFoo.test_dict")
test_dict_completion("CompletionFoo.test_customdict")
test_dict_completion("test_repl_comp_dict")
test_dict_completion("test_repl_comp_customdict")

# Issue #23004: this should not throw:
@test REPLCompletions.dict_identifier_key("test_dict_ℂ[\\", :other) isa Tuple

@testset "completion of string/cmd macros (#22577)" begin
    c, r, res = test_complete("ra")
    @test "raw\"" in c
    c, r, res = test_complete("CompletionFoo.tests")
    @test "teststr\"" in c
    c, r, res = test_complete("CompletionFoo.tϵsτs")
    @test "tϵsτstρ\"" in c
    c, r, res = test_complete("CompletionFoo.testc")
    @test "testcmd`" in c
    c, r, res = test_complete("CompletionFoo.tϵsτc")
    @test "tϵsτcmδ`" in c
end

# Test completion in context

# No CompletionFoo.CompletionFoo
let s = ""
    c, r = test_complete_foo(s)
    @test !("CompletionFoo" in c)
end

# Can see `rand()` after `using Random`
let s = "r"
    c, r = test_complete_foo(s)
    @test "rand" in c
    @test r == 1:1
    @test s[r] == "r"
end

# Can see `Test.AbstractTestSet` after `import Test`
let s = "Test.A"
    c, r = test_complete_foo(s)
    @test "AbstractTestSet" in c
    @test r == 6:6
    @test s[r] == "A"
end

# Can complete relative import
let s = "import ..M"
    c, r = test_complete_foo(s)
    @test_broken "Main" in c
    @test r == 10:10
    @test s[r] == "M"
end

let s = ""
    c, r = test_complete_foo(s)
    @test "bar" in c
    @test r === 1:0
    @test s[r] == ""
end

let s = "f"
    c, r = test_complete_foo(s)
    @test "foo" in c
    @test r == 1:1
    @test s[r] == "f"
    @test !("foobar" in c)
end

let s = "@f"
    c, r = test_complete_foo(s)
    @test "@foobar" in c
    @test r == 1:2
    @test s[r] == "@f"
    @test !("foo" in c)
end

let s = "type_test.x"
    c, r = test_complete_foo(s)
    @test "xx" in c
    @test r == 11:11
    @test s[r] == "x"
end

let s = "bar.no_val_available"
    c, r = test_complete_foo(s)
    @test length(c)==0
end

let s = "type_test.xx.y"
    c, r = test_complete_foo(s)
    @test "yy" in c
    @test r == 14:14
    @test s[r] == "y"
end

let s = ":(function foo(::Int) end).args[1].args[2]."
    c, r = test_complete_foo(s)
    @test c == Any[]
end

let s = "log(log.(x),"
    c, r = test_complete_foo(s)
    @test !isempty(c)
end

let s = "Base.return_types(getin"
    c, r = test_complete_foo(s)
    @test "getindex" in c
    @test r == 19:23
    @test s[r] == "getin"
end

let s = "using Test, Random"
    c, r = test_complete_foo(s)
    @test !("RandomDevice" in c)
end

let s = "test(1,1, "
    c, r, res = test_complete_foo(s)
    @test !res
    @test c[1] == string(first(methods(Main.CompletionFoo.test, Tuple{Int, Int})))
    @test c[2] == string(first(methods(Main.CompletionFoo.test, Tuple{})))  # corresponding to the vararg
    @test length(c) == 2
    # In particular, this checks that test(x::Real, y::Real) is not a valid completion
    # since it is strictly less specific than test(x::T, y::T) where T
    @test r == 1:4
    @test s[r] == "test"
end

let s = "test.(1,1, "
    c, r, res = test_complete_foo(s)
    @test !res
    @test length(c) == 4
    @test r == 1:4
    @test s[r] == "test"
end

let s = "prevind(\"θ\",1,"
    c, r, res = test_complete_foo(s)
    @test c[1] == string(first(methods(prevind, Tuple{String, Int})))
    @test r == 1:7
    @test s[r] == "prevind"
end

# Issue #32840
let s = "typeof(+)."
    c, r = test_complete_foo(s)
    @test length(c) == length(fieldnames(DataType))
end

let s = "test_dict[\"ab"
    c, r = test_complete_foo(s)
    @test c == Any["\"abc\"", "\"abcd\""]
end

let s = "CompletionFoo.x."
    c, r = test_complete(s)
    @test "a" in c
end

# https://github.com/JuliaLang/julia/issues/27184
let
    (test_complete("@noexist."); @test true)
    (test_complete("Main.@noexist."); @test true)
    (test_complete("@Main.noexist."); @test true)
end

let # Check that completion does not crash on (possibly invalid) macro calls
    (test_complete("@show."); @test true)
    (test_complete("@macroexpand."); @test true)
    (test_complete("@.."); @test true)
    (test_complete("CompletionFoo.@foobar."); @test true)
    (test_complete("CompletionFoo.@foobar()."); @test true)
    (test_complete("CompletionFoo.@foobar(4)."); @test true)
    (test_complete("CompletionFoo.@barfoo."); @test true)
    (test_complete("CompletionFoo.@barfoo()."); @test true)
    (test_complete("CompletionFoo.@barfoo(6)."); @test true)
    (test_complete("CompletionFoo.@error_expanding."); @test true)
    (test_complete("CompletionFoo.@error_expanding()."); @test true)
    (test_complete("CompletionFoo.@error_lowering_conditional."); @test true)
    (test_complete("CompletionFoo.@error_lowering_conditional()."); @test true)
    (test_complete("CompletionFoo.@error_lowering_conditional(3)."); @test true)
    (test_complete("CompletionFoo.@error_lowering_conditional('a')."); @test true)
    (test_complete("CompletionFoo.@error_throwing."); @test true)
    (test_complete("CompletionFoo.@error_throwing()."); @test true)
end

@testset "https://github.com/JuliaLang/julia/issues/40247" begin
    # getfield type completion can work for complicated expression

    let
        m = Module()
        @eval m begin
            struct Rs
                rs::Vector{Regex}
            end
            var = nothing
            function foo()
                global var = 1
                return Rs([r"foo"])
            end
        end

        c, r = test_complete_context("foo(#=#==#=##==#).rs[1].", m)
        @test m.var === nothing # getfield type completion should never execute `foo()`
        @test length(c) == fieldcount(Regex)
    end

    let
        m = Module()
        @eval m begin
            struct R
                r::Regex
            end
            var = nothing
            function foo()
                global var = 1
                return R(r"foo")
            end
        end

        c, r = test_complete_context("foo().r.", m)
        @test m.var === nothing # getfield type completion should never execute `foo()`
        @test length(c) == fieldcount(Regex)

        c, r = test_complete_context("foo(#=#=# =#= =#).r.", m)
        @test m.var === nothing # getfield type completion should never execute `foo()`
        @test length(c) == fieldcount(Regex)
    end
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
"""
It return the `GenericRefFE` type for a Nedelec reference FE with order `order`
on the reference polytope `p`, using as scalar value `et`, i.e., `et` can be
`Float64` or `Float32`.
"""
# @santiagobadia : Project, go to complex numbers
function NedelecRefFE(p:: Polytope, et, order::Int)

  if !(all(extrusion(p).array .== HEX_AXIS))
    @notimplemented
  end

  # 1. Prebasis
  prebasis = GradMonomialBasis(VectorValue{dim(p),et},order)

  # Nface nodes, moments, and prebasis evaluated at nodes
  nf_nodes, nf_moments, pb_moments = _initialize_arrays(prebasis,p)

  # Face nodes and moments
  fcips, fmoments = _Nedelec_edge_values(p,et,order)

  # Insert nodes and moments in the nface arrays, together with the
  # prebasis evaluated at this nodes, for all edges
  nf_nodes,nf_moments,pb_moments = _insert_nface_values!(nf_nodes,nf_moments,pb_moments,prebasis,fcips,fmoments,p,1)

  # Idem for face values
  if ( dim(p) == 3 && order > 1)

    fcips, fmoments = _Nedelec_face_values(p,et,order)
    nf_nodes,nf_moments,pb_moments = _insert_nface_values!(nf_nodes,nf_moments,pb_moments,prebasis,fcips,fmoments,p,dim(p)-1)

  end

  # Idem for cell values
  if (order > 1)

    ccips, cmoments = _Nedelec_cell_values(p,et,order)
    nf_nodes,nf_moments,pb_moments = _insert_nface_values!(nf_nodes,nf_moments,pb_moments,prebasis,ccips,cmoments,p,dim(p))

  end

  _GenericRefFE(p,prebasis,nf_nodes,nf_moments,pb_moments)

end

# It provides for every edge the nodes and the moments arrays
function _Nedelec_edge_values(p,et,order)

  # Reference facet
  dim1 = 1
  fp = ref_nface_polytope(p,dim1) # dim 1

  # geomap from ref face to polytope faces
  fgeomap = _ref_face_to_faces_geomap(p,fp)

  # Compute integration points at all polynomial edges
  degree = order*2
  fquad = Quadrature(fp,degree)
  fips = coordinates(fquad)
  wips = weights(fquad)

  c_fips, fcips, fwips = _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)

  # Edge moments, i.e., M(Ei)_{ab} = q_RE^a(xgp_REi^b) w_Fi^b t_Ei ⋅ ()
  fshfs = Gridap.RefFEs._monomial_basis(fp,et,order-1)
  fmoments = _Nedelec_edge_moments(p, fshfs, c_fips, fcips, fwips)

  return fcips, fmoments

end

function _Nedelec_edge_moments(p, fshfs, c_fips, fcips, fwips)
  ts = edge_tangents(p)
  nc = length(c_fips)
  cfshfs = ConstantCellValue(fshfs, nc)
  cvals = evaluate(cfshfs,c_fips)
  cvals = [fwips[i]'.*cvals[i] for i in 1:nc]
  # @santiagobadia : Only working for oriented meshes now
  cvals = [ _broadcast(typeof(t),t,b) for (t,b) in zip(ts,cvals)]
  return cvals
end

# It provides for every face the nodes and the moments arrays
function _Nedelec_face_values(p,et,order)

  # Reference facet
  dimf = dim(p)-1
  fp = ref_nface_polytope(p,dimf) # dim 1

  # geomap from ref face to polytope faces
  fgeomap = _ref_face_to_faces_geomap(p,fp)

  # Compute integration points at all polynomial edges
  degree = order*2
  fquad = Quadrature(fp,degree)
  fips = coordinates(fquad)
  wips = weights(fquad)

  c_fips, fcips, fwips = _nfaces_evaluation_points_weights(p, fgeomap, fips, wips)

  # Face moments, i.e., M(Fi)_{ab} = w_Fi^b q_RF^a(xgp_RFi^b) (n_Fi × ())
  fshfs = GradMonomialBasis(VectorValue{dim(fp),et},order-1)
  fmoments = _Nedelec_face_moments(p, fshfs, c_fips, fcips, fwips)

  return fcips, fmoments

end

function _Nedelec_face_moments(p, fshfs, c_fips, fcips, fwips)
  nc = length(c_fips)
  fvs = nfaces_vertices(p,dim(p)-1)
  fts = [hcat([vs[2]-vs[1]...],[vs[3]-vs[1]...]) for vs in fvs]
  cfshfs = ConstantCellValue(fshfs, nc)
  # Ref facet FE functions evaluated at the facet integration points (in ref facet)
  cvals = evaluate(cfshfs,c_fips)
  cvals = [fwips[i]'.*cvals[i] for i in 1:nc]
  fns, os = face_normals(p)
  # @santiagobadia : Temporary hack for making it work for structured hex meshes
  ft = eltype(fns)
  cvals = [ Gridap.RefFEs.GenericRefFEs._broadcast_extend(ft,Tm,b) for (Tm,b) in zip(fts,cvals)]
  cvals = [ Gridap.RefFEs.GenericRefFEs._broadcast_cross(ft,n*o,b) for (n,o,b) in zip(fns,os,cvals)]
  return cvals
end

# It provides for every cell the nodes and the moments arrays
function _Nedelec_cell_values(p,et,order)
  # Compute integration points at interior
  degree = 2*order
  iquad = Quadrature(p,degree)
  ccips = coordinates(iquad)
  cwips = weights(iquad)

  # Cell moments, i.e., M(C)_{ab} = q_C^a(xgp_C^b) w_C^b ⋅ ()
  cbasis = CurlGradMonomialBasis(VectorValue{dim(p),et},order-1)
  cmoments = _Nedelec_cell_moments(p, cbasis, ccips, cwips )

  return [ccips], [cmoments]

end

const _Nedelec_cell_moments = _RT_cell_moments

ext: jl


text:
using JuMP
import NLopt
using Ipopt

function D2(v0::Vector)

    function ExplicitFunction(uw::Vector, grad::Vector, v::Vector)
        println("called my own OFunction")
        if length(grad) > 0
            grad[1] = 0
            grad[2] = 0.5 / sqrt(x[2])
        end
        w1, w2, u1, u2, u3, u4 = uw
        return u3 * (u2 + u4 - u1 - u3 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
               u1 * (u2 + u4 - u1 - u3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
               u4 * (-4 + u1 + u3 - u2 - u4 + (u1 + u2 - u3 - u4 - 2w2 - 6w1) / (2w1 + 2w2)) +
               u2 * (-4 + u1 + u3 - u2 - u4 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2)) +
               w1 * (-3 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2 +
               (1 / 2) * (-u1 + u2 - u3 + u4)^2 +
               w2 * (-1 + (u3 + u4 - u1 - u2 + 2w2 + 6w1) / (2w1 + 2w2))^2 - w1 * v[1] - w2 * v[2]
    end
    function ImplicitFunction(uw::Vector, grad::Vector, v::Vector)

        if length(grad) > 0
            grad[1] = 0
            grad[2] = 0.5 / sqrt(x[2])
        end
        model2 = Model(Ipopt.Optimizer)
        JuMP.set_silent(model2)
        JuMP.@variable(model2, x1)
        JuMP.@variable(model2, x2)
        JuMP.@NLexpression(model2, my_expr, uw[1] * (x1 - 3)^2 + (x2 - 1)^2 + uw[2] * (x1 - 1)^2 + (x2 - 1)^2)
        JuMP.@NLobjective(model2, Min,
            uw[3] * (x1 + 2 * x2 - 2) +
            uw[4] * (x1 - 2 * x2 - 2) +
            uw[5] * (-x1 + 2 * x2 - 2) +
            uw[6] * (-x1 - 2 * x2 - 2) +
            my_expr
        )
        JuMP.optimize!(model2)
        return JuMP.objective_value(model2) - uw[1] * v[1] - uw[2] * v[2]
    end

    function constraint(uw::Vector, _grad::Vector)
        uw[1] + uw[2] - 1
    end

    opt = NLopt.Opt(:LN_COBYLA, 6) # LN_COBYLA LD_MMA
    opt.lower_bounds = repeat([0], 6)
    # opt.xtol_rel = 1e-2

    opt.max_objective = (x, g) -> ImplicitFunction(x, g, v0)
    NLopt.equality_constraint!(opt, constraint, 1e-8)

    (minf, minx, ret) = NLopt.optimize(opt, repeat([0.9], 6))
    numevals = opt.numevals # the number of function evaluations
    w1, w2, u1, u2, u3, u4 = minx
    return (value.([w1, w2]), value.([u1, u2, u3, u4]))
end


text:
include("utils/helper_functions.jl")
include("calibration/calibration_helper_functions.jl")
include("calibration/create_log_posterior.jl")
include("calibration/create_log_posterior_wide_priors.jl")

include("create_models/additional_model_components/rf_ch4_direct_fund.jl")
include("create_models/additional_model_components/rf_ch4_etminan.jl")
include("create_models/additional_model_components/rf_ch4_myhre_fair.jl")
include("create_models/additional_model_components/rf_ch4_total_fund.jl")
include("create_models/additional_model_components/rf_co2_etminan.jl")
include("create_models/additional_model_components/total_co2_emissions.jl")
include("create_models/additional_model_components/rf_total.jl")

include("create_models/create_iam_dice.jl")
include("create_models/create_iam_fund.jl")
include("create_models/create_sneasy_fairch4.jl")
include("create_models/create_sneasy_fundch4.jl")
include("create_models/create_sneasy_hectorch4.jl")
include("create_models/create_sneasy_magiccch4.jl")

include("calibration/run_climate_models/run_sneasy_fairch4.jl")
include("calibration/run_climate_models/run_sneasy_fundch4.jl")
include("calibration/run_climate_models/run_sneasy_hectorch4.jl")
include("calibration/run_climate_models/run_sneasy_magiccch4.jl")

include("climate_projections/sneasych4_baseline_case.jl")
include("climate_projections/sneasych4_outdated_forcing.jl")
include("climate_projections/sneasych4_remove_correlations.jl")
include("climate_projections/sneasych4_us_ecs.jl")

include("scch4/dice_scch4.jl")
include("scch4/fund_scch4.jl")


text:
@deprecate label_components(A::AbstractArray, region::Union{Dims, AbstractVector{Int}}, bkg = 0) label_components(A; bkg=bkg, dims=region)
@deprecate label_components(A::AbstractArray{T,N}, connectivity::Array{Bool,N}, bkg) where {T,N}   label_components(A, connectivity; bkg=bkg)
@deprecate label_components!(out::AbstractArray{Int}, A::AbstractArray, region::Union{Dims, AbstractVector{Int}}, bkg = 0)  label_components!(out, A; bkg=bkg, dims=region)
@deprecate label_components!(out::AbstractArray{Int,N}, A::AbstractArray{T,N}, connectivity::Array{Bool,N}, bkg) where {T,N}  label_components!(out, A, connectivity; bkg=bkg)

@deprecate imfill(img::AbstractArray{Bool}, interval::Tuple{Real,Real}, dims::Union{Dims, AbstractVector{Int}}) imfill(img, interval; dims=dims)

import .FeatureTransform: feature_transform
@deprecate feature_transform(img, weights; kwargs...) feature_transform(img; weights=weights, kwargs...)


text:
repo = joinpath(homedir(), "ReferenceImages", "gallery")
recordings = joinpath(@__DIR__, "test_recordings")
cp(recordings, repo, force = true)


text:
# The code is from
# http://thirld.com/blog/2015/05/30/julia-profiling-cheat-sheet/
function benchmark(func, args...; filename = "trace.def", delay = 0.01)
    # Any setup code goes here.

    # Run once, to force compilation.
    println("======================= First run:")
    @time func(args...)

    # Run a second time, with profiling.
    println("======================= Second run:")

    Profile.init(delay = delay)
    Profile.clear()
    Profile.clear_malloc_data()

    @profile @time func(args...)

    # Write profile results to profile.bin.
    f = open(filename, "w")
    Profile.print(f)
    close(f)
end


text:
using LabelledGraphs

export ising_graph, rank_vec, cluster, rank, nodes, basis_size, biases, couplings, IsingGraph, prune

const Instance = Union{String, Dict}

unique_nodes(ising_tuples) = sort(collect(Set(Iterators.flatten((i, j) for (i, j, _) ∈ ising_tuples))))

const IsingGraph = LabelledGraph{MetaGraph{Int64, Float64}, Int64}

"""
$(TYPEDSIGNATURES)

Create the Ising spin glass model.

# Details

Store extra information
"""
function ising_graph(
    instance::Instance;
    sgn::Number=1.0,
    rank_override::Dict{Int, Int}=Dict{Int, Int}()
)
    # load the Ising instance
    if instance isa String
        ising = CSV.File(instance, types = [Int, Int, Float64], header=0, comment = "#")
    else
        ising = [ (i, j, J) for ((i, j), J) ∈ instance ]
    end

    nodes = unique_nodes(ising)
    L = length(nodes)
    nodes_to_vertices = Dict(w => i for (i, w) ∈ enumerate(nodes))

    ig = LabelledGraph{MetaGraph}(nodes)

    set_prop!.(Ref(ig), vertices(ig), :h, 0)
    foreach(v -> set_prop!(ig, v, :rank, get(rank_override, v, 2)), vertices(ig))

    for (i, j, v) ∈ ising
        v *= sgn

        if i == j
            set_prop!(ig, i, :h, v)
        else
            add_edge!(ig, i, j) || throw(ArgumentError("Duplicate Egde ($i, $j)"))
            set_prop!(ig, i, j, :J, v)
        end
    end

    set_prop!(
        ig,
        :rank,
        Dict{Int, Int}(
            v => get(rank_override, v, 2) for v in vertices(ig)
        )
    )

    ig
end

rank_vec(ig::IsingGraph) = Int[get_prop((ig), v, :rank) for v ∈ vertices(ig)]
basis_size(ig::IsingGraph) = prod(prod(rank_vec(ig)))
biases(ig::IsingGraph) = get_prop.(Ref(ig), vertices(ig), :h)

function couplings(ig::IsingGraph)
    J = zeros(nv(ig), nv(ig))
    for edge in edges(ig)
        i, j = ig.reverse_label_map[src(edge)], ig.reverse_label_map[dst(edge)]
        J[i, j] = get_prop(ig, edge, :J)
    end
    J
end

cluster(ig::IsingGraph, verts) = induced_subgraph(ig, collect(verts))

function inter_cluster_edges(ig::IsingGraph, cl1::IsingGraph, cl2::IsingGraph)
    verts1, verts2 = vertices(cl1), vertices(cl2)

    outer_edges = [
        LabelledEdge(i, j)
        for i ∈ vertices(cl1), j ∈ vertices(cl2)
        if has_edge(ig, i, j)
    ]

    J = zeros(nv(cl1), nv(cl2))
    for e ∈ outer_edges
        i, j = cl1.reverse_label_map[src(e)], cl2.reverse_label_map[dst(e)]
        @inbounds J[i, j] = get_prop(ig, e, :J)
    end
    outer_edges, J
end

function prune(ig::IsingGraph)
    to_keep = vcat(
        findall(!iszero, degree(ig)),
        findall(x->iszero(degree(ig, x)) && !isapprox(get_prop(ig, x, :h), 0, atol=1e-14), vertices(ig))
    )

    gg = ig[ig.labels[to_keep]]
    labels = collect(vertices(gg.inner_graph))
    reverse_label_map = Dict(i => i for i=1:nv(gg.inner_graph))
    LabelledGraph(labels, gg.inner_graph, reverse_label_map)
end


text:
header("GRAPHICS")
@testset "graphics" begin
    str = """
    <graphics>
     <offset x="1" y="2" />
     <line  color="linecolor" shape="line" style="solid" width="1.0"/>
     <position  x="1" y="2" />
     <position  x="3" y="4" />
     <dimension x="5" y="6" />
     <offset    x="7" y="8" /><!-- override first offset -->
     <fill  color="fillcolor" gradient-color="none" gradient-rotation="horizontal"/>
     <font align="center" family="Dialog" rotation="0.0" size="11"
           style="normal" weight="normal" />
    </graphics>
    """
    doc = EzXML.parsexml(str)
    n = parse_node(root(doc); reg=PNML.IDRegistry())
    printnode(n)

    @test n.offset isa PNML.Coordinate
    @test n.dimension isa PNML.Coordinate
    @test n.position isa Vector{PNML.Coordinate}

    @test_call PNML.Coordinate(1,2)
    @test_call PNML.Coordinate(1.1,2.2)

    # There can only be one offset, last tag parsed wins.
    @test n.offset == PNML.Coordinate(7,8)
    @test n.dimension == PNML.Coordinate(5,6)
    @test length(n.position) == 2
    @test n.position == [PNML.Coordinate(1,2), PNML.Coordinate(3,4)]

    @test n.line isa PNML.Line
    @test n.line.color == "linecolor"
    @test n.line.shape == "line"
    @test n.line.style == "solid"
    @test n.line.width == "1.0"

    @test n.fill isa PNML.Fill
    @test n.fill.color == "fillcolor"
    @test n.fill.image === nothing
    @test n.fill.gradient_color == "none"
    @test n.fill.gradient_rotation === "horizontal"

    @test n.font isa PNML.Font
    @test n.font.family == "Dialog"
    @test n.font.style == "normal"
    @test n.font.weight == "normal"
    @test n.font.size == "11"
    @test n.font.decoration === nothing
    @test n.font.align == "center"
    @test n.font.rotation == "0.0"
end

@testset "tokengraphics" begin
    str0 = """
 <tokengraphics>
 </tokengraphics>
"""
    str1 = """
 <tokengraphics>
     <tokenposition x="-9" y="-2"/>
 </tokengraphics>
"""
    str2 = """
 <tokengraphics>
     <tokenposition x="-9" y="-2"/>
     <tokenposition x="2"  y="3"/>
 </tokengraphics>
"""
    str3 = """
 <tokengraphics>
     <tokenposition x="-9" y="-2"/>
     <tokenposition x="2"  y="3"/>
     <tokenposition x="-2" y="2"/>
 </tokengraphics>
"""
    str4 = """
 <tokengraphics>
     <tokenposition x="-9" y="-2"/>
     <tokenposition x="2"  y="3"/>
     <tokenposition x="-2" y="2"/>
     <tokenposition x="-2" y="-22"/>
 </tokengraphics>
"""
    str5 = """
 <tokengraphics>
     <tokenposition x="-1.2" y="-22.33"/>
 </tokengraphics>
"""
    str6 = """
 <tokengraphics>
     <tokenposition x="-1" y="-2"/>
     <tokenposition x="-1.2" y="-22.33"/>
 </tokengraphics>
"""
    @testset "tokengraphics $l tokenpositions" for (s,l) in [str0=>0,
                                                             str1=>1,
                                                             str2=>2,
                                                             str3=>3,
                                                             str4=>4,
                                                             str5=>1,
                                                             str6=>2]
        n = parse_node(root(EzXML.parsexml(s)); reg=PNML.IDRegistry())
        @test n isa PNML.TokenGraphics
        @test length(n.positions) == l
        #for p in n.positions
        #    @show p, PNML.eltype(p.x)
        #end
    end
end


text:
using HTTP, JSON

todos = Dict(     # type is Dict{Int64, String}
    1 => "Getting groceries",
    2 => "Visiting my therapist",
    3 => "Getting a haircut"
)

json_string = JSON.json(todos)      # 1
# "{\"2\":\"Visiting my therapist\",\"3\":\"Getting a haircut\",\"1\":\"Getting groceries\"}"
todos2 = JSON.parse(json_string)
# Dict{String, Any} with 3 entries:
#   "1" => "Getting groceries"
#   "2" => "Visiting my therapist"
#   "3" => "Getting a haircut"

resp = HTTP.Response(
    200,
    ["Content-Type" => "application/json"],
    body=JSON.json(todos)
)

# HTTP.Messages.Response:
# """
# HTTP/1.1 200 OK
# Content-Type: application/json

# {"2":"Visiting my therapist","3":"Getting a haircut","1":"Getting groceries"}"""

body = HTTP.payload(resp) # this gives a Vector of UInt8 bytes
io = IOBuffer(body)
todos = JSON.parse(io)
# Dict{String, Any} with 3 entries:
#   "1" => "Getting groceries"
#   "2" => "Visiting my therapist"
#   "3" => "Getting a haircut"


text:
function _isavailable_xclip()
    if !success(pipeline(`which xclip`))
        error("Please install xclip to your system")
    end
end

"""
Paste an image from clipboard using xclip
"""
function _xclip()
    _isavailable_xclip()
    img_buf = IOBuffer()

    # Pipe clipboard image to buffer
    if success(pipeline(`xclip -selection clipboard -t image/png -o`; stdout=img_buf))
        # Load image from buffer
        img = load(img_buf)
        return img
    else
        error("No image in clipboard")
    end
end

"""
Copy an image to clipboard using xclip
"""
function _xclip(img::Matrix{<:Colorant})
    _isavailable_xclip()
    img_buf = IOBuffer()

    # Save given image to buffer
    save(Stream{format"PNG"}(img_buf), img)

    # Copy to clipboard
    open(`xclip -selection clipboard -t image/png`, "w", stdout) do f
        print(f, String(take!(img_buf)))
    end
end


text:
InputFolder = './Images/RiceGrains/';
OutputFolder = './Results/Images/RiceGrains/';

@iA = '*.tif';

@fxm_mModBinWat [iA] > [M];
params.SmallHolesArea = 3;
params.DistStdRad = 3;
params.MinDistLocVar = 5;
/endf

/show iA > M;
/keep M > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
const textmime = MIME"text/plain"()

@testset "display" begin
    o1 = SetUnion([2,3])
    @test Display.hasstencil(o1)
    @test Display.displaystencil(o1) == Any[2, " ∪ ", 3]
    @test repr(o1) == "SetUnion([2, 3])"
    io = IOBuffer()
    show(io, textmime, o1)
    @test String(take!(io)) == "2 ∪ 3"

    A1 = 1:4
    A2 = 5:10
    A3 = 11:20
    c1 = CompositeObject(A1)
    c2 = CompositeObject(A2)
    c3 = CompositeObject(A3)
    @test repr(c1) == "CompositeObject(1:4)"
    show(io, textmime, c1)
    @test String(take!(io)) == "CompositeObject(1, 2, 3, 4)"
    o2 = SetUnion((c1, c2, c3))
    @test Display.hasstencil(o2)
    @test Display.displaystencil(o2) == Any[c1, " ∪ ", c2, " ∪ ", c3]
    @test repr(o2) == "SetUnion((CompositeObject(1:4), CompositeObject(5:10), CompositeObject(11:20)))"

    show(io, textmime, o2)
    @test String(take!(io))[1:32] == "CompositeObject(1, 2, 3, 4) ∪ "

    o3 = SetUnion( (SetUnion((c1,c2)),c3) )
    @test repr(o3) == "SetUnion((SetUnion((CompositeObject(1:4), CompositeObject(5:10))), CompositeObject(11:20)))"
    show(io, textmime, o3)
    @test String(take!(io))[1:6] == "S ∪ "

    m = StencilObject(cos, c1)
    @test Display.hasstencil(m)
    @test Display.displaystencil(m) == [Display.SymbolObject(cos), '(', c1, ')']
    @test repr(m) == "StencilObject(cos, CompositeObject(1:4))"
    show(io, textmime, m)
    @test String(take!(io)) == "F(CompositeObject(1, 2, 3, 4))\n\nF = cos"

    m1 = StencilObject(cos, c1)
    m2 = StencilObject(cos, c2)
    m3 = StencilObject(sin, c3)
    C = SetUnion((m1,m2,m3))
    @test repr(C) == "SetUnion((StencilObject(cos, CompositeObject(1:4)), StencilObject(cos, CompositeObject(5:10)), StencilObject(sin, CompositeObject(11:20))))"
    show(io, textmime, C)
    @test String(take!(io))[1] == 'F'

    textc = CompositeObject(["abc", "def"])
    show(io, textmime, textc)
    @test String(take!(io)) == "CompositeObject(\"abc\", \"def\")"
end


text:
function oneDArray(t::DataType, len::Int)
    return zeros(t, len)
end

function twoDArray(t::DataType, len_x::Int, len_y::Int)
    return zeros(t, len_x, len_y)
end


text:
module ADT

using MLStyle.Err
using MLStyle.Private
import MLStyle.Match: PatternDef, pattern_match

export @case, @data

function _check_components(arg)
    if !(isa(arg, Expr) && arg.head == :(::) ||
         isa(arg, Symbol)                    ||
         isa(arg, LineNumberNode))
        SyntaxError("`$arg`")
    end
end

macro case(cons)

    ty = nothing

    if isa(cons, Expr)
        if cons.head == :<:
            cons, ty = cons.args
        end
    end

    let cons =
        if isa(cons, Expr) && cons.head == :macrocall
                macroexpand(
                    __module__,
                    cons,
                    recursive=False) |> esc
        else
            cons
        end

    let args =
        if isa(cons, Expr) && cons.head == :call
            cons.args
        else
            [cons]
        end

    foreach(_check_components, args)

    let head = args[1], tail = args[2:end]

        # make struct
        if ty === nothing
            quote
                struct ($head)
                    $(tail...)
                end
            end

        else
            quote
                struct ($head) <: ($ty)
                    $(tail...)
                end
            end
        end |>
        function (ast)
            @eval __module__ $(ast)
        end

        most_union_all = get_most_union_all(head, __module__)

        # get fieldnames
        fields = fieldnames(most_union_all)

        PatternDef.App(most_union_all)  do  args, guard, tag, mod
            if length(args) != length(fields)
                DataTypeUsageError("Got patterns `$(repr(args))`, expected: `$fields`") |> throw
            end

            map(zip(fields, args)) do (field, arg)
                pattern_match(arg, nothing, :($tag.$field), mod)
            end |>
            function (last)
                reduce((a, b) -> Expr(:&&, a, b), last, init=:($isa($tag, $most_union_all)))
            end |>
            function (last)
                if guard === nothing
                    last
                else
                    :($last && $guard)
                end
            end
        end
        nothing

    end
    end
    end
end

macro data(abs_ty, cases)

    if !isa(cases, Expr) || !(cases.head in (:block, :bracescat, :braces))
        SyntaxError("$(repr(cases)).") |> throw
    end

    @eval __module__ abstract type $abs_ty end

    map(cases.args) do case
        if isa(case, LineNumberNode)
            case
        else
            :(@case $case <: $abs_ty)
        end
    end |>
    function (last)
        quote
            $(last...)
        end |> esc
    end

end

end


text:
module CodingTheory

include("linearcode.jl")

# export everything needed

end # module

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
"Parameters for right-hand side function of implicit partitioned Runge-Kutta methods."
mutable struct ParametersPRKimplicit{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::PartitionedTableau{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersPRKimplicit{DT,D}(equs::ET, tab::PartitionedTableau{TT}, Δt::TT) where {D, DT, TT, ET <: NamedTuple}
        new{DT, TT, D, tab.s, ET}(equs, tab, Δt, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

@doc raw"""
Implicit partitioned Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of q
* `p̃`: initial guess of p
* `ṽ`: initial guess of v
* `f̃`: initial guess of f
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of q
* `P`: internal stages of p
* `V`: internal stages of v
* `F`: internal stages of f
* `Y`: vector field of internal stages of q
* `Z`: vector field of internal stages of p
"""
struct IntegratorCachePRKimplicit{ST,D,S} <: PODEIntegratorCache{ST,D}
    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    Q::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    F::Vector{Vector{ST}}
    Y::Vector{Vector{ST}}
    Z::Vector{Vector{ST}}

    function IntegratorCachePRKimplicit{ST,D,S}() where {ST,D,S}
        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        Q = create_internal_stage_vector(ST, D, S)
        P = create_internal_stage_vector(ST, D, S)
        V = create_internal_stage_vector(ST, D, S)
        F = create_internal_stage_vector(ST, D, S)
        Y = create_internal_stage_vector(ST, D, S)
        Z = create_internal_stage_vector(ST, D, S)

        new(q̃, p̃, ṽ, f̃, s̃, Q, P, V, F, Y, Z)
    end
end

function IntegratorCache(params::ParametersPRKimplicit{DT,TT,D,S}; kwargs...) where {DT,TT,D,S}
    IntegratorCachePRKimplicit{DT,D,S}(; kwargs...)
end

function IntegratorCache{ST}(params::ParametersPRKimplicit{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCachePRKimplicit{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersPRKimplicit{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCachePRKimplicit{ST,D,S}

@doc raw"""
Implicit partitioned Runge-Kutta integrator solving the system
```math
\begin{aligned}
P_{n,i} &= \vartheta (Q_{n,i}, V_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= f (Q_{n,i}, V_{n,i}) , &
P_{n,i} &= p_{n} + h  \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} ,
\end{aligned}
```
Usually we are interested in Lagrangian systems, where
```math
\begin{aligned}
P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , &
F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) ,
\end{aligned}
```
and tableaus satisfying the symplecticity conditions
```math
\begin{aligned}
b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , &
\bar{b}_i &= b_i .
\end{aligned}
```
"""
struct IntegratorPRKimplicit{DT, TT, D, S, PT <: ParametersPRKimplicit{DT,TT},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessIODE{TT}} <: AbstractIntegratorIRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorPRKimplicit(params::ParametersPRKimplicit{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorPRKimplicit{DT,D}(equations::NamedTuple, tableau::PartitionedTableau{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        # create params
        params = ParametersPRKimplicit{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, D*S, params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorPRKimplicit(params, solver, iguess, caches)
    end

    function IntegratorPRKimplicit{DT,D}(v::Function, f::Function, tableau::PartitionedTableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorPRKimplicit{DT,D}(NamedTuple{(:v,:f)}((v,f)), tableau, Δt; kwargs...)
    end

    function IntegratorPRKimplicit{DT,D}(v::Function, f::Function, h::Function, tableau::PartitionedTableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorPRKimplicit{DT,D}(NamedTuple{(:v,:f,:h)}((v,f,h)), tableau, Δt; kwargs...)
    end

    function IntegratorPRKimplicit(equation::Union{IODE{DT}, LODE{DT}}, tableau::PartitionedTableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorPRKimplicit{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorPRKimplicit{DT,TT,D,S}) where {DT,TT,D,S} = D

function initialize!(int::IntegratorPRKimplicit, sol::AtomicSolutionPODE)
    sol.t̄ = sol.t - timestep(int)

    equations(int)[:v̄](sol.t, sol.q, sol.v)
    equations(int)[:f̄](sol.t, sol.q, sol.v, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̄, sol.q̄, sol.p̄, sol.v̄, sol.f̄)
end

function update_params!(int::IntegratorPRKimplicit, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.p .= sol.p
end

function initial_guess!(int::IntegratorPRKimplicit{DT}, sol::AtomicSolutionPODE{DT},
                        cache::IntegratorCachePRKimplicit{DT}=int.caches[DT]) where {DT}

    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.Q[i], cache.P[i], cache.V[i], cache.F[i],
                              tableau(int).q.c[i], tableau(int).p.c[i])
    end
    for i in eachstage(int)
        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = 0
            for j in eachstage(int)
                int.solver.x[ndims(int)*(i-1)+k] += tableau(int).q.a[i,j] * cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                        P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        params::ParametersPRKimplicit{DT,TT,D,S}) where {ST,DT,TT,D,S}
    local tqᵢ::TT
    local tpᵢ::TT
    local y1::ST
    local y2::ST

    for i in 1:S
        for k in 1:D
            # copy x to V
            V[i][k] = x[D*(i-1)+k]

            # compute Q
            y1 = 0
            y2 = 0
            for j in 1:S
                y1 += params.tab.q.a[i,j] * V[j][k]
                y2 += params.tab.q.â[i,j] * V[j][k]
            end
            Q[i][k] = params.q[k] + params.Δt * (y1 + y2)
        end

        # compute time of internal stage
        tqᵢ = params.t + params.Δt * params.tab.q.c[i]
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]

        # compute ϑ(Q,V) and f(Q,V)
        params.equs[:ϑ](tqᵢ, Q[i], V[i], P[i])
        params.equs[:f](tpᵢ, Q[i], V[i], F[i])
    end
end

# Compute stages of implicit partitioned Runge-Kutta methods.
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersPRKimplicit{DT,TT,D,S},
                          caches::CacheDict) where {ST,DT,TT,D,S}

    # temporary variables
    local z1::ST
    local z2::ST

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q, cache.V, cache.Y, cache.P, cache.F, cache.Z, params)

    # compute b = [P-p-AF]
    for i in 1:S
        for k in 1:D
            z1 = 0
            z2 = 0
            for j in 1:S
                z1 += params.tab.p.a[i,j] * cache.F[j][k]
                z2 += params.tab.p.â[i,j] * cache.F[j][k]
            end
            b[D*(i-1)+k] = cache.P[i][k] - params.p[k] - params.Δt * (z1 + z2)
        end
    end
end

function integrate_step!(int::IntegratorPRKimplicit{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCachePRKimplicit{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from atomic solution
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.Y,
                                  cache.P, cache.F, cache.Z, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
# Julia wrapper for header: /usr/include/scip/scip_branch.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeBranchrule(scip, name, desc, priority, maxdepth, maxbounddist, branchcopy, branchfree, branchinit, branchexit, branchinitsol, branchexitsol, branchexeclp, branchexecext, branchexecps, branchruledata)
    ccall((:SCIPincludeBranchrule, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, Cdouble, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_BRANCHRULEDATA}), scip, name, desc, priority, maxdepth, maxbounddist, branchcopy, branchfree, branchinit, branchexit, branchinitsol, branchexitsol, branchexeclp, branchexecext, branchexecps, branchruledata)
end

function SCIPincludeBranchruleBasic(scip, branchruleptr, name, desc, priority, maxdepth, maxbounddist, branchruledata)
    ccall((:SCIPincludeBranchruleBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_BRANCHRULE}}, Cstring, Cstring, Cint, Cint, Cdouble, Ptr{SCIP_BRANCHRULEDATA}), scip, branchruleptr, name, desc, priority, maxdepth, maxbounddist, branchruledata)
end

function SCIPsetBranchruleCopy(scip, branchrule, branchcopy)
    ccall((:SCIPsetBranchruleCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchcopy)
end

function SCIPsetBranchruleFree(scip, branchrule, branchfree)
    ccall((:SCIPsetBranchruleFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchfree)
end

function SCIPsetBranchruleInit(scip, branchrule, branchinit)
    ccall((:SCIPsetBranchruleInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchinit)
end

function SCIPsetBranchruleExit(scip, branchrule, branchexit)
    ccall((:SCIPsetBranchruleExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchexit)
end

function SCIPsetBranchruleInitsol(scip, branchrule, branchinitsol)
    ccall((:SCIPsetBranchruleInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchinitsol)
end

function SCIPsetBranchruleExitsol(scip, branchrule, branchexitsol)
    ccall((:SCIPsetBranchruleExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchexitsol)
end

function SCIPsetBranchruleExecLp(scip, branchrule, branchexeclp)
    ccall((:SCIPsetBranchruleExecLp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchexeclp)
end

function SCIPsetBranchruleExecExt(scip, branchrule, branchexecext)
    ccall((:SCIPsetBranchruleExecExt, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchexecext)
end

function SCIPsetBranchruleExecPs(scip, branchrule, branchexecps)
    ccall((:SCIPsetBranchruleExecPs, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Ptr{Cvoid}), scip, branchrule, branchexecps)
end

function SCIPfindBranchrule(scip, name)
    ccall((:SCIPfindBranchrule, libscip), Ptr{SCIP_BRANCHRULE}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetBranchrules(scip)
    ccall((:SCIPgetBranchrules, libscip), Ptr{Ptr{SCIP_BRANCHRULE}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNBranchrules(scip)
    ccall((:SCIPgetNBranchrules, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetBranchrulePriority(scip, branchrule, priority)
    ccall((:SCIPsetBranchrulePriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Cint), scip, branchrule, priority)
end

function SCIPsetBranchruleMaxdepth(scip, branchrule, maxdepth)
    ccall((:SCIPsetBranchruleMaxdepth, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Cint), scip, branchrule, maxdepth)
end

function SCIPsetBranchruleMaxbounddist(scip, branchrule, maxbounddist)
    ccall((:SCIPsetBranchruleMaxbounddist, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_BRANCHRULE}, Cdouble), scip, branchrule, maxbounddist)
end

function SCIPgetLPBranchCands(scip, lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars)
    ccall((:SCIPgetLPBranchCands, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), scip, lpcands, lpcandssol, lpcandsfrac, nlpcands, npriolpcands, nfracimplvars)
end

function SCIPgetNLPBranchCands(scip)
    ccall((:SCIPgetNLPBranchCands, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioLPBranchCands(scip)
    ccall((:SCIPgetNPrioLPBranchCands, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetExternBranchCands(scip, externcands, externcandssol, externcandsscore, nexterncands, nprioexterncands, nprioexternbins, nprioexternints, nprioexternimpls)
    ccall((:SCIPgetExternBranchCands, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), scip, externcands, externcandssol, externcandsscore, nexterncands, nprioexterncands, nprioexternbins, nprioexternints, nprioexternimpls)
end

function SCIPgetNExternBranchCands(scip)
    ccall((:SCIPgetNExternBranchCands, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioExternBranchCands(scip)
    ccall((:SCIPgetNPrioExternBranchCands, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioExternBranchBins(scip)
    ccall((:SCIPgetNPrioExternBranchBins, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioExternBranchInts(scip)
    ccall((:SCIPgetNPrioExternBranchInts, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioExternBranchImpls(scip)
    ccall((:SCIPgetNPrioExternBranchImpls, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioExternBranchConts(scip)
    ccall((:SCIPgetNPrioExternBranchConts, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPaddExternBranchCand(scip, var, score, solval)
    ccall((:SCIPaddExternBranchCand, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble), scip, var, score, solval)
end

function SCIPclearExternBranchCands(scip)
    ccall((:SCIPclearExternBranchCands, libscip), Cvoid, (Ptr{SCIP_},), scip)
end

function SCIPcontainsExternBranchCand(scip, var)
    ccall((:SCIPcontainsExternBranchCand, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetPseudoBranchCands(scip, pseudocands, npseudocands, npriopseudocands)
    ccall((:SCIPgetPseudoBranchCands, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{Ptr{SCIP_VAR}}}, Ptr{Cint}, Ptr{Cint}), scip, pseudocands, npseudocands, npriopseudocands)
end

function SCIPgetNPseudoBranchCands(scip)
    ccall((:SCIPgetNPseudoBranchCands, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioPseudoBranchCands(scip)
    ccall((:SCIPgetNPrioPseudoBranchCands, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioPseudoBranchBins(scip)
    ccall((:SCIPgetNPrioPseudoBranchBins, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioPseudoBranchInts(scip)
    ccall((:SCIPgetNPrioPseudoBranchInts, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetNPrioPseudoBranchImpls(scip)
    ccall((:SCIPgetNPrioPseudoBranchImpls, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetBranchScore(scip, var, downgain, upgain)
    ccall((:SCIPgetBranchScore, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble), scip, var, downgain, upgain)
end

function SCIPgetBranchScoreMultiple(scip, var, nchildren, gains)
    ccall((:SCIPgetBranchScoreMultiple, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint, Ptr{Cdouble}), scip, var, nchildren, gains)
end

function SCIPgetBranchingPoint(scip, var, suggestion)
    ccall((:SCIPgetBranchingPoint, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, suggestion)
end

function SCIPcalcNodeselPriority(scip, var, branchdir, targetvalue)
    ccall((:SCIPcalcNodeselPriority, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BRANCHDIR, Cdouble), scip, var, branchdir, targetvalue)
end

function SCIPcalcChildEstimate(scip, var, targetvalue)
    ccall((:SCIPcalcChildEstimate, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble), scip, var, targetvalue)
end

function SCIPcreateChild(scip, node, nodeselprio, estimate)
    ccall((:SCIPcreateChild, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_NODE}}, Cdouble, Cdouble), scip, node, nodeselprio, estimate)
end

function SCIPbranchVar(scip, var, downchild, eqchild, upchild)
    ccall((:SCIPbranchVar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{Ptr{SCIP_NODE}}, Ptr{Ptr{SCIP_NODE}}, Ptr{Ptr{SCIP_NODE}}), scip, var, downchild, eqchild, upchild)
end

function SCIPbranchVarHole(scip, var, left, right, downchild, upchild)
    ccall((:SCIPbranchVarHole, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cdouble, Ptr{Ptr{SCIP_NODE}}, Ptr{Ptr{SCIP_NODE}}), scip, var, left, right, downchild, upchild)
end

function SCIPbranchVarVal(scip, var, val, downchild, eqchild, upchild)
    ccall((:SCIPbranchVarVal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Ptr{Ptr{SCIP_NODE}}, Ptr{Ptr{SCIP_NODE}}, Ptr{Ptr{SCIP_NODE}}), scip, var, val, downchild, eqchild, upchild)
end

function SCIPbranchVarValNary(scip, var, val, n, minwidth, widthfactor, nchildren)
    ccall((:SCIPbranchVarValNary, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cdouble, Cint, Cdouble, Cdouble, Ptr{Cint}), scip, var, val, n, minwidth, widthfactor, nchildren)
end

function SCIPbranchLP(scip, result)
    ccall((:SCIPbranchLP, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RESULT}), scip, result)
end

function SCIPbranchExtern(scip, result)
    ccall((:SCIPbranchExtern, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RESULT}), scip, result)
end

function SCIPbranchPseudo(scip, result)
    ccall((:SCIPbranchPseudo, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_RESULT}), scip, result)
end


text:
# furnsh("src/kernels/de440s.bsp")
# furnsh("src/kernels/naif0012.tls")

"""
    rot2inert(rv, θ, μ; origin=:barycenter)

Convert state vector `rv = [r; v]` {NON; NON} from rotating (synodic) frame to inertial
frame in the normalized CR3BP given time `θ` {NON} and mass parameter `μ` {NON}.
"""
function rot2inert(rv, θ, μ; origin=:barycenter)
    if origin == :barycenter
    elseif origin == :prim
        rv -= [-μ, 0, 0, 0, 0, 0]
    elseif origin == :sec
        rv -= [1-μ, 0, 0, 0, 0, 0]
    end
    cθ, sθ = cos(θ), sin(θ)
    A =  [cθ -sθ  0;
          sθ  cθ  0;
           0   0  1]
    B = [-sθ -cθ  0;
          cθ -sθ  0;
           0   0  0]
    C = [A zeros(3,3);
         B A]
    return C*rv
end

"""
    rot2inert(rv, θ, sys::System; origin=:barycenter)

Convert state vector `rv = [r; v]` {NON; NON} from rotating (synodic) frame to inertial
frame in the normalized CR3BP given time `θ` {NON} and CR3BP system `sys`.
"""
rot2inert(rv, θ, sys::System; origin=:barycenter) = rot2inert(rv, θ, sys.μ, origin=origin)

"""
    rot2inert(rv, θ, p::Array; origin=:barycenter)

Convert state vector `rv = [r; v]` {NON; NON} from rotating (synodic) frame to inertial
frame in the non-normalized CR3BP given time `θ` {NON} and p = [μ₁;μ₂;d]`
{km³/s²; km³/s²; km}, which contains the gravitational parameters of the first and second
primary bodies as well as the distance between them.
"""
function rot2inert(rv, θ, p::Array; origin=:barycenter)
    μ₁,μ₂,d = p # parameters

    if origin == :barycenter
    elseif origin == :prim
        rv -= [-d*μ₂/(μ₁ + μ₂), 0, 0, 0, 0, 0]
    elseif origin == :sec
        rv -= [d*(1-μ₂/(μ₁ + μ₂)), 0, 0, 0, 0, 0]
    end

    ωₛ = sqrt((μ₁ + μ₂)/d^3);
    cθ, sθ = cos(ωₛ*θ), sin(ωₛ*θ)
    A = [cθ -sθ  0;
         sθ  cθ  0;
          0   0  1]
    B = [-sθ -cθ  0;
          cθ -sθ  0;
           0   0  0]
    C = [A zeros(3,3);
         B A]
    return C*rv
end

"""
    rot2inert!(rv, θ, μ)

In-place version of rot2inert(rv, θ, μ)
"""
function rot2inert!(rv, θ, μ; origin=:barycenter)
    rv[1:6] = rot2inert(rv, θ, μ, origin=origin)
    return nothing
end

"""
    rot2inert!(rv, θ, sys::System)

In-place version of rot2inert(rv, θ, sys::System)
"""
function rot2inert!(rv, θ, sys::System; origin=:barycenter)
    rv[1:6] = rot2inert(rv, θ, sys, origin=origin)
    return nothing
end

"""
    rot2inert!(rv, θ, p::Array)

In-place version of rot2inert(rv, θ, p::Array)
"""
function rot2inert!(rv, θ, p::Array; origin=0)
    rv[1:6] = rot2inert(rv, θ, p, origin=origin)
    return nothing
end

"""
    inert2rot(rv, θ, μ)

Inertial frame to rotating (synodic) frame
"""
function inert2rot(rv, t, μ; origin=:barycenter)
    ct, st = cos(t), sin(t)
    A =  [ct  st  0;
         -st  ct  0;
           0   0  1]
    # A = rotz(-t)
    B = [-st  ct  0;
         -ct -st  0;
           0   0  0]
    C = [A zeros(3,3);
         B A]
    if origin == :barycenter
    elseif origin == :prim
        rv += [-μ, 0, 0, 0, 0, 0]
    elseif origin == :sec
        rv += [1-μ, 0, 0, 0, 0, 0]
    end
    return C*rv
end

"""
    inert2rot(rv, θ, sys::System)

Inertial frame to rotating (synodic) frame
"""
inert2rot(rv, t, sys::System; origin=:barycenter) = inert2rot(rv, t, sys.μ, origin=origin)

"""
    inert2rot(rv, θ, p::Array)

Inertial frame to rotating (synodic) frame
"""
function inert2rot(rv, t, p::Array; origin=:barycenter)
# function S2I!(rv,t,p::Array)
    μ₁,μ₂,d = p # parameters

    ωₛ = sqrt((μ₁ + μ₂)/d^3);
    ct, st = cos(ωₛ*t), sin(ωₛ*t)
    A = [ct  st  0;
        -st  ct  0;
          0   0  1]
    # A = rotz(ωₛ*t)
    B = [-st  ct  0;
         -ct -st  0;
           0   0  0]
    C = [A zeros(3,3);
         B A]

    if origin == :barycenter
    elseif origin == :prim
        rv += [-d*μ₂/(μ₁ + μ₂), 0, 0, 0, 0, 0]
    elseif origin == :sec
        rv += [d*(1-μ₂/(μ₁ + μ₂)), 0, 0, 0, 0, 0]
    end

    return C*rv
end

"""
    inert2rot!(rv, t, μ)

In-place version of `inert2rot(rv, t, μ)`.
"""
function inert2rot!(rv, t, μ; origin=:barycenter)
    # function inert2rot!(rv, t)
    rv[1:6] = inert2rot(rv, t, μ, origin=origin)
    return nothing
end

"""
    inert2rot!(rv, t, sys::System)

In-place version of `inert2rot(rv, t, sys::System)`.
"""
function inert2rot!(rv, t, sys::System; origin=:barycenter)
    rv[1:6] = inert2rot(rv, t, sys, origin=origin)
    return nothing
end

"""
    inert2rot!(rv, t, p::Array)

In-place version of `inert2rot(rv, t, p::Array)`.
"""
function inert2rot!(rv, t, p::Array; origin=0)
    rv[1:6] = inert2rot(rv, t, p, origin=origin)
    return nothing
end

"""
    ecef2eci(rv_ecef, θ; ω=7.292115373194e-5, ang_unit::Symbol=:rad)

Convert state vector `rv = [r; v]` {km; km/s} from Earth-Centered Earth-Fixed (ECEF) frame
to Earth-Centered Inertial (ECI) frame given Greenwich Mean Sidereal Time (GMST) `θ` {rad}.
"""
function ecef2eci(rv_ecef, θ; ω=7.292115373194e-5, ang_unit::Symbol=:rad)
    if ang_unit == :deg
        θ = deg2rad(θ)
    elseif ang_unit == :rad
    else
        error("ang_unit should be :rad or :deg")
    end

    R = rotz(θ)
    O = zeros(3,3)
    R_dot = ω*[-sin(θ) -cos(θ) 0;
                cos(θ) -sin(θ) 0;
                     0       0 0]

    rv_eci = [R O; R_dot R]*rv_ecef

    return rv_eci
end

"""
    ecef2eci(rv_ecef, θ; ω=7.292115373194e-5, ang_unit::Symbol=:rad)

Convert state vector `rv = [r; v]` {km; km/s} from Earth-Centered Inertial (ECI) frame to
Earth-Centered Earth-Fixed (ECEF) frame given Greenwich Mean Sidereal Time (GMST) `θ` {rad}.
"""
function eci2ecef(rv_eci, θ; ω=7.292115373194e-5, ang_unit::Symbol=:rad)
    if ang_unit == :deg
        θ = deg2rad(θ)
    elseif ang_unit == :rad
    else
        error("ang_unit should be :rad or :deg")
    end

    R = rotz(-θ)
    O = zeros(3,3)
    R_dot = ω*[sin(-θ) cos(-θ) 0;
              -cos(-θ) sin(-θ) 0;
                     0       0 0]

    rv_ecef = [R O; R_dot R]*rv_eci

    return rv_ecef
end

"""
    eci2sci(rv_eci, rv_sun_eci; ε=23.43929)
"""
function eci2sci(rv_eci, rv_sun_eci; ε=23.439292)
    R = rotx(-ε)
    O = zeros(3,3)
    rv_sci = [R O; O R]*(rv_eci-rv_sun_eci)
    return rv_sci
end

"""
    sci2eci(rv_sci, rv_sun_eci; ε=23.439292)
"""
function sci2eci(rv_sci, rv_sun_eci; ε=23.439292)
    R = rotx(ε)
    O = zeros(3,3)
    rv_eci = [R O; O R]*rv_sci + rv_sun_eci
    return rv_eci
end

"""
    enu2ecef(rv_enu, ϕ, λ, h=0; geodetic=true, ang_unit=:deg)
Converts state from ENU frame to inertial frame
"""
function enu2ecef(rv_enu, ϕ, λ, h=0; geodetic=true, ang_unit::Symbol=:deg, e_earth=0.0818, r_earth=6.378136e3)
    d = h + r_earth # [km] radius of the earth plus altitude
    if geodetic # geodetic latitude
        if ang_unit == :rad
            N = r_earth/sqrt(1-(e_earth^2)*sin(λ)^2)
        elseif ang_unit == :deg
            N = r_earth/sqrt(1-(e_earth^2)*sind(λ)^2)
        else
            error("ang_unit should be :deg or :rad")
        end
        r_radar_ecef = latlon2cart(ϕ, λ, 1, ang_unit=ang_unit)
        geodetic_factors = [N + h; N + h; N*(1-e_earth^2) + h]
        r_radar_ecef = r_radar_ecef.*geodetic_factors
    else # geocentric latitude and longitude
        r_radar_ecef = latlon2cart(ϕ, λ, d, ang_unit=ang_unit)
    end

    R = rotlatlon(ϕ, λ, ang_unit=ang_unit)

    r_ecef = R*rv_enu[1:3] + r_radar_ecef
    v_ecef = R*rv_enu[4:6]
    rv_ecef = [r_ecef; v_ecef]
    return rv_ecef
end

"""
    ecef2enu(rv_ecef, ϕ, λ, h=0; geodetic=true, ang_unit::Symbol=:deg, e_earth=0.0818, r_earth=6.378136e3)
Converts state from ENU frame to inertial frame
"""
function ecef2enu(rv_ecef, ϕ, λ, h=0; geodetic=true, ang_unit::Symbol=:deg, e_earth=0.0818, r_earth=6.378136e3)
    d = h + r_earth # [km] radius of the earth plus altitude

    R = rotlatlon(ϕ, λ, ang_unit=ang_unit)

    if geodetic # geodetic latitude
        if ang_unit == :rad
            N = r_earth/sqrt(1-(e_earth^2)*sin(λ)^2)
        elseif ang_unit == :deg
            N = r_earth/sqrt(1-(e_earth^2)*sind(λ)^2)
        else
            error("ang_unit should be :deg or :rad")
        end
        r_radar_ecef = latlon2cart(ϕ, λ, 1, ang_unit=ang_unit)
        geodetic_factors = [N + h; N + h; N*(1-e_earth^2) + h]
        r_radar_ecef = r_radar_ecef.*geodetic_factors
    else # geocentric latitude and longitude
        r_radar_ecef = latlon2cart(ϕ, λ, d, ang_unit=ang_unit)
    end
    r_enu = R'*(rv_ecef[1:3] - r_radar_ecef)
    v_enu = R'*rv_ecef[4:6]
    rv_enu = [r_enu; v_enu]
    return rv_enu
end

"""
    latlon2cart(ϕ, λ, d; ang_unit::Symbol=:deg)
"""
function latlon2cart(ϕ, λ, d; ang_unit::Symbol=:deg)
    if ang_unit == :rad
        r = d*[cos(ϕ)*cos(λ), cos(ϕ)*sin(λ), sin(ϕ)]
    elseif ang_unit == :deg
        r = d*[cosd(ϕ)*cosd(λ), cosd(ϕ)*sind(λ), sind(ϕ)]
    else
        error("ang_unit should be :deg or :rad")
    end
    return r
end

"""
    cart2latlon(r::Array)
"""
function cart2latlon(r::Array, ang_unit::Symbol=:deg)
    x,y,z = r
    xy = sqrt(x^2 + y^2)
    d = norm(r)
    if ang_unit == :rad
        ϕ = asin(z/d)
        λ = atan(y,x)
    elseif ang_unit == :deg
        ϕ = asind(z/d)
        λ = atand(y,x)
    else
        error("ang_unit should be :deg or :rad")
    end
    return ϕ, λ, d
end
# y   r
# |  /
# | /  lon
# |/------- x

function eci2CR3BP()
end

function CR3BP2eci()
end

function CR3BP2cart()
end

function sci2CR3BP()
end

function CR3BP2sci()
end

"""
    dimensionalize(rv, sys::System)

Add dimensions to normalized state vector `rv`
"""
function dimensionalize(rv, sys::System)
    return vcat(rv[1:3]*sys.RUNIT, rv[4:6]*sys.VUNIT)
end

"""
    dimensionalize!(rv, sys::System)

In-place version of dimensionalize(rv, sys::System)
"""
function dimensionalize!(rv, sys::System)
    rv[1:6] = dimensionalize(rv,sys)
    return nothing
end

"""
    nondimensionalize(rv, sys::System)

Normalize state vector `rv`
"""
function nondimensionalize(rv, sys::System)
    return vcat(rv[1:3]/sys.RUNIT, rv[4:6]/sys.VUNIT)
end

"""
    nondimensionalize!(rv, sys::System)

In-place version of nondimensionalize(rv, sys::System)
"""
function nondimensionalize!(rv, sys::System)
    rv[1:6] = nondimensionalize(rv,sys)
    return nothing
end


text:
using GeoJSONTables
using JSON3
import GeoInterface
using Tables
using Test

# copied from the GeoJSON.jl test suite
include("geojson_samples.jl")
featurecollections = [g, multipolygon, realmultipolygon, polyline, point, pointnull,
    poly, polyhole, collection, osm_buildings]

@testset "GeoJSONTables.jl" begin
    # only FeatureCollection supported for now
    @testset "Not FeatureCollections" begin
        @test_throws ArgumentError GeoJSONTables.read(a)
        @test_throws ArgumentError GeoJSONTables.read(b)
        @test_throws ArgumentError GeoJSONTables.read(c)
        @test_throws ArgumentError GeoJSONTables.read(d)
        @test_throws ArgumentError GeoJSONTables.read(e)
        @test_throws ArgumentError GeoJSONTables.read(f)
        @test_throws ArgumentError GeoJSONTables.read(h)
    end

    @testset "Read not crash" begin
        for featurecollection in featurecollections
            GeoJSONTables.read(featurecollection)
        end
    end

    @testset "FeatureCollection of one MultiPolygon" begin
        t = GeoJSONTables.read(g)
        @test Tables.istable(t)
        @test Tables.rows(t) === t
        @test Tables.columns(t) isa Tables.CopiedColumns
        @test t isa GeoJSONTables.FeatureCollection{<:JSON3.Array{JSON3.Object}}
        @test Base.propertynames(t) == (:json,)  # override this?
        @test Tables.rowtable(t) isa Vector{<:NamedTuple}
        @test Tables.columntable(t) isa NamedTuple

        f1, _ = iterate(t)
        @test f1 isa GeoJSONTables.Feature{<:JSON3.Object}
        @test all(Base.propertynames(f1) .== [:cartodb_id, :addr1, :addr2, :park])
        @test f1 == t[1]
        @test GeoJSONTables.geometry(f1) isa JSON3.Object
        @test GeoJSONTables.geometry(f1).type === "MultiPolygon"
        @test GeoJSONTables.geometry(f1).coordinates isa JSON3.Array
        @test length(GeoJSONTables.geometry(f1).coordinates[1][1]) == 4
        @test GeoJSONTables.geometry(f1).coordinates[1][1][1] == [-117.913883,33.96657]
        @test GeoJSONTables.geometry(f1).coordinates[1][1][2] == [-117.907767,33.967747]
        @test GeoJSONTables.geometry(f1).coordinates[1][1][3] == [-117.912919,33.96445]
        @test GeoJSONTables.geometry(f1).coordinates[1][1][4] == [-117.913883,33.96657]

        @testset "GeoInterface" begin
            @test GeoInterface.geotype(t) === :FeatureCollection
            @test GeoInterface.geotype(f1) === :Feature
            gi_mp = GeoInterface.geometry(f1)
            @test gi_mp isa GeoInterface.MultiPolygon
            @test GeoInterface.geotype(gi_mp) === :MultiPolygon
            properties = GeoInterface.properties(f1)
            @test properties isa Dict{String, Any}
            @test properties["addr2"] === "Rowland Heights"
            @test_throws MethodError GeoInterface.bbox(t)
            @test GeoInterface.bbox(f1) === nothing
            coordinates = GeoInterface.coordinates(f1)
            @test coordinates == Vector{Vector{Vector{Float64}}}[[[
                [-117.913883, 33.96657],
                [-117.907767, 33.967747],
                [-117.912919, 33.96445],
                [-117.913883, 33.96657],
            ]]]
            gi_f = GeoInterface.Feature(f1)
            @test gi_f isa GeoInterface.Feature
            gi_fc = GeoInterface.FeatureCollection(t)
            @test gi_fc isa GeoInterface.FeatureCollection
        end
    end
end


text:
abstract LogReg <: RegERM

methodname(::LogReg) = "Logistic Regression"
regularizer(LogReg::LogReg, w::Vector) = L2reg(w, LogReg.params.λ)

## binomial

immutable BinomialLogReg <: LogReg
    X::Matrix               # n x m matrix of n m-dimensional training examples
    y::Vector               # 1 x n vector with training classes
    n::Int                  # number of training examples
    m::Int                  # number of features
    kernel::Symbol          # kernel function
    regression_type::Symbol # ordinal, binomial, multinomial
    params::Hyperparameters # hyperparameters (e.g. λ)
end
function BinomialLogReg(X::Matrix, y::Vector; kernel::Symbol=:linear, λ::Float64=0.1)
    check_arguments(X, y, :binomial)
    BinomialLogReg(X, y, size(X)..., kernel, :binomial, RegularizationParameters(λ))
end

loss(::BinomialLogReg) = LogisticLoss()

## multinomial

immutable MultinomialLogReg <: LogReg
    X::Matrix               # n x m matrix of n m-dimensional training examples
    y::Vector               # 1 x n vector with training classes
    n::Int                  # number of training examples
    m::Int                  # number of features
    kernel::Symbol          # kernel function
    regression_type::Symbol # ordinal, binomial, multinomial
    params::Hyperparameters # hyperparameters (e.g. λ)
end
function MultinomialLogReg(X::Matrix, y::Vector; kernel::Symbol=:linear, λ::Float64=0.1)
    check_arguments(X, y, :multinomial)
    MultinomialLogReg(X, y, size(X)..., kernel, :multinomial, RegularizationParameters(λ))
end

loss(::MultinomialLogReg) = MultinomialLogisticLoss()