text:
using FoundationDB
using Test

function timewaittask(watchtask, ref_twatch)
    f = () -> begin
        t1 = time()
        fetch(watchtask)
        ref_twatch[] = time() - t1
    end
    @async f()
end

function arraycmp(a, b, cmp)
    cmp(a, b)
end

@testset "start network" begin
    @test !is_client_running()
    @test start_client() === nothing
    @test is_client_running()
end # testset "start network"

try
    @testset "basic julia apis" begin
        chk_closed = Vector{Any}()

        open(FDBCluster()) do cluster
            push!(chk_closed, cluster)
            @test cluster.ptr !== C_NULL

            open(FDBDatabase(cluster)) do db
                push!(chk_closed, db)
                @test db.name == "DB"
                @test db.ptr !== C_NULL

                open(FDBTransaction(db)) do tran
                    @test clearkeyrange(tran, UInt8[0], UInt8[2]) == nothing
                end

                key = UInt8[0,1,2]
                val = UInt8[9, 9, 9]

                open(FDBTransaction(db)) do tran
                    push!(chk_closed, tran)
                    @test tran.ptr !== C_NULL

                    @test clearkey(tran, key) == nothing
                    @test getval(tran, key) == nothing
                    @test setval(tran, key, val) == nothing
                    @test getval(tran, key) == val
                end

                open(FDBTransaction(db)) do tran
                    push!(chk_closed, tran)
                    @test clearkey(tran, key) == nothing
                    @test getval(tran, key) == nothing
                end

                open(FDBTransaction(db)) do tran
                    push!(chk_closed, tran)
                    @test getval(tran, key) == nothing
                end
            end
        end

        for item in chk_closed
            @test !isopen(item)
        end
    end # testset "basic julia apis"

    @testset "auto commit" begin
        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                key = UInt8[0,1,2]
                val = UInt8[9, 9, 9]
                open(FDBTransaction(db)) do tran
                    @test clearkey(tran, key) == nothing
                    @test getval(tran, key) == nothing
                    @test setval(tran, key, val) == nothing
                    @test getval(tran, key) == val
                    @test commit(tran)
                    @test_throws FDBError commit(tran)
                end

                open(FDBTransaction(db)) do tran
                    @test clearkey(tran, key) == nothing
                    @test getval(tran, key) == nothing
                    @test commit(tran)
                end

                open(FDBTransaction(db)) do tran
                    @test getval(tran, key) == nothing
                    @test commit(tran)
                end
            end
        end
    end # testset "auto commit"

    @testset "cancel" begin
        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                key = UInt8[0,1,2]
                val = UInt8[9, 9, 9]
                open(FDBTransaction(db)) do tran
                    @test setval(tran, key, val) == nothing
                    @test cancel(tran) == nothing
                end
                open(FDBTransaction(db)) do tran
                    @test getval(tran, key) == nothing
                end
                open(FDBTransaction(db)) do tran
                    @test setval(tran, key, val) == nothing
                    @test commit(tran)
                    @test cancel(tran) == nothing
                end
                open(FDBTransaction(db)) do tran
                    @test getval(tran, key) == val
                    @test clearkey(tran, key) == nothing
                end
            end
        end
    end # testset "cancel"

    @testset "parallel updates" begin
        sumval = 0
        function do_updates()
            open(FDBCluster()) do cluster
                open(FDBDatabase(cluster)) do db
                    key = UInt8[0,1,2]
                    valarr = Int[1]
                    val = FoundationDB.unsafe_wrap(Array, convert(Ptr{UInt8}, pointer(valarr)), (sizeof(Int),), own=false)
                    for valint in 1:100
                        valarr[1] = valint
                        open(FDBTransaction(db)) do tran
                            @test setval(tran, key, val) == nothing
                        end
                        open(FDBTransaction(db)) do tran
                            @test clearkey(tran, key) == nothing
                        end
                        sleep(rand()/100)
                        open(FDBTransaction(db)) do tran
                            valnow = getval(tran, key)
                            if valnow != nothing
                                sumval += reinterpret(Int, valnow)[1]
                            end
                        end
                    end
                end
            end
        end

        @sync begin
            @async do_updates()
            @async do_updates()
        end

        println("sumval in parallel updates = ", sumval)
        @test 0 <= sumval <= (5050 * 1.5) # series sum of 1:100 = 5050
    end # testset "parallel updates"

    @testset "large key value" begin
        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                key = ones(UInt8, 10000)
                val = ones(UInt8, 100000)
                open(FDBTransaction(db)) do tran
                    @test clearkey(tran, key) == nothing
                    @test getval(tran, key) == nothing
                    @test setval(tran, key, val) == nothing
                    @test getval(tran, key) == val
                    @test commit(tran)
                    @test reset(tran) == nothing
                    @test commit(tran)              # test that commit is allowed after a reset
                end

                open(FDBTransaction(db)) do tran
                    @test clearkey(tran, key) == nothing
                    @test getval(tran, key) == nothing
                end

                open(FDBTransaction(db)) do tran
                    @test getval(tran, key) == nothing
                    @test reset(tran) == nothing
                end
            end
        end
    end # testset "large key value"

    @testset "watch" begin
        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                key = UInt8[0,1,2]
                val1 = UInt8[0, 0, 0]
                val2 = UInt8[0, 0, 0]
                open(FDBTransaction(db)) do tran
                    ref_twatch = Ref(0.0)
                    watchtask = watchkey(tran, key)
                    timetask = timewaittask(watchtask, ref_twatch)
                    sleep(0.5)
                    @test setval(tran, key, val2) == nothing
                    sleep(0.5)
                    @test clearkey(tran, key) == nothing
                    fetch(timetask)
                    @test ref_twatch[] > 0.4

                    watchhandle = FDBFuture()
                    watchtask = watchkey(tran, key; handle=watchhandle)
                    sleep(0.2)
                    @test istaskstarted(watchtask)
                    @test !istaskdone(watchtask)
                    cancel(watchhandle)
                    sleep(0.2)
                    @test istaskdone(watchtask)
                    try
                        fetch(watchtask)
                        error("watchtask should have failed!")
                    catch ex
                        @test isa(ex, FDBError)
                        @test ex.code == 1101 # Asynchronous operation cancelled
                    end
                end
            end
        end

        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                key = UInt8[0,1,2]
                val1 = UInt8[0, 0, 0]
                val2 = UInt8[0, 0, 0]

                # set an initial value
                open(FDBTransaction(db)) do tran
                    @test setval(tran, key, val1) == nothing
                end

                ref_twatch = Ref(0.0)
                # start a watch
                watchtask = open(FDBTransaction(db)) do tran
                    watchkey(tran, key)
                end

                timetask = timewaittask(watchtask, ref_twatch)

                open(FDBTransaction(db)) do tran
                    sleep(0.5)
                    @test setval(tran, key, val2) == nothing
                    sleep(0.5)
                    @test clearkey(tran, key) == nothing
                end

                fetch(timetask)
                @test ref_twatch[] > 0.4
            end
        end
    end # testset "watch"

    @testset "get key" begin
        keys = [UInt8[0,1,x] for x in 1:20]
        val = UInt8[0]

        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                # setup all keys
                open(FDBTransaction(db)) do tran
                    for key in keys
                        @test setval(tran, key, val) == nothing
                    end
                end

                # get
                open(FDBTransaction(db)) do tran
                    emptykey = UInt8[]
                    key = UInt8[0,0,0]
                    @test getkey(tran, keysel(FDBKeySel.last_less_or_equal, key)) == emptykey
                    @test getkey(tran, keysel(FDBKeySel.last_less_than, key)) == emptykey
                    @test getkey(tran, keysel(FDBKeySel.first_greater_than, key)) == UInt8[0,1,1]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_or_equal, key)) == UInt8[0,1,1]

                    key = UInt8[0,1,1]
                    @test getkey(tran, keysel(FDBKeySel.last_less_or_equal, key)) == UInt8[0,1,1]
                    @test getkey(tran, keysel(FDBKeySel.last_less_than, key)) == emptykey
                    @test getkey(tran, keysel(FDBKeySel.first_greater_than, key)) == UInt8[0,1,2]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_or_equal, key)) == UInt8[0,1,1]

                    key = UInt8[0,1,10]
                    @test getkey(tran, keysel(FDBKeySel.last_less_or_equal, key)) == UInt8[0,1,10]
                    @test getkey(tran, keysel(FDBKeySel.last_less_than, key)) == UInt8[0,1,9]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_than, key)) == UInt8[0,1,11]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_or_equal, key)) == UInt8[0,1,10]

                    key = UInt8[0,1,20]
                    @test getkey(tran, keysel(FDBKeySel.last_less_or_equal, key)) == UInt8[0,1,20]
                    @test getkey(tran, keysel(FDBKeySel.last_less_than, key)) == UInt8[0,1,19]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_than, key)) != UInt8[0,1,20]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_or_equal, key)) != UInt8[0,1,19]

                    key = UInt8[0,2,0]
                    @test getkey(tran, keysel(FDBKeySel.last_less_or_equal, key)) == UInt8[0,1,20]
                    @test getkey(tran, keysel(FDBKeySel.last_less_than, key)) == UInt8[0,1,20]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_than, key)) != UInt8[0,2,0]
                    @test getkey(tran, keysel(FDBKeySel.first_greater_or_equal, key)) != UInt8[0,2,0]
                end

                # clear all keys
                open(FDBTransaction(db)) do tran
                    for key in keys
                        @test clearkey(tran, key) == nothing
                    end
                end
            end
        end
    end # testset "get key"

    @testset "get range" begin
        keys = [UInt8[0,1,x] for x in 1:20]
        val = UInt8[1]

        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                # setup all keys
                open(FDBTransaction(db)) do tran
                    for key in keys
                        @test setval(tran, key, val) == nothing
                    end
                end

                # get
                open(FDBTransaction(db)) do tran
                    kvs, more = getrange(tran, keysel(FDBKeySel.first_greater_or_equal, keys[1]), keysel(FDBKeySel.last_less_or_equal, keys[10]))
                    @test length(kvs) == 9
                    @test !more
                    for kv in kvs
                        @test kv[1] in keys
                        @test kv[2] == val
                    end

                    kvs, more = getrange(tran, keysel(FDBKeySel.first_greater_or_equal, keys[1]), keysel(FDBKeySel.last_less_or_equal, keys[10]); limit=1)
                    @test length(kvs) == 1
                    @test more
                end

                # clear all keys
                open(FDBTransaction(db)) do tran
                    @test clearkeyrange(tran, keys[1], UInt8[0,1,21]) == nothing
                    @test getval(tran, keys[10]) == nothing
                end
            end
        end
    end # testset "get range"

    @testset "atomic" begin
        key    = UInt8[0,1,2]
        keymax = UInt8[0,1,3]
        keymin = UInt8[0,1,4]
        keyand = UInt8[0,1,5]
        keyor  = UInt8[0,1,6]
        keyxor = UInt8[0,1,7]

        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                @testset "atomic add" begin
                    # clear key to start with
                    open(FDBTransaction(db)) do tran
                        @test clearkey(tran, key) == nothing
                    end
                    @sync begin
                        for idx in 1:10
                            @async open(FDBTransaction(db)) do tran
                                atomic_add(tran, key, 1)
                            end
                        end
                        for idx in 1:5
                            @async open(FDBTransaction(db)) do tran
                                atomic_add(tran, key, -1)
                            end
                        end
                    end
                    open(FDBTransaction(db)) do tran
                        @test atomic_integer(Int, getval(tran, key)) == 5
                        @test clearkey(tran, key) == nothing
                    end
                end # testset "atomic add"

                @testset "atomic integer min max" begin
                    # clear key to start with
                    open(FDBTransaction(db)) do tran
                        @test clearkey(tran, keymax) == nothing
                        @test clearkey(tran, keymin) == nothing
                    end
                    vals = rand(Cuint, 10)
                    @sync begin
                        for idx in 1:10
                            @async open(FDBTransaction(db)) do tran
                                atomic_max(tran, keymax, vals[idx])
                                atomic_min(tran, keymin, vals[idx])
                            end
                        end
                    end
                    open(FDBTransaction(db)) do tran
                        @test atomic_integer(Cuint, getval(tran, keymin)) == minimum(vals)
                        @test clearkey(tran, keymin) == nothing

                        @test atomic_integer(Cuint, getval(tran, keymax)) == maximum(vals)
                        @test clearkey(tran, keymax) == nothing
                    end
                end # testset "atomic integer min max"

                @testset "atomic bytes min max and or xor" begin
                    # clear key to start with
                    open(FDBTransaction(db)) do tran
                        @test setval(tran, keymax, zeros(UInt8, 10)) == nothing
                        @test setval(tran, keymin, ones(UInt8, 10) * typemax(UInt8)) == nothing
                    end
                    open(FDBTransaction(db)) do tran
                        @test getval(tran, keymax) == zeros(UInt8, 10)
                        @test getval(tran, keymin) == ones(UInt8, 10) * typemax(UInt8)
                    end
                    vals = [rand(UInt8,10) for x in 1:10]
                    maxval = zeros(UInt8, 10)
                    minval = ones(UInt8, 10) * typemax(UInt8)
                    andval = vals[1]
                    orval = vals[1]
                    xorval = vals[1]
                    open(FDBTransaction(db)) do tran
                        @test setval(tran, keyand, andval) == nothing
                        @test setval(tran, keyor, orval) == nothing
                        @test setval(tran, keyxor, xorval) == nothing
                    end
                    for val in vals
                        for idx in 1:length(val)
                            arraycmp(val, maxval, >) && (maxval = val)
                            arraycmp(val, minval, <) && (minval = val)
                            if idx > 1
                                andval = andval .& val
                                orval = orval .| val
                                xorval = xorval .⊻ val
                            end
                        end
                    end
                    @sync begin
                        for idx in 1:10
                            @async open(FDBTransaction(db)) do tran
                                atomic_max(tran, keymax, vals[idx])
                                atomic_min(tran, keymin, vals[idx])
                                atomic_and(tran, keyand, vals[idx])
                                atomic_or(tran, keyor, vals[idx])
                                atomic_xor(tran, keyxor, vals[idx])
                            end
                        end
                    end
                    open(FDBTransaction(db)) do tran
                        for (k,v) in zip([keymin, keymax, keyand, keyor, keyxor], [minval, maxval, andval, orval, xorval])
                            @test getval(tran, k) == v
                            @test clearkey(tran, k) == nothing
                        end
                    end
                end # testset "atomic bytes min max and or xor"
            end
        end
    end # testset "atomic"

    @testset "conflicts" begin
        function test_conflict(db, my_trigger::Channel{Nothing}, other_trigger::Channel{Nothing})
            open(FDBTransaction(db)) do tran
                # add conflict range
                conflict(tran, UInt8[0,1,2], UInt8[0,1,3], FDBConflictRangeType.READ)
                conflict(tran, UInt8[0,1,2], UInt8[0,1,3], FDBConflictRangeType.WRITE)

                # read key, ensure it is there
                if getval(tran, UInt8[0,1,2]) == nothing
                    error("not found")
                end

                # send trigger to other transaction
                put!(other_trigger, nothing)

                # wait for other transaction to also read the key
                wait(my_trigger)

                # delete the key and try to commit
                clearkey(tran, UInt8[0,1,2])
            end
        end

        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                open(FDBTransaction(db)) do tran
                    setval(tran, UInt8[0,1,2], UInt8[0,1,2])
                end

                c1 = Channel{Nothing}(1)
                c2 = Channel{Nothing}(1)

                try
                    @sync begin
                        @async test_conflict(db, c1, c2)
                        @async test_conflict(db, c2, c1)
                    end
                    error("must throw not found")
                catch ex
                    @test isa(ex, CompositeException)
                    @test length(ex.exceptions) == 1
                    @test isa(ex.exceptions[1], CapturedException)
                    @test ex.exceptions[1].ex.msg == "not found"
                end

                open(FDBTransaction(db)) do tran
                    clearkey(tran, UInt8[0,1,2])
                end
            end
        end
    end # testset "conflicts"

    @testset "versionstamp" begin
        open(FDBCluster()) do cluster
            open(FDBDatabase(cluster)) do db
                open(FDBTransaction(db; trackversionstamp=true)) do tran
                    @test tran.versionstamp === nothing
                    atomic_setval(tran, zeros(UInt8, 8), ones(UInt8, 16), FDBMutationType.SET_VERSIONSTAMPED_VALUE)
                    commit(tran)
                    @test tran.versionstamp !== nothing
                    v1 = tran.versionstamp

                    reset(tran)
                    clearkey(tran, zeros(UInt8, 8))
                    commit(tran)
                    @test tran.versionstamp !== nothing
                    @test tran.versionstamp == v1

                    reset(tran)
                    key = UInt8[0,1,2]
                    prep_atomic_key!(key, 2)
                    atomic_setval(tran, key, zeros(UInt8, 5), FDBMutationType.SET_VERSIONSTAMPED_KEY)
                    commit(tran)
                    @test tran.versionstamp !== nothing
                    @test tran.versionstamp != v1
                    v2 = tran.versionstamp

                    reset(tran)
                    key = zeros(UInt8, 13)
                    key[2:11] = v2
                    key[12:13] = UInt8[1,2]
                    @test getval(tran, key) == zeros(UInt8, 5)
                    clearkey(tran, key)
                end
            end
        end
    end

    @testset "stop network" begin
        @test is_client_running()
        @test stop_client() === nothing
        @test !is_client_running()
        @test_throws Exception start_client()
    end
finally
    stop_client()
end


text:
alpha = 0.85
list = "test/fixtures/dangling.adj"
order = 6

M, order, max_ordinal = fastadj(Float64, Int64, alpha, list, 8)

expected = sparse([
 0.0  0.85     0.0      0.0      0.0      0.0
 0.0  0.0      0.425    0.425    0.0      0.0
 0.0  0.28333  0.0      0.28333  0.28333  0.0
 0.0  0.28333  0.28333  0.0      0.28333  0.0
 0.0  0.0      0.28333  0.28333  0.0      0.28333
 0.0  0.0      0.0      0.0      0.0      0.0
])

difference = (max(dense(expected) - dense(M)))

@assert isequal(6, order)
@assert isequal(6, max_ordinal)
@assert isequal(0, difference)

ext: jl


text:
using DocPkg
using Test

@testset "DocPkg.jl" begin
    @test 2*1 == 2
end


text:
export PhylogeneticTree, phylogenetic_tree, calculate_tree, newick_string, leaves, isleaf,
    set_property!, get_property, annotate!

struct PhylogeneticTree
    graph::MetaDiGraph
end

graph(tree::PhylogeneticTree) = getfield(tree, :graph)

function phylogenetic_tree(newick_string::String)
    graph = MetaDiGraph()
    add_vertex!(graph)
    set_prop!(graph, 1, :name, "root")

    components = newick_components(newick_string)
    for component in components
        add_to_graph!(graph, 1, component)
    end

    return PhylogeneticTree(graph)
end

function phylogenetic_tree(data::AbstractHLAData; model::String = "LG+G+I",
        verbose::Bool = false)
    !ismissing(data.tree) && return data.tree

    temp_prefix = tempname()
    temp_fasta = temp_prefix * ".fasta"

    write_numbered_fasta(data, temp_fasta)

    if verbose
        Base.run(`$(raxml_ng()) --msa $temp_fasta --model $model --threads auto
            --workers auto`)
    else
        @suppress Base.run(`$(raxml_ng()) --msa $temp_fasta --model $model --threads auto
        --workers auto`)
    end
    tree = phylogenetic_tree(readline(temp_fasta * ".raxml.bestTree"))

    raxml_files = ["bestModel", "bestTree", "log", "mlTrees", "rba", "startTree"]
    rm(temp_fasta)
    rm.(temp_fasta .* ".raxml." .* raxml_files)

    return tree
end

function write_numbered_fasta(data::AbstractHLAData, filepath::String)
    records = Escape.records(data)

    write_numbered_fasta(records, filepath)
end

function write_numbered_fasta(records::Vector{FASTX.FASTA.Record}, filepath::String)
    open(FASTA.Writer, filepath) do writer
        for (i, record) in enumerate(records)
            identifier = string(i)
            sequence = FASTA.sequence(record)
            new_record = FASTA.Record(identifier, sequence)

            write(writer, new_record)
        end
    end
end

function add_to_graph!(graph::AbstractMetaGraph, vertex::Int, s::String)
    name = extract_name(s)
    length = extract_length(s)

    add_vertex!(graph)
    set_prop!(graph, nv(graph), :name, name)

    add_edge!(graph, vertex, nv(graph))
    set_prop!(graph, Edge(vertex, nv(graph)), :branch_length, length)

    current_vertex = nv(graph)
    if occursin(')', s)
        components = newick_components(s)
        for c in components
            add_to_graph!(graph, current_vertex, c)
        end
    end
end

function extract_name(s::String)
    substring = occursin(')', s) ? s[findlast(isequal(')'), s) + 1 : end] : s
    substring != "" || return missing

    name = String(split(substring, ":")[1])
    return name != "" ? name : missing
end

function extract_length(s::String)
    substring = occursin(')', s) ? s[findlast(isequal(')'), s) + 1 : end] : s
    occursin(':', substring) || return missing

    length = parse(Float64, split(substring, ":")[end])
end

function newick_components(s::String)
    components = Vector{String}()
    component = ""

    start = occursin('(', s) ? findfirst(isequal('('), s) + 1 : 1
    stop = occursin(')', s) ? findlast(isequal(')'), s) - 1 : length(s)

    for c in s[start:stop]
        if c == ',' && count(x -> x == '(', component) == count(x -> x == ')', component)
            push!(components, strip(component))
            component = ""
        else
            component *= c
        end
    end

    push!(components, strip(component))
    return components
end

function newick_string(tree::PhylogeneticTree)
    graph = Escape.graph(tree)
    n = Graphs.neighbors(graph, 1)
    length(n) == 0 && return "();"
    length(n) == 1 && return "(" * newick_string(graph, n[1]) * ");"

    newick = "(" * newick_string(graph, Graphs.neighbors(graph, 1)[1])
    for n in Graphs.neighbors(graph, 1)[2:end]
        newick *= ","
        newick *= newick_string(graph, n)
    end
    newick *= ");"

    return newick
end

function newick_string(graph::AbstractMetaGraph, v::Int)
    name = get_prop(graph, v, :name)
    name = ismissing(name) ? "" : name
    branch_length = get_prop(graph, inneighbors(graph, v)[1], v, :branch_length)
    branch_length = ismissing(branch_length) ? "" : string(branch_length)

    if length(Graphs.neighbors(graph, v)) == 0
        return string(name, branch_length != "" ? ":" : "", branch_length)
    else
        newick = "(" * newick_string(graph, Graphs.neighbors(graph, v)[1])
        for n in Graphs.neighbors(graph, v)[2:end]
            newick *= ","
            newick *= newick_string(graph, n)
        end
        newick *= ")" * name
        newick = !ismissing(length) ? newick * ":$branch_length" : newick
    end

    return newick
end

function isleaf(v::Int, tree::PhylogeneticTree)
    v ∈ leaves(tree) && return true
    return false
end

function leaves(tree::PhylogeneticTree)
    graph = Escape.graph(tree)
    leaves = map(x -> x[1], attracting_components(graph))

    return leaves
end

function get_property(tree::PhylogeneticTree, v::Int, property::Symbol)
    graph = Escape.graph(tree)
    has_prop(graph, v, property) || return missing

    return get_prop(graph, v, property)
end

function set_property!(tree::PhylogeneticTree, v::Int, property::Symbol, value)
    graph = Escape.graph(tree)

    set_prop!(graph, v, property, value)
end

function annotate!(tree::PhylogeneticTree, data::AbstractHLAData, replacement::Replacement)
    name(data) == protein(replacement) || error("Replacement does not match HLA data.")
    records = Escape.records(data)

    for (i, record) in enumerate(records)
        sequence = FASTA.sequence(record)
        symbol = Char(sequence[position(replacement)])

        v = filter(x -> get_property(tree, x, :name) == string(i), leaves(tree))[1]

        if !negated(replacement)
            if symbol == Escape.replacement(replacement)
                set_property!(tree, v, :state, "1")
            elseif symbol in ('-', 'X')
                set_property!(tree, v, :state, missing)
            else
                set_property!(tree, v, :state, "0")
            end
        elseif negated(replacement)
            if symbol == Escape.replacement(replacement)
                set_property!(tree, v, :state, "0")
            elseif symbol in ('-', 'X')
                set_property!(tree, v, :state, missing)
            else
                set_property!(tree, v, :state, "1")
            end
        end
    end
end

function matching(tree::PhylogeneticTree, data::AbstractHLAData)
    records = Escape.records(data)
    matching = Escape.matching(tree, records)

    return matching
end

function matching(tree::PhylogeneticTree, records::Vector{FASTX.FASTA.Record})
    n_leaves = length(leaves(tree))

    msg = string(
        "Number of leaves ($n_leaves) does not match number of sequences ",
        "($(length(records)))"
    )
    n_leaves == length(records) || return DimensionMismatch(msg)

    msg = """:name property of leaves must be "1", "2", ..., "$n_leaves"."""
    leaf_names = Set(map(x -> get_property(tree, x, :name), leaves(tree)))
    leaf_names == Set(string.(1:n_leaves)) || return ErrorException(msg)

    return true
end

function states(tree::PhylogeneticTree)
    states = Vector{String}(undef, length(leaves(tree)))

    for leaf in leaves(tree)
        state = get_property(tree, leaf, :state)
        name = get_property(tree, leaf, :name)

        states[parse(Int, name)] = ismissing(state) ? "-" : state
    end

    return states
end


text:
const CACHE = GlobalCache()

"""
Get the raw docstrings associated with all documented objects in a module.
"""
getraw(m::Module) = getraw(CACHE, m)

"""
Get the raw docstring associated with a documented object ``obj`` in module ``m``.
"""
getraw(m::Module, obj) = getraw(CACHE, m, obj)

"""
Get the parsed docstrings associated with all documented objects in a module.
"""
getparsed(m::Module) = getparsed(CACHE, m)

"""
Get the parsed form of a docstring associated with an object ``obj``.

Automatically parses the raw docstring on demand when first called.
Subsequent calls will return the parsed docstring that has been cached.
"""
getparsed(m::Module, obj) = getparsed(CACHE, m, obj)

"""
Get the metadata dictionaries for all documented objects in a module.
"""
getmeta(m::Module) = getmeta(CACHE, m)

"""
Get the ``Dict{Symbol, Any}`` containing an object's metadata.
"""
getmeta(m::Module, obj) = getmeta(CACHE, m, obj)

"""
Get the ``PackageData`` object associated with a module ``m``.
"""
getpackage(m::Module) = getpackage(CACHE, m)

"""
Get the ``ModuleData`` object associated with a module ``m``.
"""
getmodule(m::Module) = getmodule(CACHE, m)

"""
Empty the documentation cache of all data.
"""
clear!() = empty!(CACHE)

"""
Return all documented objects found in a module ``m``.
"""
objects(m::Module) = objects(CACHE, m)

"""
Turn on documenting of ``Base`` and it's submodules. Off by default.
"""
togglebase() = togglebase(CACHE)

"""
Find the metadata for ``key`` associated with an object ``obj`` in module ``m``.

Returns a ``Nullable{T}`` object. ``isnull`` must be called to determine whether
an object was actually found or not.

When ``obj`` does not contain the field ``key`` then the module's metadata
and all it's parents are searched in turn. Finally the package's metadata is
searched for ``key``.
"""
findmeta(m::Module, obj, key::Symbol, T) = findmeta(CACHE, m, obj, key, T)

"""
Returns the set of all loaded modules.
"""
loadedmodules() = loadedmodules(CACHE)


text:
function apply_task_metric_chart(pn, vn, wn, task::Task{<:TaskMap{M,N}},
        C::Chart{I,N}) where {M,N,I}
    G = metric_chart(pn, task, C)
    vn'*G*wn
end

task_riemannian_norm_chart(pn, vn, task::Task{<:TaskMap{M,N}}, C::Chart{I,N}) where {M,N,I} =
    sqrt(apply_task_metric(pn, vn, vn, task, C))

function metric_chart_transition(p1, G1, C1::Chart{I,M}, C2::Chart{J,M}) where {M,I,J}
    p2 = chart_transition(p1, C1, C2)
    ∂p1_∂p2 = chart_transition_jacobian(p2, C2, C1)
    G2 = ∂p1_∂p2'*G1*∂p1_∂p2
end

function metric_from_home_chart(pni, task::Task{<:TaskMapT}, Ci)
    Ch = home_task_chart(task, Ci)
    pnh = chart_transition(pni, T(Ci), T(Ch))
    Gh = metric_chart(pnh, task, Ch)
    Gi = metric_chart_transition(pnh, Gh, T(Ch), T(Ci))
end

function metric_from_home_chart(pni, task::Task{<:BaseTaskMap}, Ci)
    Ch = home_task_chart(task, Ci)
    pnh = chart_transition(pni, Ci, Ch)
    Gh = metric_chart(pnh, task, Ch)
    Gi = metric_chart_transition(pnh, Gh, Ch, Ci)
end

function metric_from_emb(pn, task::Task{<:BaseTaskMap}, C)
    pne = chart_to_emb(pn, C)
    Ge = metric_emb(pne, task)
    ∂pne_∂pn = chart_to_emb_jacobian(pn, C)
    G = ∂pne_∂pn'*Ge*∂pne_∂pn
end

metric_coord_rep(task::Task{<:Union{BaseTaskMap,TaskMapT}}) = default_coord_rep(task)
metric_chart(pni, task::Task{<:Union{BaseTaskMap,TaskMapT}}, Ci::Chart) =
    metric_chart(metric_coord_rep(task), pni, task, Ci)
metric_chart(::EmbRep, pni, task, Ci) = metric_from_emb(pni, task, Ci)
metric_chart(::ChartRep, pni, task, Ci) = metric_from_home_chart(pni, task, Ci)
metric_emb(pne, task) = throw("Not defined!")

default_metric(pni, task::Task{<:TaskMap{M,ℝ{n}}}, CN1::Chart{1,ℝ{n}}) where {M,n} =
    SMatrix{n,n,eltype(pni)}(I)
default_metric(pni, task::Task{<:TaskMap{M,S1}}, CN1::Chart{<:AngleChart,S1}) where M =
    SMatrix{1,1,eltype(pni)}(I)
default_metric(pne, task::Task{<:TaskMap{M,𝕊{n}}}) where {M,n} =
    SMatrix{n+1,n+1,eltype(pne)}(I)

## Weight metrics
weight_metric_coord_rep(task::Task{<:Union{BaseTaskMap,TaskMapT}}) = default_coord_rep(task)

function weight_metric_from_home_chart(pni, wni, task::Task{<:BaseTaskMap}, Ci)
    Ch = home_task_chart(task, Ci)
    pnh, wnh = chart_transition_differential(pni, wni, Ci, Ch)
    Wh = weight_metric_chart(pnh, wnh, task, Ch)
    Wi = metric_chart_transition(pnh, Wh, Ch, Ci)
end

function weight_metric_from_emb(pn, wn, task::Task{<:BaseTaskMap}, C)
    pne, wne = chart_to_emb_differential(pn, wn, C)
    We = weight_metric_emb(pne, wne, task)
    ∂pne_∂pn = chart_to_emb_jacobian(pn, C)
    W = ∂pne_∂pn'*We*∂pne_∂pn
end

weight_metric_chart(pn, wn, task::Task{<:Union{BaseTaskMap,TaskMapT}}, C::Chart) =
    weight_metric_chart(weight_metric_coord_rep(task), pn, wn, task, C)
weight_metric_chart(::EmbRep, pn, wn, task, C) = weight_metric_from_emb(pn, wn, task, C)
weight_metric_chart(::ChartRep, pn, wn, task, C) =
    weight_metric_from_home_chart(pn, wn, task, C)
weight_metric_emb(pne, wne, task) = throw("Not defined!")

# function active_weight_position_from_home_chart(pn, task::Task{<:BaseTaskMap}, C)
#     Ch = home_task_chart(task, C)
#     pnh = chart_transition(pni, Ci, Ch)
#     active_weight_position_chart(pnh, task, Ch)
# end

# function active_weight_position_from_emb(pn, task::Task{<:BaseTaskMap}, C)
#     pne = chart_to_emb(pn, C)
#     active_weight_position_emb(pne, task)
# end

# active_weight_position_chart(pn, task::Task{<:Union{BaseTaskMap,TaskMapT}}, C::Chart) =
#     active_weight_position_chart(weight_metric_coord_rep(task), pn, C)
# active_weight_position_chart(::EmbRep, pn, C) = active_weight_position_from_emb(pn, wn, task, C)
# active_weight_position_chart(::ChartRep, pn, C) =
#     active_weight_position_from_home_chart(pn, task, C)
# active_weight_position_emb(pne, task) = throw("Not defined!")

active_weight_position_chart(pn, task::Task{<:Union{BaseTaskMap,TaskMapT}}, C::Chart) = true

default_weight_metric(pn, wn, task::Task{<:TaskMap{M,ℝ{n}}}, CN::Chart{1,ℝ{n}}) where {M,n} =
    SMatrix{n,n,eltype(pn)}(I)
default_weight_metric(pn, wn, task::Task{<:TaskMap{M,S1}}, CN::Chart{<:AngleChart,S1}) where M =
    SMatrix{1,1,eltype(pn)}(I)
default_weight_metric(pne, vne, task::Task{<:TaskMap{M,𝕊{n}}}) where {M,n} =
    SMatrix{n+1,n+1,eltype(pne)}(I)

## GDS metrics
function metric_from_home_chart(xni, vni, task, Ci::Chart{I,N}) where {N,I}
    Ch = home_task_chart(task, Ci)
    xnh, vnh = chart_transition_differential(xni, vni, Ci, Ch)
    Gh = metric_chart(xnh, vnh, task, Ch)

    # Since the block of the metric used for GDS is on the horizontal bundle, cannot perform a chart
    # transition without considering the full metric
    n, nt = dim(N), dim(T{N})
    Gh_full = SMatrix{nt,nt,eltype(Gh)}([Gh zeros(n,n); zeros(n,n) I])
    pnh = [xnh; vnh]
    Gi_full = metric_chart_transition(pnh, Gh, T(Ch), T(Ci))
    Gi = SMatrix{n,n,eltype(Gh)}(Gi_full[1:n,1:n])

    # If projection was onto vertical bundle:
    # Gi = metric_chart_transition(xnh, Gh, Ch, Ci)
end

function metric_from_emb(xn, vn, task, C::Chart{J,N}) where {N,J}
    xne, vne = chart_to_emb_differential(xn, vn, C)
    Ge = metric_emb(xne, vne, task)
    # Since the block of the metric used for GDS is on the horizontal bundle, cannot perform a chart
    # transition without considering the full metric
    n, ne, nte = dim(N), embdim(N), embdim(T{N})
    Ge_full = SMatrix{nte,nte,eltype(Ge)}([Ge zeros(ne,ne); zeros(ne,ne) I])
    pn = [xn; vn]
    ∂pne_∂pn = chart_to_emb_jacobian(pn, T(C))
    G_full = ∂pne_∂pn'*Ge_full*∂pne_∂pn
    G = SMatrix{n,n,eltype(Ge)}(G_full[1:n,1:n])

    # If projection was onto vertical bundle:
    # ∂xne_∂xn = chart_to_emb_jacobian(xn, C)
    # G = ∂xne_∂xn'*Ge*∂xne_∂xn
end

metric_coord_rep(task::TaskGDS) = default_coord_rep(task)
metric_chart(xni, vni, task::TaskGDS, Ci) = metric_chart(metric_coord_rep(task), xni, vni, task, Ci)
metric_chart(::EmbRep, xni, vni, task, Ci) = metric_from_emb(xni, vni, task, Ci)
metric_chart(::ChartRep, xni, vni, task, Ci) = metric_from_home_chart(xni, vni, task, Ci)
metric_emb(xne, vne, task) = throw("Not defined!")

default_metric(xn, vn, task::TaskGDS{<:TaskMapT{M,ℝ{n}}}, CN::Chart{1,ℝ{n}}) where {M,n} =
    SMatrix{n,n,eltype(xn)}(I)


text:
# Some utility functions for optimizing linear algebra operations that aren't specific
# to any particular rule definition

# F .* (X - X'), overwrites X if possible
function _mulsubtrans!!(X::AbstractMatrix{T}, F::AbstractMatrix{T}) where T<:Real
    k = size(X, 1)
    @inbounds for j = 1:k, i = 1:j  # Iterate the upper triangle
        if i == j
            X[i,i] = zero(T)
        else
            X[i,j], X[j,i] = F[i,j] * (X[i,j] - X[j,i]), F[j,i] * (X[j,i] - X[i,j])
        end
    end
    return X
end
_mulsubtrans!!(X::AbstractZero, F::AbstractZero) = X
_mulsubtrans!!(X::AbstractZero, F::AbstractMatrix{<:Real}) = X
_mulsubtrans!!(X::AbstractMatrix{<:Real}, F::AbstractZero) = F

# I - X, overwrites X
function _eyesubx!(X::AbstractMatrix)
    n, m = size(X)
    @inbounds for j = 1:m, i = 1:n
        X[i,j] = (i == j) - X[i,j]
    end
    return X
end

_extract_imag(x) = complex(0, imag(x))


text:
module FractionalCalculus

using QuadGK
using SpecialFunctions: gamma
using LinearAlgebra, InvertedIndices

include("./Derivative/Derivative.jl")
include("./Integral/Integral.jl")

export FracDiffAlg, FracIntAlg
export fracdiff, fracint

#Export fractional derivative releating API
export Caputo, GL, RLDiff
export Caputo_Direct, Caputo_Direct_First_Diff_Known, Caputo_Direct_First_Second_Diff_Known
export Caputo_Piecewise
export Caputo_Diethelm
export GL_Direct

export GL_Multiplicative_Additive, GL_Lagrange_Three_Point_Interp
export RLDiff_Approx
export RLDiff_Matrix

#Export fractional integral relating API
export RLInt
export RL_Direct, RL_Direct_First_Diff_Known, RL_Piecewise
export RLInt_Approx, RL_LinearInterp, RLInt_Matrix
end


text:
# Create a model with shared parameters ****************
#
export Model2Fit_Shared_Parameters
export get_model, get_model_θ

@doc raw"""

Share a set of parameters

# Example

Let's assume that we have an initial model `model` with `θ` as parameter vector.
If we want θ1, θ3, θ5 to share a same value θshared then:

```julia
indices_to_share = [1,3,5]

model_with_shared_params = Model2Fit_Shared_Parameters(model, indices_to_share)
```

Do not forget to udpate parameter vector θ, this can be done as
follows:

```julia
deleteat!(θ, indices_to_share)
push!(θ, θshared)
```

# Also see
- [`Model2Fit_Mapped_Parameters`](@ref)
"""
struct Model2Fit_Shared_Parameters{MODEL <: Abstract_Model2Fit,
                                   INDICES <: AbstractVector{Int}}  <: Abstract_Model2Fit
    _model::MODEL
    _indices::INDICES

    function Model2Fit_Shared_Parameters(model::MODEL,
                                         indices::INDICES) where {MODEL <: Abstract_Model2Fit,
                                                                  INDICES <: AbstractVector{Int}}
        @assert (length(indices) == 0) || ( (1 ≤ minimum(indices)) && (maximum(indices) ≤ parameter_size(model)) )
        @assert allunique(indices)

        new{MODEL,INDICES}(model, indices)

    end

end

# Convenience methods ================
#
@doc raw"""
```julia
get_model(mp::Model2Fit_Shared_Parameters) -> Absatrct_Model2Fit
```

Get back the wrapped model
"""
get_model(mp::Model2Fit_Shared_Parameters) = mp._model

@doc raw"""
```julia
get_model_θ(mp::Model2Fit_Shared_Parameters,θ::AbstractVector) -> θ::AbstractVector
```

Retrieve the parameter vector θ associated to the wrapped model [`get_model`](@ref).
"""
function get_model_θ(model::Model2Fit_Shared_Parameters,θ::AbstractVector)
    # Count unshared parameters
    #
    θ_not_shared_size = parameter_size(model._model) - length(model._indices)
    θ_not_shared = @view θ[1:θ_not_shared_size]

    # Retrieve shared value ================
    #
    # This is the last parameter
    #
    @assert θ_not_shared_size+1 == length(θ)
    θ_shared = θ[end]

    # Create the new θ with θ_shared set a [indices]
    #
    insert_some_elements(θ_not_shared, model._indices, θ_shared)
end

# Visit  ================
#
# note: no need to check
#   @assert 1 ≤ submodel_idx ≤ visit_submodel_size(model)
#   @assert length(X) == length(Y)
# as this is done higher in the visit() method
#
visit_submodel_size(model::Model2Fit_Shared_Parameters) = 1
visit_get_submodel(model::Model2Fit_Shared_Parameters,submodel_idx::Int) = model._model
visit_get_X(model::Model2Fit_Shared_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = X
visit_get_Y(model::Model2Fit_Shared_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = Y
visit_get_θ(model::Model2Fit_Shared_Parameters,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = get_model_θ(model,θ)

# Interface ================
#
parameter_size(mp::Model2Fit_Shared_Parameters) = parameter_size(mp._model) - length(mp._indices) + 1

function accumulate_y!(mp::Model2Fit_Shared_Parameters,Y::AbstractVector,X::AbstractVector,θ::AbstractVector)
    @assert length(θ) == parameter_size(mp)

    θ_model = get_model_θ(mp,θ)
    model = get_model(mp)

    accumulate_y!(model,Y,X,θ_model)
end


text:
fun(x::Nanosecond) = x/10^6
fun(vf::VideoFile) = VideoFile(vf, (duration = fun(vf.duration), ))

fun(v::WholeVideo) = WholeVideo(v, (duration = fun(v.file.duration), ))
fun(v::FragmentedVideo) = FragmentedVideo(v, (files = [fun(file) for file in files(v)], ))
fun(v::DisjointVideo) = DisjointVideo(v, (files = [fun(file) for file in files(v)], ))

fun(t::Instantaneous) = Instantaneous(fun(t.anchor))
fun(t::Prolonged) = Prolonged(fun(t.anchor), fun(t.duration))

fun(x::Temporal) = Temporal(fun(x.video), fun(x.time), x.comment)

fun(c::Calibration) = Calibration(fun(c.intrinsic), fun(c.extrinsic), c.board, c.comment)

fun(_::Missing) = missing

function change2ms()

    mv(joinpath(sourcefolder, "video"), joinpath(sourcefolder, "video_ns"), force=true)
    vs = deserialize(joinpath(sourcefolder, "video_ns"))
    map!(fun, vs, vs)
    serialize(joinpath(sourcefolder, "video"), vs)

    mv(joinpath(sourcefolder, "calibration"), joinpath(sourcefolder, "calibration_ns"), force=true)
    cs = deserialize(joinpath(sourcefolder, "calibration_ns"))
    map!(fun, cs, cs)
    serialize(joinpath(sourcefolder, "calibration"), cs)

end


text:
"""
    UnitVector{T}(N, k)

Helper vector type of length `N` where the `k`th element is `one(T)`
and all the others `zero(T)`.
"""
struct UnitVector{T} <: AbstractVector{T}
    N::Int
    k::Int
end

Base.size(e::UnitVector) = (e.N,)
Base.getindex(e::UnitVector{T}, i::Int) where T = i == e.k ? one(T) : zero(T)
Base.show(io::IO, e::UnitVector{T}) where T = write(io, "ê{$T}($(e.k)|$(e.N))")


text:
using DiffEqBase, Test
using Distributions

@test DiffEqBase.promote_tspan((0.0,1.0)) == (0.0,1.0)
@test DiffEqBase.promote_tspan((0,1.0)) == (0.0,1.0)
@test DiffEqBase.promote_tspan(1.0) == (0.0,1.0)
@test DiffEqBase.promote_tspan(nothing) == (nothing,nothing)

prob = ODEProblem((u,p,t)->u,(p,t0)->p[1],(p)->(0.0,p[2]),(2.0,1.0))
prob2 = DiffEqBase.get_concrete_problem(prob,nothing)

@test prob2.u0 == 2.0
@test prob2.tspan == (0.0,1.0)

prob = ODEProblem((u,p,t)->u,(p,t)->Normal(p,1),(0.0,1.0),1.0)
prob2 = DiffEqBase.get_concrete_problem(prob,nothing)
@test typeof(prob2.u0) == Float64

prob = ODEProblem((u,p,t)->u,1.0,(0,1))
@test_logs (:warn,
            "Integer time values are incompatible with adaptive integrators. Utilize " *
            "floating point numbers instead of integers in this case, i.e. (0.0,1.0) " *
            "instead of (0,1).") DiffEqBase.adaptive_warn(prob.u0,prob.tspan)


text:
function load_data(case::Int64)::Tuple{ms_type,mp_type,ps_type,pp_type,u_type}

    df_sets = CSV.read("$(pwd())/data/df_sets.csv" ;delim=",",use_mmap=false)
    P  = 1:df_sets[1,:P]
    G  = 1:df_sets[1,:G]
    B  = 1:df_sets[1,:B]
    R  = 1:df_sets[1,:R]
    S  = 1:df_sets[1,:S]
    Hs = 1:df_sets[1,:H]
    G0 = 1:df_sets[1,:G0]
    gn = df_sets[1,:gn]
    b0 = df_sets[1,:b0]
    r0 = df_sets[1,:r0]

    df_unc_sets = CSV.read("$(pwd())/data/df_unc_sets_case$(case).csv" ;delim=",",use_mmap=false)
    NI1 = df_unc_sets[1,:NΩ1]
    NI2 = df_unc_sets[1,:NΩ2]
    NI3 = df_unc_sets[1,:NΩ3]
    NI0 = NI1 + NI2
    NI  = NI2 + NI3

    df_unc_params = CSV.read("$(pwd())/data/df_unc_params_case$(case).csv" ;delim=",",use_mmap=false)
    H      = zeros(NI,2)
    C      = zeros(NI,2)
    H[:,1] = df_unc_params[:,:ρco2 ]
    H[:,2] = -df_unc_params[:,:νD   ]
    C[:,1] = df_unc_params[:,:c_co2]
    C[:,2] = df_unc_params[:,:c_ur ]

    df_inv_params_P = CSV.read("$(pwd())/data/df_inv_params_P.csv" ;delim=",",use_mmap=false)
    c_inv_0   = df_inv_params_P[:,:c_inv_0]
    c_inv_5   = df_inv_params_P[:,:c_inv_5  ]
    c_fix     = df_inv_params_P[:,:c_fix    ].*10^3
    x_hist_5  = df_inv_params_P[:,:x_hist_5 ]
    x_hist_10 = df_inv_params_P[:,:x_hist_10]
    x_max     = df_inv_params_P[:,:x_max    ].*10^-3
    c_inv     = hcat(repeat(c_inv_0 ,1,NI1),repeat(c_inv_5  ,1,NI2)).*10^3
    x_hist    = hcat(repeat(x_hist_5,1,NI2),repeat(x_hist_10,1,NI3)).*10^-3
    πi0       = vcat(1/NI1*ones(NI1),1/NI2*ones(NI2))
    πi        = vcat(1/NI2*ones(NI2),1/NI3*ones(NI3))
    itoi0     = Vector{Array{Int64,1}}()  # [ωM,..] = ΩStoM[ωS]
    for i in 1:NI2 push!(itoi0,[1]) end
    for i in (NI2+1):(NI2+NI3) push!(itoi0,[1,Int64(((i-NI2-1)-(i-NI2-1)%NI2)/NI2+1)+1]) end

    df_inv_params_O = CSV.read("$(pwd())/data/df_inv_params_O.csv" ;delim=",",use_mmap=false)
    κ = df_inv_params_O[1,:κ]

    df_oper_params_G = CSV.read("$(pwd())/data/df_oper_params_G.csv" ;delim=",",use_mmap=false)
    c_OMvarG = df_oper_params_G[:,:c_varOM]
    c_fuelG  = df_oper_params_G[:,:c_fuel ]
    em_co2G  = df_oper_params_G[:,:em_co2 ]
    ηG       = df_oper_params_G[:,:η      ]
    rampG    = df_oper_params_G[:,:ramp   ]

    df_oper_params_B = CSV.read("$(pwd())/data/df_oper_params_B.csv" ;delim=",",use_mmap=false)
    ηB = df_oper_params_B[:,:η]
    PB = df_oper_params_B[:,:P]

    df_oper_params_R = CSV.read("$(pwd())/data/df_oper_params_R.csv" ;delim=",",use_mmap=false)
    PR = zeros(R[end],S[end],Hs[end])
    PR[1,1,:] = df_oper_params_R[:,:P1_S1]
    PR[1,2,:] = df_oper_params_R[:,:P1_S2]
    PR[1,3,:] = df_oper_params_R[:,:P1_S3]
    PR[1,4,:] = df_oper_params_R[:,:P1_S4]
    PR[2,1,:] = df_oper_params_R[:,:P2_S1]
    PR[2,2,:] = df_oper_params_R[:,:P2_S2]
    PR[2,3,:] = df_oper_params_R[:,:P2_S3]
    PR[2,4,:] = df_oper_params_R[:,:P2_S4]
    PR[3,1,:] = df_oper_params_R[:,:P3_S1]
    PR[3,2,:] = df_oper_params_R[:,:P3_S2]
    PR[3,3,:] = df_oper_params_R[:,:P3_S3]
    PR[3,4,:] = df_oper_params_R[:,:P3_S4]

    df_oper_params_D = CSV.read("$(pwd())/data/df_oper_params_D.csv" ;delim=",",use_mmap=false)
    PD = zeros(S[end],Hs[end])
    PD[1,:] = df_oper_params_D[:,:P_S1]
    PD[2,:] = df_oper_params_D[:,:P_S2]
    PD[3,:] = df_oper_params_D[:,:P_S3]
    PD[4,:] = df_oper_params_D[:,:P_S4]

    df_oper_params_O = CSV.read("$(pwd())/data/df_oper_params_O.csv" ;delim=",",use_mmap=false)
    νD      = df_oper_params_O[1,:νD     ]
    c_shed  = df_oper_params_O[1,:c_shed ]
    α       = df_oper_params_O[1,:α      ]
    c_co2   = df_oper_params_O[1,:c_co2  ]
    co2_lim = df_oper_params_O[1,:co2_lim]
    ρco2    = df_oper_params_O[1,:ρco2   ]

    ms  = ms_type(P,1:NI0,1:NI)
    mp  = mp_type(κ,x_hist,x_max,c_inv,c_fix,πi0,πi,itoi0)
    ps  = ps_type(P,G,B,R,S,Hs,G0,gn,b0,r0)
    pp  = pp_type(c_OMvarG,c_fuelG,em_co2G,ηG,rampG,ηB,PB,PR,PD,νD,c_shed,c_co2,co2_lim,ρco2,α)
    unc = u_type(NI,P[end]+size(H)[2],P[end],size(H)[2],size(C)[2],H,C)

    return ms,mp,ps,pp,unc
end


text:
using smallPT
using Test

@testset "smallPT.jl" begin
    # Write your own tests here.
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots
h(x, y) = x^(x^y) - y

x = range(0, 0.11; length=1000)
y = range(0, 1.1; length=1000)
contour(x, y, h; levels=zeros(1), color=:red, colorbar=false)
plot!(; xlim=extrema(x), ylim=extrema(y), size=(400, 400))

# %%
using Plots

function f(x, y, n)
    z = x^y
    for i in 1:n
        z = x^z
    end
    z
end

x = range(0, 2.1, length=1000)
y = range(0, 10.1; length=1000)
plot(; colorbar=false)
for n in 1:4
    contour!(x, y, (x, y) -> f(x, y, n) - y; levels=zeros(1), c=n, ls=:auto)
end
plot!(; xlim=extrema(x), ylim=extrema(y), size=(400, 400))

# %%

max_stars_repo_path: 0022/implicit plot.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0022/implicit plot.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0022/implicit plot.jl

max_forks_repo_name: genkuroki/public


text:
using OrdinaryDiffEq, ODEInterfaceDiffEq, Sundials, DiffEqDevTools, LSODA
using LinearAlgebra
using Plots
gr()

const T = Float64
abstract type AbstractFilamentCache end
abstract type AbstractMagneticForce end
abstract type AbstractInextensibilityCache end
abstract type AbstractSolver end
abstract type AbstractSolverCache end

struct FerromagneticContinuous <: AbstractMagneticForce
    ω :: T
    F :: Vector{T}
end

mutable struct FilamentCache{
        MagneticForce        <: AbstractMagneticForce,
        InextensibilityCache <: AbstractInextensibilityCache,
        SolverCache          <: AbstractSolverCache
            } <: AbstractFilamentCache
    N  :: Int
    μ  :: T
    Cm :: T
    x  :: SubArray{T,1,Vector{T},Tuple{StepRange{Int,Int}},true}
    y  :: SubArray{T,1,Vector{T},Tuple{StepRange{Int,Int}},true}
    z  :: SubArray{T,1,Vector{T},Tuple{StepRange{Int,Int}},true}
    A  :: Matrix{T}
    P  :: InextensibilityCache
    F  :: MagneticForce
    Sc :: SolverCache
end

struct NoHydroProjectionCache <: AbstractInextensibilityCache
    J         :: Matrix{T}
    P         :: Matrix{T}
    J_JT      :: Matrix{T}
    J_JT_LDLT :: LinearAlgebra.LDLt{T, SymTridiagonal{T}}
    P0        :: Matrix{T}

    NoHydroProjectionCache(N::Int) = new(
        zeros(N, 3*(N+1)),          # J
        zeros(3*(N+1), 3*(N+1)),    # P
        zeros(N,N),                 # J_JT
        LinearAlgebra.LDLt{T,SymTridiagonal{T}}(SymTridiagonal(zeros(N), zeros(N-1))),
        zeros(N, 3*(N+1))
    )
end

struct DiffEqSolverCache <: AbstractSolverCache
    S1 :: Vector{T}
    S2 :: Vector{T}

    DiffEqSolverCache(N::Integer) = new(zeros(T,3*(N+1)), zeros(T,3*(N+1)))
end

function FilamentCache(N=20; Cm=32, ω=200, Solver=SolverDiffEq)
    InextensibilityCache = NoHydroProjectionCache
    SolverCache = DiffEqSolverCache
    tmp = zeros(3*(N+1))
    FilamentCache{FerromagneticContinuous, InextensibilityCache, SolverCache}(
        N, N+1, Cm, view(tmp,1:3:3*(N+1)), view(tmp,2:3:3*(N+1)), view(tmp,3:3:3*(N+1)),
        zeros(3*(N+1), 3*(N+1)), # A
        InextensibilityCache(N), # P
        FerromagneticContinuous(ω, zeros(3*(N+1))),
        SolverCache(N)
    )
end

function stiffness_matrix!(f::AbstractFilamentCache)
    N, μ, A = f.N, f.μ, f.A
    @inbounds for j in axes(A, 2), i in axes(A, 1)
      A[i, j] = j == i ? 1 : 0
    end
    @inbounds for i in 1 : 3
        A[i,i] =    1
        A[i,3+i] = -2
        A[i,6+i] =  1

        A[3+i,i]   = -2
        A[3+i,3+i] =  5
        A[3+i,6+i] = -4
        A[3+i,9+i] =  1

        A[3*(N-1)+i,3*(N-3)+i] =  1
        A[3*(N-1)+i,3*(N-2)+i] = -4
        A[3*(N-1)+i,3*(N-1)+i] =  5
        A[3*(N-1)+i,3*N+i]     = -2

        A[3*N+i,3*(N-2)+i]     =  1
        A[3*N+i,3*(N-1)+i]     = -2
        A[3*N+i,3*N+i]         =  1

        for j in 2 : N-2
            A[3*j+i,3*j+i]     =  6
            A[3*j+i,3*(j-1)+i] = -4
            A[3*j+i,3*(j+1)+i] = -4
            A[3*j+i,3*(j-2)+i] =  1
            A[3*j+i,3*(j+2)+i] =  1
        end
    end
    rmul!(A, -μ^4)
    nothing
end

function update_separate_coordinates!(f::AbstractFilamentCache, r)
    N, x, y, z = f.N, f.x, f.y, f.z
    @inbounds for i in 1 : length(x)
        x[i] = r[3*i-2]
        y[i] = r[3*i-1]
        z[i] = r[3*i]
    end
    nothing
end

function update_united_coordinates!(f::AbstractFilamentCache, r)
    N, x, y, z = f.N, f.x, f.y, f.z
    @inbounds for i in 1 : length(x)
        r[3*i-2] = x[i]
        r[3*i-1] = y[i]
        r[3*i]   = z[i]
    end
    nothing
end

function update_united_coordinates(f::AbstractFilamentCache)
    r = zeros(T, 3*length(f.x))
    update_united_coordinates!(f, r)
    r
end

function initialize!(initial_conf_type::Symbol, f::AbstractFilamentCache)
    N, x, y, z = f.N, f.x, f.y, f.z
    if initial_conf_type == :StraightX
        x .= range(0, stop=1, length=N+1)
        y .= 0
        z .= 0
    else
        error("Unknown initial configuration requested.")
    end
    update_united_coordinates(f)
end

function magnetic_force!(::FerromagneticContinuous, f::AbstractFilamentCache, t)
    # TODO: generalize this for different magnetic fields as well
    N, μ, Cm, ω, F = f.N, f.μ, f.Cm, f.F.ω, f.F.F
    F[1]         = -μ * Cm * cos(ω*t)
    F[2]         = -μ * Cm * sin(ω*t)
    F[3*(N+1)-2] =  μ * Cm * cos(ω*t)
    F[3*(N+1)-1] =  μ * Cm * sin(ω*t)
    nothing
end

struct SolverDiffEq <: AbstractSolver end

function (f::FilamentCache)(dr, r, p, t)
    @views f.x, f.y, f.z = r[1:3:end], r[2:3:end], r[3:3:end]
    jacobian!(f)
    projection!(f)
    magnetic_force!(f.F, f, t)
    A, P, F, S1, S2 = f.A, f.P.P, f.F.F, f.Sc.S1, f.Sc.S2

    # implement dr = P * (A*r + F) in an optimized way to avoid temporaries
    mul!(S1, A, r)
    S1 .+= F
    mul!(S2, P, S1)
    copyto!(dr, S2)
    return dr
end

function jacobian!(f::FilamentCache)
    N, x, y, z, J = f.N, f.x, f.y, f.z, f.P.J
    @inbounds for i in 1 : N
        J[i, 3*i-2]     = -2 * (x[i+1]-x[i])
        J[i, 3*i-1]     = -2 * (y[i+1]-y[i])
        J[i, 3*i]       = -2 * (z[i+1]-z[i])
        J[i, 3*(i+1)-2] =  2 * (x[i+1]-x[i])
        J[i, 3*(i+1)-1] =  2 * (y[i+1]-y[i])
        J[i, 3*(i+1)]   =  2 * (z[i+1]-z[i])
    end
    nothing
end

function projection!(f::FilamentCache)
    # implement P[:] = I - J'/(J*J')*J in an optimized way to avoid temporaries
    J, P, J_JT, J_JT_LDLT, P0 = f.P.J, f.P.P, f.P.J_JT, f.P.J_JT_LDLT, f.P.P0
    mul!(J_JT, J, J')
    LDLt_inplace!(J_JT_LDLT, J_JT)
    ldiv!(P0, J_JT_LDLT, J)
    mul!(P, P0', J)
    subtract_from_identity!(P)
    nothing
end

function subtract_from_identity!(A)
    lmul!(-1, A)
    @inbounds for i in 1 : size(A,1)
        A[i,i] += 1
    end
    nothing
end

function LDLt_inplace!(L::LinearAlgebra.LDLt{T,SymTridiagonal{T}}, A::Matrix{T}) where {T<:Real}
    n = size(A,1)
    dv, ev = L.data.dv, L.data.ev
    @inbounds for (i,d) in enumerate(diagind(A))
        dv[i] = A[d]
    end
    @inbounds for (i,d) in enumerate(diagind(A,-1))
        ev[i] = A[d]
    end
    @inbounds @simd for i in 1 : n-1
        ev[i]   /= dv[i]
        dv[i+1] -= abs2(ev[i]) * dv[i]
    end
    L
end

function run(::SolverDiffEq; N=20, Cm=32, ω=200, time_end=1., solver=TRBDF2(autodiff=false), reltol=1e-6, abstol=1e-6)
    f = FilamentCache(N, Solver=SolverDiffEq, Cm=Cm, ω=ω)
    r0 = initialize!(:StraightX, f)
    stiffness_matrix!(f)
    prob = ODEProblem(ODEFunction(f, jac=(J, u, p, t)->(mul!(J, f.P.P, f.A); nothing)), r0, (0., time_end))
    sol = solve(prob, solver, dense=false, reltol=reltol, abstol=abstol)
end

sol = run(SolverDiffEq())
plot(sol,vars = (0,25))

N=20
f = FilamentCache(N, Solver=SolverDiffEq)
r0 = initialize!(:StraightX, f)
stiffness_matrix!(f)
prob = ODEProblem(f, r0, (0., 0.01))

sol = solve(prob, Vern9(), reltol=1e-14, abstol=1e-14)
test_sol = TestSolution(sol);

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => Rosenbrock23(autodiff=false)),
    Dict(:alg => Rodas4(autodiff=false)),
    Dict(:alg => radau()),
    Dict(:alg=>Exprb43(autodiff=false)),
    Dict(:alg=>Exprb32(autodiff=false)),
    Dict(:alg=>ImplicitEulerExtrapolation(autodiff=false)),
    Dict(:alg=>ImplicitDeuflhardExtrapolation(autodiff=false)),
    Dict(:alg=>ImplicitHairerWannerExtrapolation(autodiff=false)),
    ];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

sol = solve(prob, EPIRK4s3B(autodiff=false), dt=2^-3)

abstols=1 ./10 .^(3:5)
reltols=1 ./10 .^(3:5)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => HochOst4(),:dts=>2.0.^(-3:-1:-5)),
    Dict(:alg => EPIRK4s3B(),:dts=>2.0.^(-3:-1:-5)),
    Dict(:alg => EXPRB53s3(),:dts=>2.0.^(-3:-1:-5)),
    ];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => BS3()),
    Dict(:alg => Tsit5()),
    Dict(:alg => ImplicitEuler(autodiff=false)),
    Dict(:alg => Trapezoid(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => rodas()),
    Dict(:alg => dop853()),
    Dict(:alg => lsoda()),
    Dict(:alg => ROCK2()),
    Dict(:alg => ROCK4()),
    Dict(:alg => ESERK5())
    ];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => ImplicitEuler(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => ABDF2(autodiff=false)),
    Dict(:alg => QNDF(autodiff=false)),
    Dict(:alg => RadauIIA5(autodiff=false)),
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => CVODE_BDF(linear_solver=:GMRES)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false,linsolve=LinSolveGMRES())),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false,linsolve=LinSolveGMRES())),
];

names = [
    "CVODE-BDF",
    "CVODE-BDF (GMRES)",
    "TRBDF2",
    "TRBDF2 (GMRES)",
    "KenCarp4",
    "KenCarp4 (GMRES)",
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; names=names, appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => Trapezoid(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => rodas()),
    Dict(:alg => lsoda()),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => ROCK2()),
    Dict(:alg => ROCK4()),
    Dict(:alg => ESERK5())
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => ROCK2()),
    Dict(:alg => ROCK4()),
    Dict(:alg => ESERK5())
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false, dense_errors = true, error_estimate=:L2)
plot(wp)

abstols=1 ./10 .^(6:12)
reltols=1 ./10 .^(6:12)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => Vern7()),
    Dict(:alg => Vern9()),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => dop853()),
    Dict(:alg => ROCK4())
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(6:12)
reltols=1 ./10 .^(6:12)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => radau()),
    Dict(:alg => RadauIIA5(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno5(autodiff=false)),
    Dict(:alg => KenCarp5(autodiff=false)),
    Dict(:alg => lsoda()),
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                                    maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(6:12)
reltols=1 ./10 .^(6:12)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => radau()),
    Dict(:alg => RadauIIA5(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno5(autodiff=false)),
    Dict(:alg => KenCarp5(autodiff=false)),
    Dict(:alg => lsoda()),
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false, error_estimate = :l2)
plot(wp)

abstols=1 ./10 .^(6:12)
reltols=1 ./10 .^(6:12)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => radau()),
    Dict(:alg => RadauIIA5(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno5(autodiff=false)),
    Dict(:alg => KenCarp5(autodiff=false)),
    Dict(:alg => lsoda()),
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false, dense_errors=true, error_estimate = :L2)
plot(wp)

N=20
f = FilamentCache(N, Solver=SolverDiffEq)
r0 = initialize!(:StraightX, f)
stiffness_matrix!(f)
prob = ODEProblem(ODEFunction(f, jac=nothing), r0, (0., 0.01))

sol = solve(prob, Vern9(), reltol=1e-14, abstol=1e-14)
test_sol = TestSolution(sol.t, sol.u);

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => BS3()),
    Dict(:alg => Tsit5()),
    Dict(:alg => ImplicitEuler(autodiff=false)),
    Dict(:alg => Trapezoid(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => rodas()),
    Dict(:alg => dop853()),
    Dict(:alg => lsoda())
    ];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => BS3()),
    Dict(:alg => Tsit5()),
    Dict(:alg => ImplicitEuler(autodiff=false)),
    Dict(:alg => Trapezoid(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => rodas()),
    Dict(:alg => dop853()),
    Dict(:alg => lsoda()),
    Dict(:alg => ROCK2()),
    Dict(:alg => ROCK4()),
    Dict(:alg => ESERK5())
    ];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(3:8)
reltols=1 ./10 .^(3:8)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => CVODE_BDF(linear_solver=:GMRES)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false,linsolve=LinSolveGMRES())),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false,linsolve=LinSolveGMRES())),
];

names = [
    "CVODE-BDF",
    "CVODE-BDF (GMRES)",
    "TRBDF2",
    "TRBDF2 (GMRES)",
    "KenCarp4",
    "KenCarp4 (GMRES)",
];

wp = WorkPrecisionSet(prob, abstols, reltols, setups; names=names, appxsol=test_sol,
                      maxiters=Int(1e6), verbose = false)
plot(wp)

abstols=1 ./10 .^(6:12)
reltols=1 ./10 .^(6:12)
setups = [
    Dict(:alg => CVODE_BDF()),
    Dict(:alg => radau()),
    Dict(:alg => RadauIIA5(autodiff=false)),
    Dict(:alg => TRBDF2(autodiff=false)),
    Dict(:alg => Kvaerno3(autodiff=false)),
    Dict(:alg => KenCarp3(autodiff=false)),
    Dict(:alg => Kvaerno4(autodiff=false)),
    Dict(:alg => KenCarp4(autodiff=false)),
    Dict(:alg => Kvaerno5(autodiff=false)),
    Dict(:alg => KenCarp5(autodiff=false)),
    Dict(:alg => lsoda()),
];
wp = WorkPrecisionSet(prob, abstols, reltols, setups; appxsol=test_sol,
                                    maxiters=Int(1e6), verbose = false)
plot(wp)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using Test, Printf, ForwardDiff

using CLIMA.TurbulenceConvection.FiniteDifferenceGrids

@testset "Grid interface" begin
  for FT in (Float32, Float64)
    for n_ghost in (1, 2)
      n_elems = 12
      n_elems_real = n_elems-2*n_ghost
      elem_indexes = 1:n_elems
      elem_indexes_real = elem_indexes[1+n_ghost:end-n_ghost]
      Δz = FT(1/n_elems_real)
      grid = Grid(FT(0.0), FT(1.0), n_elems_real, n_ghost)
      @test n_hat(Zmin()) == -1
      @test n_hat(Zmax()) == 1
      @test binary(Zmin()) == 0
      @test binary(Zmax()) == 1
      @test ghost_dual(Zmin()) == [1, 0]
      @test ghost_dual(Zmax()) == [0, 1]
      @test eltype(grid) == FT
      @test all(over_elems(grid) .== elem_indexes)
      @test all(over_elems_real(grid) .== elem_indexes_real)
      @test length(over_elems(grid)) == n_elems
      @test length(over_elems_real(grid)) == n_elems_real
      @test first_interior(grid, Zmax()) == n_elems-n_ghost
      @test first_interior(grid, Zmin()) == 1+n_ghost
      @test boundary(grid, Zmin()) == 1+n_ghost
      @test boundary(grid, Zmax()) == n_elems-n_ghost
      @test over_elems_ghost(grid) == [(1:n_ghost)..., (n_elems+1-n_ghost:n_elems)...]
      @test grid.zn_min ≈ FT(0.0)
      @test grid.zn_max  ≈ FT(1.0)
      sprint(show, grid)
    end
  end
end

function ∇(f, x::T) where {T}
    tag = typeof(ForwardDiff.Tag(f, T))
    y = f(ForwardDiff.Dual{tag}(x,one(x)))
    ForwardDiff.partials(tag, y, 1)
end

@testset "Grid operators" begin
  for FT in (Float32, Float64)
    n_ghost = 1
    n_elems = 12
    n_elems_real = n_elems-2*n_ghost
    elem_indexes = 1:n_elems
    elem_indexes_real = elem_indexes[1+n_ghost:end-n_ghost]
    Δz = FT(1/n_elems_real)
    grid = Grid(FT(0.0), FT(1.0), n_elems_real, n_ghost)
    f = FT[1, 5, 3]
    w = FT[2, 3, 8]
    @test grad(f[1:2], grid) ≈ FT(40)
    @test grad(f, grid) ≈ FT(10)
    @test_throws AssertionError grad(f[1], grid)
    @test_throws AssertionError advect(f[1:2], w[1:2], grid)
    @test advect(f, w, grid) ≈ FT(40)
    @test advect(f, -w, grid) ≈ FT(-20)
    @test ∇_pos(f, grid) ≈ FT(40)
    @test ∇_neg(f, grid) ≈ FT(-20)
  end
end

@testset "Grid operators convergence" begin
  z_min = 0.0
  z_max = 2.0*π
  tol = 1/2
  ϕ_(x) = 5+sin(x)
  ∇ϕ_(x) = cos(x) # computed manually, since ForwardDiff doesn't like gradient functions as inputs.
  ψ_(x) = 5+cos(x)
  w_(x) = 5+sin(x)
  K_(x) = 5+sin(x)
  ϕψ_(x) = ψ_(x)*ϕ_(x)
  wϕ_(x) = w_(x)*ϕ_(x)
  K∇ϕ_(x) = K_(x)*∇ϕ_(x)

  for n_elems_real in [2^k for k in 3:8]
    grid = Grid(z_min, z_max, n_elems_real)
    z = grid.zc

    ϕ     = ϕ_.(z)
    ψ     = ψ_.(z)
    ϕψ    = ϕψ_.(z)
    w     = w_.(z)
    K     = K_.(z)

    ∇ϕ  = ∇.(ϕ_, z)
    ∇ψ  = ∇.(ψ_, z)
    ∇ϕψ = ∇.(ϕψ_, z)
    ∇w  = ∇.(w_, z)
    ∇K∇ϕ = ∇.(K∇ϕ_, z)
    ∇ϕ  = ∇.(ϕ_, z)
    ∇wϕ  = ∇.(wϕ_, z)

    del_err = [abs(∇ϕ[k] - ∇_z_upwind(ϕ[k-1:k+1], w[k-1:k+1], grid)) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz for e in del_err])

    del_err = [abs((∇ϕ[k-1]+∇ϕ[k])/2 - ∇_z_flux(ϕ[k-1:k], grid)) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz^2 for e in del_err])
    del_err = [abs((∇ϕ[k]+∇ϕ[k+1])/2 - ∇_z_flux(ϕ[k:k+1], grid)) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz^2 for e in del_err])

    del_err = [abs(∇ϕ[k] - ∇_z_centered(ϕ[k-1:k+1], grid)) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz^2 for e in del_err])

    del_err = [abs((∇ϕ[k-1]+∇ϕ[k])/2 - ∇_z_dual(ϕ[k-1:k+1], grid)[1]) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz^2 for e in del_err])
    del_err = [abs((∇ϕ[k]+∇ϕ[k+1])/2 - ∇_z_dual(ϕ[k-1:k+1], grid)[2]) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz^2 for e in del_err])

    Lap_err = [abs(∇ψ[k] - Δ_z(ϕ[k-1:k+1], grid)) for k in over_elems_real(grid)]
    @test all([e<tol*grid.Δz^2 for e in Lap_err])

    Lap_err = [abs(∇K∇ϕ[k] - Δ_z(ϕ[k-1:k+1], grid, K[k-1:k+1])) for k in over_elems_real(grid)]
    @test all([e<2*tol*grid.Δz^2 for e in Lap_err])

    Lap_err = [abs(∇K∇ϕ[k] - Δ_z_dual(ϕ[k-1:k+1], grid, 0.5*(K[k-1:k].+K[k:k+1]))) for k in over_elems_real(grid)]
    @test all([e<2*tol*grid.Δz^2 for e in Lap_err])
    Δt = 1.0

    # ∇•(wϕ)
    for conservatove_scheme in (
                                UpwindCollocated(),
                                OneSidedUp(),
                                OneSidedDn(),
                                CenteredUnstable(),
                                )
      adv_err = [abs(∇wϕ[k] - advect_old(ϕ[k-1:k+1],
                                     0.5*(ϕ[k-1:k].+ϕ[k:k+1]),
                                     w[k-1:k+1],
                                     0.5*(w[k-1:k].+w[k:k+1]),
                                     grid, conservatove_scheme, Δt)) for k in over_elems_real(grid)]
      @test all([e<12*tol*grid.Δz for e in adv_err])
    end
    # w•∇ϕ
    for non_conservative_scheme in (UpwindAdvective(),)
      adv_err = [abs(w[k]*∇ϕ[k] - advect_old(ϕ[k-1:k+1],
                                         0.5*(ϕ[k-1:k].+ϕ[k:k+1]),
                                         w[k-1:k+1],
                                         0.5*(w[k-1:k].+w[k:k+1]),
                                         grid, non_conservative_scheme, Δt)) for k in over_elems_real(grid)]
      @test all([e<6*tol*grid.Δz for e in adv_err])
    end

  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
export CubicSplineKernel, TabulatedCubicSplineKernel

"Cubic spline smoothing kernel."
struct CubicSplineKernel{N} <: KernelFormula{N}
    σ::Float

    CubicSplineKernel{1}() = new(2 / 3)
    CubicSplineKernel{2}() = new(10 / (7π))
    CubicSplineKernel{3}() = new(1 / π)
end

Base.show(io::IO, ::Type{CubicSplineKernel{N}}) where {N} =
    print(io, "CubicSplineKernel")

"""
    (W::CubicSplineKernel)(q::Number) -> Float

Evaluate the kernel with unit width at relative distance `q`.
"""
(W::CubicSplineKernel)(q::Number) =
    @dbgassert(0 ≤ q ≤ 2) &&
    W.σ * (q ≤ 1 ? (1 - 1.5 * q^2 * (1 - 0.5 * q)) : (0.25 * (2 - q)^3))

"""
    ddr(W::CubicSplineKernel, q::Number) -> Float

Evaluate the derivative of the kernel with unit width at relative distance `q`.
The derivative is with respect to `r`.
"""
ddr(W::CubicSplineKernel, q::Number) =
    @dbgassert(0 ≤ q ≤ 2) &&
    W.σ * (q ≤ 1 ? q * (2.25 * q - 3) : -0.75 * (2 - q)^2)

"""
    ddh(W::CubicSplineKernel, q::Number) -> Float

Evaluate the derivative of the kernel with unit width at relative distance `q`.
The derivative is with respect to `h`.
"""
ddh(W::CubicSplineKernel{N}, q::Number) where {N} =
    @dbgassert(0 ≤ q ≤ 2) &&
    W.σ * (
        q ≤ 1 ? ((3 + 1.5N) * q^2 - (2.25 + 0.75N) * q^3 - N) :
        (-0.5 * (2 - q)^2 * (N - (1.5 + 0.5 * N) * q))
    )

get_support(W::CubicSplineKernel) = 2.0

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
function My_sub2ind(n::Array{Int64},sub::Array{Int64})
if length(sub)==2
	return sub2ind(tuple(n...),sub[1],sub[2]);
else
	return sub2ind(tuple(n...),sub[1],sub[2],sub[3]);
end
end


text:
let
    num_notebooks_tested = 0
    notebookdir = joinpath(@__DIR__, "..", "examples")
    excludedirs = [".ipynb_checkpoints"]
    excludefiles = String[]
    push!(excludefiles, "6. Symbolics using SymPy.ipynb") # Manifest used for 1.1 doesn't work for 1.0.
    for (root, dir, files) in walkdir(notebookdir)
        basename(root) in excludedirs && continue
        for file in files
            file in excludefiles && continue
            name, ext = splitext(file)
            lowercase(ext) == ".ipynb" || continue
            path = joinpath(root, file)
            @eval module $(gensym()) # Each notebook is run in its own module.
            using Test
            using NBInclude
            @testset "Notebook: $($name)" begin
                # Note: use #NBSKIP in a cell to skip it during tests.
                @nbinclude($path; regex = r"^((?!\#NBSKIP).)*$"s)
            end
            end # module
            num_notebooks_tested += 1
        end
    end
    @test num_notebooks_tested > 0
end


text:
module DataTransform

  using ...Statistics
  #assume that you have a matrix m-n

  function convertToFloat(A,datatype=Float64)
    for i=1:length(A)
      val = A[i]
      #println(typeof(val))

      if typeof(val)==Int64
        A[i]=convert(datatype,val)
      elseif typeof(val)==ASCIIString || typeof(val)==SubString{ASCIIString}
        A[i]=parse(datatype,val)
      end
    end
    return convert(Array{datatype,ndims(A)}, A)

  end

  function addProxyColumns(A,proxyColumns)
    columns=[]
    for i in proxyColumns

      #for every column, create the proxy columns, find indicies for each, and yeah.
      col = i["column"]
      push!(columns,col)
      proxyentries = i["proxy"]

      proxy = zeros(Int, size(A)[1],length(proxyentries))
      column = A[:,col]
      for j=1:length(proxyentries)
        proxy[find(column.==proxyentries[j]),j]=1
      end
      A=[A proxy];
    end
    allcolumns = collect(1:size(A)[2])
    return A[:,setdiff(allcolumns,columns)];
  end

  function zScale(A,mean=nothing, std=nothing)
    if mean==nothing || std==nothing
      mean=Statistics.meanColumn(A)
      std = Statistics.stdColumn(A)
    end
    std= 1./std
    #B= convert(Array{Float64,2}, A)
    for i=1:size(A)[1]
        A[i,:]= (A[i,:]-mean).*std
    end
    return A

  end

  function zScaleColumn(A,col,average=nothing, standardeviation=nothing)
    column = A[:,col]
    if average==nothing || standardeviation==nothing
      average=mean(column)
      standardeviation = std(column)
    end
    A[:,col]= (column - average)/standardeviation
    return A
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using Yao
using Yao.Intrinsics
using Compat

######################### CONTINUOUS ################################
uniform_state(num_bit) = register(ones(Complex128, 1<<num_bit)/sqrt(1<<num_bit))

######################### PEAKS ################################
"""
GHZ wave function.
"""
function ghz(num_bit::Int; x::Int=0)
    v = zeros(ComplexF64, 1<<num_bit)
    v[x+1] = 1/sqrt(2)
    v[flip(x, bmask(1:num_bit))+1] = 1/sqrt(2)
    return v
end

"""
onehot wave function.
"""
function onehot(num_bit::Int, x::Int)
    v = zeros(ComplexF64, 1<<num_bit)
    v[x+1] = 1
    return v
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
using Test

macro test_error(ex)
    @gensym err
    quote
        let $err = nothing
            $Test.@test try
                $ex
                false
            catch _err
                $err = _err
                true
            end
            $err
        end
    end |> esc
end

const ⊏ = occursin


text:
function Clifford(x::Float64, y::Float64, a::Float64, b::Float64, c::Float64, d::Float64)
    sin(a * y) + c * cos(a * x), sin(b * x) + d * cos(b * y)
end
function De_Jong(x::Float64, y::Float64, a::Float64, b::Float64, c::Float64, d::Float64)
    sin(a * y) - cos(b * x),sin(c * x) - cos(d * y)
end
function Svensson(x::Float64, y::Float64, a::Float64, b::Float64, c::Float64, d::Float64)
    d * sin(a * x) - sin(b * y), c * cos(a * x) + cos(b * y)
end
function Bedhead(x::Float64, y::Float64, a::Float64, b::Float64, c::Float64, d::Float64)
    sin(x*y/b)*y + cos(a*x-y), x + sin(y)/b
end
function Fractal_Dream(x::Float64, y::Float64, a::Float64, b::Float64, c::Float64, d::Float64)
    sin(y*b)+c*sin(x*b), sin(x*a) + d*sin(y*a)
end
function trajectory(fn, x0::Float64, y0::Float64, a::Float64, b::Float64, c::Float64, d::Float64,
        dθ::Float64,  n::Int64)
    x, y, θ = zeros(n), zeros(n), 0.0
    x[1], y[1] = x0, y0
    for i = 1:n
        xd, yd = fn(x[i], y[i], a, b, c, d)
        @inbounds x[i+1], y[i+1] = xd*cos(θ), yd*cos(θ)
        θ += dθ
    end
    x, y
end
#Examples
cliffordExamples = [
    [-1.3, -1.3, -1.8, -1.9],
    [-1.4, 1.6, 1.0, 0.7],
    [1.7, 1.7, 0.6, 1.2],
    [1.7, 0.7, 1.4, 2.0],
    [-1.6, 1.6, 0.7, -1.0],
    [-1.32, -1.65, 0.74, 1.81]]
De_JongExamples = [
    [-1.244, -1.251, -1.815, -1.908],
    [1.4, -2.3, 2.4, -2.1],
    [1.4, 1.56, 1.4, -6.56],
    [2.01, -2.53, 1.61, -0.33],
    [-0.827, -1.637, 1.659, -0.943],
    [-0.709, 1.638, 0.452, 1.740]]
SvenssonExamples = [
    [1.5, -1.8, 1.6, 0.9],
    [-1.78, 1.29, -0.09, -1.18],
    [-0.91, -1.29, -1.97, -1.56],
    [1.4, 1.56, 1.4, -6.56],
    [-0.827, -1.637, 1.659, -0.943],
    [1.7, 1.7, 0.6, 1.2]]
BedheadExamples = [
    [-0.81, -0.92, 0.0, 0.0],
    [-0.64, 0.76, 0.0, 0.0],
    [0.06, 0.98, 0.0, 0.0],
    [-0.67, 0.83, 0.0, 0.0],
    [0.65343, 0.7345345, 0.0, 0.0],
    [1.7, 1.7, 0.0, 0.0]]

attractorExamples = OrderedDict(Clifford => cliffordExamples,
    De_Jong => De_JongExamples, Svensson => SvenssonExamples,
    Bedhead => BedheadExamples)

function pkgsVersion()
    println("Julia == " * string(VERSION))
    for (key, version) ∈ sort(collect(Pkg.installed()))
        try
            isa(eval(Symbol(key)), Module) && println(key * " == " * string(version))
        catch
        end
    end
end

max_stars_repo_name: lazarusA/CodeSnippets

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
# An implementation of CRMR for the solution of the
# (under/over-determined or square) linear system
#
#  Ax = b.
#
# The method seeks to solve the minimum-norm problem
#
#  min ‖x‖²  s.t. Ax = b,
#
# and is equivalent to applying the conjugate residual method
# to the linear system
#
#  AA'y = b.
#
# This method is equivalent to Craig-MR, described in
#
# M. Arioli and D. Orban, Iterative Methods for Symmetric
# Quasi-Definite Linear Systems, Part I: Theory.
# Cahier du GERAD G-2013-32, GERAD, Montreal QC, Canada, 2013.
#
# D. Orban, The Projected Golub-Kahan Process for Constrained
# Linear Least-Squares Problems. Cahier du GERAD G-2014-15,
# GERAD, Montreal QC, Canada, 2014.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Montreal, QC, April 2015.

export crmr

"""Solve the consistent linear system

  Ax + √λs = b

using the Conjugate Residual (CR) method, where λ ≥ 0 is a regularization
parameter. This method is equivalent to applying CR to the normal equations
of the second kind

  (AA' + λI) y = b

but is more stable. When λ = 0, this method solves the minimum-norm problem

  min ‖x‖₂  s.t.  x ∈ argmin ‖Ax - b‖₂.

When λ > 0, this method solves the problem

  min ‖(x,s)‖₂  s.t. Ax + √λs = b.

CGMR produces monotonic residuals ‖r‖₂.
It is formally equivalent to CRAIG-MR, though can be slightly less accurate,
but simpler to implement. Only the x-part of the solution is returned.

A preconditioner M may be provided in the form of a linear operator.
"""
function crmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(A,1)), λ :: Float64=0.0,
              atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
              itmax :: Int=0, verbose :: Bool=false) where T <: Number

  m, n = size(A);
  size(b, 1) == m || error("Inconsistent problem size");
  verbose && @printf("CRMR: system of %d equations in %d variables\n", m, n);

  x = zeros(T, n) # initial estimation x = 0
  r = copy(M * b) # initial residual r = M * (b - Ax) = M * b
  bNorm = @knrm2(m, r)  # norm(b - A * x0) if x0 ≠ 0.
  bNorm == 0 && return x, SimpleStats(true, false, [0.0], [0.0], "x = 0 is a zero-residual solution");
  rNorm = bNorm;  # + λ * ‖x0‖ if x0 ≠ 0 and λ > 0.
  λ > 0 && (s = copy(r));
  Ar = A' * r;  # - λ * x0 if x0 ≠ 0.
  p  = copy(Ar);
  γ  = @kdot(n, Ar, Ar)  # Faster than γ = dot(Ar, Ar);
  λ > 0 && (γ += λ * rNorm * rNorm);
  iter = 0;
  itmax == 0 && (itmax = m + n);

  ArNorm = sqrt(γ);
  rNorms = [rNorm;];
  ArNorms = [ArNorm;];
  ɛ_c = atol + rtol * rNorm;   # Stopping tolerance for consistent systems.
  ɛ_i = atol + rtol * ArNorm;  # Stopping tolerance for inconsistent systems.
  verbose && @printf("%5s  %8s  %8s\n", "Aprod", "‖A'r‖", "‖r‖")
  verbose && @printf("%5d  %8.2e  %8.2e\n", 1, ArNorm, rNorm);

  status = "unknown";
  solved = rNorm ≤ ɛ_c;
  inconsistent = (rNorm > 1.0e+2 * ɛ_c) && (ArNorm ≤ ɛ_i);
  tired = iter ≥ itmax;

  while ! (solved || inconsistent || tired)
    q = A * p;
    λ > 0 && @kaxpy!(m, λ, s, q)  # q = q + λ * s;
    Mq = M * q
    α = γ / @kdot(m, q, Mq)    # Compute qᵗ * M * q
    @kaxpy!(n,  α, p, x)       # Faster than  x =  x + α *  p;
    @kaxpy!(m, -α, Mq, r)      # Faster than  r =  r - α * Mq;
    rNorm = @knrm2(m, r)       # norm(r);
    Ar = A' * r;
    γ_next = @kdot(n, Ar, Ar)  # Faster than γ_next = dot(Ar, Ar);
    λ > 0 && (γ_next += λ * rNorm * rNorm);
    β = γ_next / γ;

    @kaxpby!(n, 1.0, Ar, β, p)  # Faster than  p = Ar + β * p
    if λ > 0
      @kaxpby!(m, 1.0, r, β, s) # s = r + β * s
    end

    γ = γ_next;
    ArNorm = sqrt(γ);
    push!(rNorms, rNorm);
    push!(ArNorms, ArNorm);
    iter = iter + 1;
    verbose && @printf("%5d  %8.2e  %8.2e\n", 1 + 2 * iter, ArNorm, rNorm);
    solved = rNorm ≤ ɛ_c;
    inconsistent = (rNorm > 1.0e+2 * ɛ_c) && (ArNorm ≤ ɛ_i);
    tired = iter ≥ itmax;
  end

  status = tired ? "maximum number of iterations exceeded" : (inconsistent ? "system probably inconsistent but least squares/norm solution found" : "solution good enough given atol and rtol")
  stats = SimpleStats(solved, inconsistent, rNorms, ArNorms, status);
  return (x, stats);
end


text:
module BugReporting

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) BugReporting

export replay, make_interactive_report

using Base.Filesystem: uperm
using rr_jll
using Zstd_jll
using HTTP, JSON
using AWS, AWSS3
using Tar
using Pkg

# https://github.com/JuliaLang/julia/pull/29411
if isdefined(Base, :exit_on_sigint)
    using Base: exit_on_sigint
else
    exit_on_sigint(on::Bool) =
        ccall(:jl_exit_on_sigint, Cvoid, (Cint,), on)
end

const WSS_ENDPOINT = "wss://53ly7yebjg.execute-api.us-east-1.amazonaws.com/test"
const GITHUB_APP_ID = "Iv1.c29a629771fe63c4"
const TRACE_BUCKET = "julialang-dumps"

function check_rr_available()
    if !isdefined(rr_jll, :rr_path)
        error("RR not available on this platform")
    end
end

# Values that are initialized in `__init__()`
global_record_flags = String[]
ignore_child_status = false

struct InvalidPerfEventParanoidError <: Exception
    value
end

function Base.showerror(io::IO, err::InvalidPerfEventParanoidError)
    println(io, "InvalidPerfEventParanoidError")
    print(io, """
    rr needs /proc/sys/kernel/perf_event_paranoid <= 1, but it is $(err.value).
    Change it to 1, or use `JULIA_RR_RECORD_ARGS=-n julia --bug-report=rr` (slow).
    Consider putting 'kernel.perf_event_paranoid = 1' in /etc/sysctl.conf
    or change it temporarily by
        echo 1 | sudo tee /proc/sys/kernel/perf_event_paranoid
    """)
end

# `path` used for testing
function check_perf_event_paranoid(path = "/proc/sys/kernel/perf_event_paranoid"; record_flags = global_record_flags)
    isempty(intersect(["-n", "--no-syscall-buffer"], record_flags)) || return
    isfile(path) || return  # let `rr` handle this
    value = tryparse(Int, read(path, String))
    value === nothing && return  # let `rr` handle this
    value <= 1 && return
    throw(InvalidPerfEventParanoidError(value))
end

function collect_inner_traces(trace_directory)
    # If this is already a trace directory, return just it
    if isfile(joinpath(trace_directory, "version"))
        return [trace_directory]
    end
    # If this is a directory containing traces, return those
    traces = String[]
    for f in readdir(trace_directory, join=true)
        if isfile(joinpath(f, "version"))
            push!(traces, f)
        end
    end
    return unique(realpath.(traces))
end

function find_latest_trace(trace_directory)
    # What _RR_TRACE_DIR calls a "trace directory" is not what `rr pack` calls a "trace directory"
    # This function allows us to normalize to the inner "latest-trace" directory, if necessary.
    latest_symlink = joinpath(trace_directory, "latest-trace")
    if !isfile(joinpath(trace_directory, "version")) && islink(latest_symlink)
        return realpath(latest_symlink)
    end
    return trace_directory
end

function rr_pack(trace_directory)
    check_rr_available()

    rr() do rr_path
        for dir in collect_inner_traces(trace_directory)
            @debug("rr pack'ing $(dir)")
            run(`$rr_path pack $(dir)`)
        end
    end
end

# Helper function to compress a trace into a `.tar.zst` file
function compress_trace(trace_directory::String, output_file::String)
    # Ensure it's packed
    rr_pack(trace_directory)

    # Start up our `zstdmt` process to write out to that file
    proc = zstdmt() do zstdp
        open(pipeline(`$(zstdp) - -o $(output_file)`, stderr=stderr), "r+")
    end

    # Feed the tarball into that waiting process
    Tar.create(trace_directory, proc)

    # Ensure everything closes nicely
    close(proc.in)
    wait(proc)
    return nothing
end

function rr_record(args...; trace_dir=nothing)
    check_rr_available()
    check_perf_event_paranoid()

    rr() do rr_path
        new_env = copy(ENV)
        if trace_dir !== nothing
            new_env["_RR_TRACE_DIR"] = trace_dir
        end
        # Intersperse all given arguments with spaces, then splat:
        rr_cmd = `$(rr_path) record $(global_record_flags)`
        for arg in args
            rr_cmd = `$(rr_cmd) $(arg)`
        end
        return run(ignorestatus(setenv(rr_cmd, new_env)))
    end
end

function decompress_rr_trace(trace_file, out_dir)
    proc = zstdmt() do zstdp
        open(`$zstdp --stdout -d $trace_file`, "r+")
    end
    Tar.extract(proc, out_dir)
    return nothing
end

function decompress_rr_trace(trace_file)
    artifact_hash = Pkg.create_artifact() do dir
        decompress_rr_trace(trace_file, dir)
    end
    return Pkg.artifact_path(artifact_hash)
end

function download_rr_trace(trace_url; verbose=true)
    isdefined(Pkg.PlatformEngines, :probe_platform_engines) && Pkg.PlatformEngines.probe_platform_engines!()
    mktempdir() do dl_dir
        # Download into temporary directory, unpack into artifact directory
        local_path = joinpath(dl_dir, "trace.tar.zst")
        Pkg.PlatformEngines.download(trace_url, local_path; verbose=verbose)
        decompress_rr_trace(local_path)
    end
end

function replay(trace_url)
    if startswith(trace_url, "s3://")
        trace_url = string("https://s3.amazonaws.com/julialang-dumps/", trace_url[6:end])
    end
    if startswith(trace_url, "https://")
        trace_url = download_rr_trace(trace_url)
    end

    # If it's a file, try to decompress it
    if isfile(trace_url)
        trace_url = decompress_rr_trace(trace_url)
    end

    # If it's not a directory by now, we don't know what to do
    if !isdir(trace_url)
        error("Invalid trace location: $(trace_url)")
    end

    rr() do rr_path
        run(`$(rr_path) replay $(find_latest_trace(trace_url))`)
    end
end

function handle_child_error(p::Base.Process)
    # If the user has requested that we ignore child status, do so
    if ignore_child_status
        return
    end

    if !success(p)
        @info("Debugged process failed. Unless you see rr errors above, the trace likely completed.", exitcode=p.exitcode, termsignal=p.termsignal)

        # Return the exit code if that is nonzero
        if p.exitcode != 0
            exit(p.exitcode)
        end

        # If the child instead signalled, we recreate the same signal in ourselves
        ccall(:signal, Ptr{Cvoid}, (Cint, Ptr{Cvoid}), p.termsignal, C_NULL)
        ccall(:raise, Cint, (Cint,), p.termsignal)
    end
end

function make_interactive_report(report_type, ARGS=[])
    default_julia_args = `--history-file=no`
    if report_type == "rr-local"
        proc = rr_record(`$(Base.julia_cmd()) $default_julia_args`, ARGS)
        handle_child_error(proc)
        return
    elseif report_type == "rr"
        exit_on_sigint(false)  # throw InterruptException on Ctrl-C
        local proc
        artifact_hash = Pkg.create_artifact() do trace_dir
            proc = rr_record(`$(Base.julia_cmd()) $default_julia_args`, ARGS; trace_dir=trace_dir)
            @info "Preparing trace directory for upload (if your trace is large this may take a few minutes)"
            rr_pack(trace_dir)
        end
        upload_rr_trace(Pkg.artifact_path(artifact_hash))
        handle_child_error(proc)
        return
    elseif report_type == "help"
        show(stdout, "text/plain", @doc(BugReporting))
        println()
        return
    end
    error("Unknown report type: $report_type")
end

const S3_CHUNK_SIZE = 25 * 1024 * 1024 # 25 MB

include("sync_compat.jl")

function upload_rr_trace(trace_directory)
    # Auto-pack this trace directory if it hasn't already been:
    sample_directory = joinpath(trace_directory, "latest-trace")
    if isdir(sample_directory) && uperm(sample_directory) & 0x2 == 0
        @info "`$sample_directory` not writable. Skipping `rr pack`."
    else
        rr_pack(trace_directory)
    end

    c = Channel()
    t = @async HTTP.WebSockets.open(WSS_ENDPOINT) do ws
        write(ws, "Hello Server, if it's not too much trouble, please send me S3 credentials")
        x = readavailable(ws)
        put!(c, JSON.parse(String(x))["connectionId"])
        # This will block until the user has completed the authentication flow
        x = readavailable(ws)
        push!(c, JSON.parse(String(x)))
    end
    bind(c, t)
    connectionId = take!(c)
    println()
    println("""
    ### IMPORTANT =============================================================
    You are about to upload a trace directory to a publicly accessible location.
    Such traces contain any information that was accessed by the traced
    executable during its execution. This includes any code loaded, any
    secrets entered, the contents of any configuration files, etc.
    DO NOT proceed, if you do not wish to make this information publicly available.
    By proceeding you explicitly agree to waive any privacy interest in the
    uploaded information.
    ### =======================================================================
    """)
    println("To upload a trace, please authenticate, by visiting:\n")
    println("\thttps://github.com/login/oauth/authorize?client_id=$GITHUB_APP_ID&state=$(HTTP.escapeuri(connectionId))")
    println("\n", "You can cancel upload by `Ctrl-C`.")
    flush(stdout)
    s3creds = try
        take!(c)
    catch err
        if err isa InterruptException
            println()
            println("Cancel uploading the trace.")
            return
        end
        rethrow()
    end

    println()
    @info "Uploading Trace directory"

    creds = AWS.AWSCredentials(
        s3creds["AWS_ACCESS_KEY_ID"],
        s3creds["AWS_SECRET_ACCESS_KEY"],
        s3creds["AWS_SESSION_TOKEN"])
    aws = AWS.AWSConfig(creds = creds, region="us-east-1")

    # Tar it up
    proc = zstdmt() do zstdp
        open(`$zstdp -`, "r+")
    end

    t = @async begin try
        upload = s3_begin_multipart_upload(aws, TRACE_BUCKET, s3creds["UPLOAD_PATH"])
        tags = Vector{String}()
        i = 1
        @Base.Experimental.sync begin
            while isopen(proc)
                buf = Vector{UInt8}(undef, S3_CHUNK_SIZE)
                n = readbytes!(proc, buf)
                n < S3_CHUNK_SIZE && resize!(buf, n)
                resize!(tags, i)
                let partno = i, buf=buf
                    @async begin
                        try
                            tags[partno] = s3_upload_part(aws, upload, partno, buf)
                        catch e
                            close(proc)
                            rethrow(e)
                        end
                    end
                end
                i += 1
            end
        end
        s3_complete_multipart_upload(aws, upload, tags)
    catch e
        Base.showerror(stderr, e)
    end
    end

    # Start the Tar creation process, the file will be uploaded as it's created
    Tar.create(trace_directory, proc)
    close(proc.in)

    wait(t)
    println("Uploaded to https://s3.amazonaws.com/$TRACE_BUCKET/$(s3creds["UPLOAD_PATH"])")
end

function __init__()
    # Read in environment variable settings
    global_record_flags = split(get(ENV, "JULIA_RR_RECORD_ARGS", ""), ' ', keepempty=false)
    ignore_child_status = parse(Bool, get(ENV, "JULIA_RR_IGNORE_STATUS", "false"))
end

end # module