text:
using Compat
using Random
using Ripserer
using StaticArrays
using Test
using TupleTools

using Ripserer:
    _one_hot, _cubemap, _from_cubemap, _to_cubemap, nv, coboundary, boundary, edges, births

include("interfacetest.jl")

@testset "CubeMap" begin
    @testset "cubemap" begin
        @testset "1d" begin
            data = Float64[9, 1, 8, 2, 1, 3]
            @test _cubemap(data) == [9, 9, 1, 8, 8, 8, 2, 2, 1, 3, 3]
        end
        @testset "2d" begin
            data = [1 2; 3 4]
            @test _cubemap(data) == [1 2 2; 3 4 4; 3 4 4]

            data = [1; 2; 3]
            @test _cubemap(data) == [1; 2; 2; 3; 3]
        end
        @testset "nd" begin
            for D in 3:6
                data = ones(ntuple(identity, D)...)
                map = _cubemap(data)

                @test map == ones(ntuple(i -> 2i - 1, D)...)
            end
        end
    end

    @testset "from/to_cubemap" begin
        @test _from_cubemap(CartesianIndex(2, 2, 2), Val(8)) == (
            CartesianIndex(1, 1, 1),
            CartesianIndex(2, 1, 1),
            CartesianIndex(1, 2, 1),
            CartesianIndex(2, 2, 1),
            CartesianIndex(1, 1, 2),
            CartesianIndex(2, 1, 2),
            CartesianIndex(1, 2, 2),
            CartesianIndex(2, 2, 2),
        )

        for vertices in (
            (CartesianIndex(5, 4),),
            (CartesianIndex(1024, 1072), CartesianIndex(1025, 1072)),
            (
                CartesianIndex(1024, 1072, 15),
                CartesianIndex(1025, 1072, 15),
                CartesianIndex(1024, 1073, 15),
                CartesianIndex(1025, 1073, 15),
                CartesianIndex(1024, 1072, 14),
                CartesianIndex(1025, 1072, 14),
                CartesianIndex(1024, 1073, 14),
                CartesianIndex(1025, 1073, 14),
            ),
        )
            @test _from_cubemap(_to_cubemap(vertices), Val(length(vertices))) ==
                TupleTools.sort(vertices)
            root = _to_cubemap(vertices)
            @test begin
                _from_cubemap(root, Val(length(vertices)))
                true
            end
        end

        @test_throws ArgumentError _from_cubemap(CartesianIndex(2, 2, 1), Val(8))
    end
end

@testset "Cube" begin
    @testset "Constructors" begin
        @test Cube{1}([CartesianIndex(1, 1), CartesianIndex(1, 2)], 1) ==
            Cube{1}(CartesianIndex(1, 2), 1)
        @test Cube{1}((CartesianIndex(2, 1), CartesianIndex(1, 1)), 1) ==
            Cube{1}(CartesianIndex(2, 1), 1)
        @test Cube{2}([(1, 1), (1, 2), (2, 2), (2, 1)], 1) ==
            Cube{2}(CartesianIndex(2, 2), 1)
    end

    @testset "Randomized tests for AbstractSimplex interface" begin
        for (D, K) in ((0, 5), (1, 1), (2, 3), (4, 4)), T in (Float32, Int)
            @testset "Cube{$D, $T, $K}" begin
                d = rand(T)
                root = CartesianIndex{K}(ntuple(_ -> 2rand(1:100) - 1, Val(K)))
                dirs = Tuple(sort!(shuffle(1:K)[1:D]))
                for i in dirs
                    root += _one_hot(i, Val(K))
                end
                c = Cube{D}(root, d)

                @testset "Type stuff" begin
                    @test c ≡ Cube{D,T,K}(root, d)
                    @test typeof(c) ≡ Cube{D,T,K}
                    D > 0 && @test_throws DomainError Cube{-D}(root, d)
                    @test eltype(c) == CartesianIndex{K}
                end

                @testset "Signs, getters and equality" begin
                    @test index(c) == root
                    @test c == c
                    @test birth(c) ≡ d
                    @test sign(c) == 1
                    @test -c == c
                    @test abs(c) == c
                    @test dim(c) == D
                end

                @testset "Ordering" begin
                    root2 = CartesianIndex{K}(ntuple(_ -> 2rand(1:100) - 1, Val(K)))
                    @test c < Cube{D}(root2, d + 1)
                    @test c > Cube{D}(root2, d - 1)

                    @test c < Cube{D}(root - 2_one_hot(rand(1:K), Val(K)), d)
                    @test c > Cube{D}(root + 2_one_hot(rand(1:K), Val(K)), d)
                end

                @testset "Vertices, iteration, indexing" begin
                    verts = vertices(c)

                    @test eltype(c) == eltype(verts)
                    @test length(c) == length(verts) == 2^D
                    @test Cube{D}(verts, birth(c)) == c
                    @test tuple(verts...) == verts
                    @test Base.to_index(c) == SVector(verts)

                    for (i, v) in enumerate(c)
                        @test v == verts[i]
                    end

                    @test begin
                        @inferred vertices(c)
                        true
                    end
                end
            end
        end
    end

    @testset "Vertices" begin
        @test vertices(Cube{0}(((1, 1),), 1)) == (CartesianIndex(1, 1),)
        @test vertices(Cube{1}(((1, 1), (1, 2)), 1)) ==
            (CartesianIndex(1, 1), CartesianIndex(1, 2))
        @test vertices(Cube{4}(CartesianIndex(10, 14, 15, 22, 2), 1)) == map(
            CartesianIndex,
            (
                (5, 7, 8, 11, 1),
                (6, 7, 8, 11, 1),
                (5, 8, 8, 11, 1),
                (6, 8, 8, 11, 1),
                (5, 7, 8, 12, 1),
                (6, 7, 8, 12, 1),
                (5, 8, 8, 12, 1),
                (6, 8, 8, 12, 1),
                (5, 7, 8, 11, 2),
                (6, 7, 8, 11, 2),
                (5, 8, 8, 11, 2),
                (6, 8, 8, 11, 2),
                (5, 7, 8, 12, 2),
                (6, 7, 8, 12, 2),
                (5, 8, 8, 12, 2),
                (6, 8, 8, 12, 2),
            ),
        )
    end
end

@testset "Cubical" begin
    for K in 1:4
        data = rand(Float64, ntuple(x -> 10 + x, Val(K)))
        filtration = Cubical(data)

        @test nv(filtration) == length(data)
        @test births(filtration) == data
        @test size(filtration.cubemap) == size(data) .* 2 .- 1
        @test vertices(filtration) == CartesianIndices(data)
    end
end

@testset "Coboundary" begin
    @testset "Edges have no coboundary in 1d" begin
        data = cos.(range(0, 4π; length=1000))

        cob = Cube{2,Float64,1}[]
        flt = Cubical(data)
        cub = simplex(flt, Val(1), (CartesianIndex(3), CartesianIndex(2)))
        for c in coboundary(flt, cub)
            push!(cob, c)
        end
        @test isempty(cob)
    end

    @testset "Coboundary of edges in 2d" begin
        data = [
            2 2 2 0 1 1 1 1 1 1 1
            2 1 2 0 1 2 2 2 2 2 1
            2 1 2 0 1 2 3 3 3 2 1
            2 0 2 0 1 2 3 2 3 2 1
            2 1 2 0 1 2 3 3 3 2 1
            2 1 2 0 1 2 2 2 2 2 1
            2 2 2 0 1 1 1 1 1 1 1
        ]
        cob = []
        flt = Cubical(data)
        cub = Cube{1}([(3, 2), (2, 2)], 1)
        for c in coboundary(flt, cub)
            push!(cob, c)
        end
        @test length(cob) == 2
        @test cob[1] == Cube{2}([(3, 3), (2, 3), (3, 2), (2, 2)], 2)
        @test sign(cob[1]) == 1
        @test cob[2] == Cube{2}(((3, 1), (2, 1), (3, 2), (2, 2)), 2)
        @test sign(cob[2]) == 1
        @test issorted(cob; by=index, rev=true)

        cocob = []
        for c in coboundary(flt, cob[1])
            push!(cocob, c)
        end
        for c in coboundary(flt, cob[2])
            push!(cocob, c)
        end
        @test isempty(cocob)
    end

    @testset "births of 2-cubelet coboundary in 2×2×2 3d data" begin
        data = zeros(2, 2, 2)
        data[:, 1, 1] .= 1
        data[:, 2, 1] .= 2
        data[:, 1, 2] .= 3
        data[:, 2, 2] .= 4

        @testset "Birth increasing, cofacets are supersets of original" begin
            for cube in edges(Cubical(data))
                cob = []
                for c in coboundary(Cubical(data), cube)
                    push!(cob, c)
                end
                @test issorted(cob; by=index, rev=true)
                @test all(c -> issubset(cube, c), cob)
                @test all(c -> birth(cube) ≤ birth(c), cob)
            end
        end
        @testset "All 2-cubes have the same coboundary" begin
            for cube in [
                Cube{2}(((1, 1, 1), (2, 1, 1), (1, 2, 1), (2, 2, 1)), 2.0),
                Cube{2}(((1, 1, 1), (2, 1, 1), (1, 1, 2), (2, 1, 2)), 4.0),
                Cube{2}(((1, 1, 1), (1, 2, 1), (1, 1, 2), (1, 2, 2)), 3.0),
                Cube{2}(((2, 1, 1), (2, 2, 1), (2, 1, 2), (2, 2, 2)), 4.0),
                Cube{2}(((1, 2, 1), (2, 2, 1), (1, 2, 2), (2, 2, 2)), 4.0),
                Cube{2}(((1, 1, 2), (2, 1, 2), (1, 2, 2), (2, 2, 2)), 4.0),
            ]
                cofacet = only(coboundary(Cubical(data), cube))
                @test abs(cofacet) == Cube{3}(
                    (
                        (2, 2, 2),
                        (1, 2, 2),
                        (2, 1, 2),
                        (1, 1, 2),
                        (2, 2, 1),
                        (1, 2, 1),
                        (2, 1, 1),
                        (1, 1, 1),
                    ),
                    4.0,
                )
            end
        end
    end
    @testset "Counting cofacets in 3d" begin
        data = reshape(fill(1.0, 1000), (10, 10, 10))
        cob = []
        flt = Cubical(data)
        cub = Cube{1}(((1, 5, 1), (1, 5, 2)), 1)
        for c in coboundary(flt, cub)
            push!(cob, c)
        end
        @test length(cob) == 3
        @test issorted(cob; by=index, rev=true)

        cocob = []
        for c in coboundary(flt, cob[1])
            push!(cocob, c)
        end
        @test length(cocob) == 1
    end
end

@testset "Boundary" begin
    data = [
        1 2
        3 4
    ]
    cube = Cube{2}(
        (
            CartesianIndex(1, 1),
            CartesianIndex(2, 1),
            CartesianIndex(1, 2),
            CartesianIndex(2, 2),
        ),
        4,
    )

    bnd = []
    for f in boundary(Cubical(data), cube)
        push!(bnd, f)
    end
    @test length(bnd) == 4
    @test issorted(bnd; by=index)
    @test sort(birth.(bnd)) == [2, 3, 4, 4]
end

@testset "ripserer" begin
    @testset "1D" begin
        data = [1, 0, 1, 2, 3, 4, 3, 2, 3, 2, 1, 2]
        d0, d1, d2 = ripserer(Cubical, data; dim_max=2)

        @test d0 == [(2, 3), (1, 4), (0, Inf)]
        @test d1 == []
        @test d2 == []
    end
    @testset "1D representatives" begin
        n = 1000
        x = range(0, 1; length=n)
        curve = sin.(2π * 5x) .* x

        d0, _ = ripserer(Cubical, curve; reps=true)

        for int in d0
            birth_sx = birth_simplex(int)
            @test curve[only(birth_sx)] == birth(int) == birth(birth_sx)
        end
        @test sort!(only.(birth_simplex.(d0))) ==
            CartesianIndex.([1, 157, 354, 552, 752, 951])
    end
    @testset "2D image" begin
        data = [
            0 0 0 0 0
            0 2 2 2 0
            0 2 1 2 0
            0 2 2 2 0
            0 0 0 0 0
        ]

        d0, d1, d2 = ripserer(Cubical, data; reps=true, dim_max=2)

        @test d0 == [(1, 2), (0, Inf)]
        @test d1 == [(0, 2)]
        @test d2 == []

        @test vertices(only(representative(d0[1]))) == (CartesianIndex(3, 3),)
        @test sort(vertices.(representative(d0[2]))) ==
            sort(tuple.(vec(CartesianIndices(data))))
    end
    @testset "3D image" begin
        # Cube with hole in the middle.
        data = zeros(5, 5, 5)
        data[2, 2:4, 2:4] .= 1
        data[3, :, :] .= [0 0 0 0 0; 0 1 1 1 0; 0 1 0 1 0; 0 1 1 1 0; 0 0 0 0 0]
        data[4, 2:4, 2:4] .= 1

        d0, d1, d2 = ripserer(Cubical, data; dim_max=2)

        @test d0 == [(0, 1), (0, Inf)]
        @test d1 == []
        @test d2 == [(0, 1)]
    end
    @testset "Thresholding" begin
        data = [
            1 1 1
            1 2 1
            1 1 1
        ]
        d0, d1 = ripserer(Cubical, data; threshold=1)
        @test d0 == [(1, Inf)]
        @test d1 == [(1, Inf)]
    end
    @testset "Homology and explicit cohomology" begin
        data = zeros(5, 5, 5)
        data[2, 2:4, 2:4] .= 1
        data[3, :, :] .= [0 0 0 0 0; 0 1 1 1 0; 0 1 0 1 0; 0 1 1 1 0; 0 0 0 0 0]
        data[4, 2:4, 2:4] .= 1

        test_filtration(Cubical, data; dim_max=2)
    end
end


text:
using Plots

x = 0:0.1:2π
y = cos.(x)

plot(x, y)

function foo(x, y)
    z = 22(x + y)
    return sin(z)
end

ext: jl


text:
highlighter_exists(config::CthulhuConfig) =
    Sys.which(config.highlighter.exec[1]) !== nothing

__init__() = CONFIG.enable_highlighter = highlighter_exists(CONFIG)

function highlight(io, x, lexer, config::CthulhuConfig)
    _print = endswith(x, '\n') ? print : println
    config.enable_highlighter || return _print(io, x)
    if !highlighter_exists(config)
        @warn "Highlighter command $(config.highlighter.exec[1]) does not exist."
        return _print(io, x)
    end
    cmd = `$(config.highlighter) $lexer`
    open(pipeline(cmd; stdout=io, stderr=stderr), "w") do io
        _print(io, x)
    end
end

function cthulhu_llvm(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig,
                      dump_module = false)
    dump = InteractiveUtils._dump_function_linfo_llvm(
        mi, params.world,
        #=wrapper=# false, #=strip_ir_metadata=# true,
        dump_module,
        optimize, debuginfo ? :source : :none, Base.CodegenParams())
    highlight(io, dump, "llvm", config)
end

function cthulhu_native(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    dump = InteractiveUtils._dump_function_linfo_native(
        mi, params.world,
        #=wrapper=# false, #=syntax=# config.asm_syntax,
        debuginfo ? :source : :none)
    highlight(io, dump, "asm", config)
end

function cthulhu_ast(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    meth = mi.def
    ast = definition(Expr, meth)
    if ast!==nothing
        if !config.pretty_ast
            dump(io, ast; maxdepth=typemax(Int))
        else
            show(io, ast)
            # Meta.show_sexpr(io, ast)
            # Could even highlight the above as some kind-of LISP
        end
    else
        @info "Could not retrieve AST. AST display requires Revise.jl to be loaded." meth
    end
end

function cthulhu_source(io::IO, mi, optimize, debuginfo, params, config::CthulhuConfig)
    meth = mi.def
    src, line = definition(String, meth)
    highlight(io, src, "julia", config)
end

cthulhu_warntype(args...) = cthulhu_warntype(stdout, args...)
function cthulhu_warntype(io::IO, src, rettype, debuginfo, stable_code)
    debuginfo = Base.IRShow.debuginfo(debuginfo)
    lineprinter = Base.IRShow.__debuginfo[debuginfo]
    lambda_io::IOContext = io
    if hasfield(typeof(src), :slotnames) && src.slotnames !== nothing
        slotnames = Base.sourceinfo_slotnames(src)
        lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES => slotnames)
        show_variables(io, src, slotnames)
    end
    print(io, "Body")
    InteractiveUtils.warntype_type_printer(io, rettype, true)
    println(io)
    if isa(src, IRCode)
        show(io, src)
        # XXX this doesn't properly show warntype
    else
        ir_printer = stable_code ? Base.IRShow.show_ir : show_ir
        ir_printer(lambda_io, src, lineprinter(src), InteractiveUtils.warntype_type_printer)
    end
    return nothing
end

function cthulu_typed(io::IO, debuginfo_key, CI, rettype, mi, iswarn, stable_code)
    println(io)
    println(io, "│ ─ $(string(Callsite(-1, MICallInfo(mi, rettype), :invoke)))")

    if iswarn
        cthulhu_warntype(io, CI, rettype, debuginfo_key, stable_code)
    elseif isa(CI, IRCode)
        show(io, CI)
    else
        show(io, CI, debuginfo = debuginfo_key)
    end
    println(io)
end

function show_variables(io, src, slotnames)
    println(io, "Variables")
    slottypes = src.slottypes
    for i = 1:length(slotnames)
        print(io, "  ", slotnames[i])
        if isa(slottypes, Vector{Any})
            InteractiveUtils.warntype_type_printer(io, slottypes[i], true)
        end
        println(io)
    end
    println(io)
end

# These are standard code views that don't need any special handling,
# This namedtuple maps toggle::Symbol to function
const codeviews = (;
    llvm=cthulhu_llvm,
    native=cthulhu_native,
    ast=cthulhu_ast,
    source=cthulhu_source,
)

"""
    Cthulhu.Bookmark

A `Cthulhu.Bookmark` remembers a method marked by `b` key during a descent.
It can be used with the following functions:

* `descend(::Bookmark)`, `descend_code_typed(::Bookmark)`,
  `descend_code_warntype(::Bookmark)`: continue the descent.
* `code_typed(::Bookmark)`, `code_warntype([::IO,] ::Bookmark)`: show typed IR
* `code_llvm([::IO,] ::Bookmark)`: pretty-print LLVM IR
* `code_native([::IO,] ::Bookmark)`: pretty-print native code
"""
struct Bookmark
    mi::MethodInstance
    params::CompilerParams
end

"""
    Cthulhu.BOOKMARKS :: Vector{Bookmark}

During a descent, methods can be "bookmarked" by pressing `b` key.  It
pushes a [`Cthulhu.Bookmark`](@ref) into `Cthulhu.BOOKMARKS`.  This can be
used to, e.g., continue descending by `descend(Cthulhu.BOOKMARKS[end])`.
See [`Cthulhu.Bookmark`](@ref) for other usages.
"""
const BOOKMARKS = Bookmark[]

# Default `show` is broken for `Core.Compiler.Params`. Trying not invoke it.
Base.show(io::IO, b::Bookmark) =
    print(io, "Cthulhu.Bookmark(", b.mi, ", ::", CompilerParams, ")")

# Turn off `optimize` and `debuginfo` for default `show` so that the
# output is smaller.
function Base.show(io::IO, ::MIME"text/plain", b::Bookmark;
                   optimize = false, debuginfo = :none, iswarn=false)
    CI, rt = InteractiveUtils.code_typed(b, optimize = optimize)
    if get(io, :typeinfo, Any) === Bookmark  # a hack to check if in Vector etc.
        print(io, Callsite(-1, MICallInfo(b.mi, rt)), :invoke)
        print(io, " (world: ", b.params.world, ")")
        return
    end
    println(io, "Cthulhu.Bookmark (world: ", b.params.world, ")")
    cthulu_typed(io, debuginfo, CI, rt, b.mi, iswarn)
end

function InteractiveUtils.code_typed(b::Bookmark; optimize = true)
    (CI, rt, slottypes) = do_typeinf_slottypes(b.mi, optimize, b.params)
    preprocess_ci!(CI, b.mi, optimize, CONFIG)
    return CI => rt
end

InteractiveUtils.code_warntype(b::Bookmark; kw...) =
    InteractiveUtils.code_warntype(stdout, b; kw...)
function InteractiveUtils.code_warntype(io::IO, b::Bookmark; debuginfo = :source, kw...)
    CI, rt = InteractiveUtils.code_typed(b; kw...)
    cthulhu_warntype(io, CI, rt, debuginfo)
end

InteractiveUtils.code_llvm(b::Bookmark) = InteractiveUtils.code_llvm(stdout, b)
InteractiveUtils.code_llvm(io::IO, b::Bookmark; optimize = true, debuginfo = :source,
                           dump_module = false, config = CONFIG) =
    cthulhu_llvm(io, b.mi, optimize, debuginfo == :source, b.params, config, dump_module)

InteractiveUtils.code_native(b::Bookmark; kw...) =
    InteractiveUtils.code_native(stdout, b; kw...)
InteractiveUtils.code_native(io::IO, b::Bookmark; optimize = true, debuginfo = :source,
                             config = CONFIG) =
    cthulhu_native(io, b.mi, optimize, debuginfo == :source, b.params, config)

@nospecialize
using Base.IRShow: compute_basic_blocks, scan_ssa_use!, should_print_ssa_type, print_stmt, GotoIfNot, GotoNode, PhiNode, block_for_inst
function show_ir(io::IO, code::Core.CodeInfo, line_info_preprinter, line_info_postprinter)
    cols = displaysize(io)[2]
    used = BitSet()
    stmts = code.code
    types = code.ssavaluetypes
    cfg = compute_basic_blocks(stmts)
    max_bb_idx_size = length(string(length(cfg.blocks)))
    for stmt in stmts
        scan_ssa_use!(push!, used, stmt)
    end
    bb_idx = 1

    if isempty(used)
        maxlength_idx = 0
    else
        maxused = maximum(used)
        maxlength_idx = length(string(maxused))
    end
    for idx in eachindex(stmts)
        if !isassigned(stmts, idx)
            # This is invalid, but do something useful rather
            # than erroring, to make debugging easier
            printstyled(io, "#UNDEF\n", color=:red)
            continue
        end
        stmt = stmts[idx]
        show_type = types isa Vector{Any} && should_print_ssa_type(stmt)
        if types isa Vector{Any} # ignore types for pre-inference code
            if isassigned(types, idx) && show_type
                typ = types[idx]
                if (idx in used) && typ isa Type && (!Base.isdispatchelem(typ) || typ == Core.Box)
                else
                    continue
                end
            else
                continue
            end
        end
        # Compute BB guard rail
        if bb_idx > length(cfg.blocks)
            # If invariants are violated, print a special leader
            linestart = " "^(max_bb_idx_size + 2) # not inside a basic block bracket
            inlining_indent = line_info_preprinter(io, linestart, idx)
            printstyled(io, "!!! ", "─"^max_bb_idx_size, color=:light_black)
        else
            bbrange = cfg.blocks[bb_idx].stmts
            bbrange = bbrange.start:bbrange.stop
            # Print line info update
            linestart = idx == first(bbrange) ? "  " : sprint(io -> printstyled(io, "│ ", color=:light_black), context=io)
            linestart *= " "^max_bb_idx_size
            inlining_indent = line_info_preprinter(io, linestart, idx)
            if idx == first(bbrange)
                bb_idx_str = string(bb_idx)
                bb_pad = max_bb_idx_size - length(bb_idx_str)
                bb_type = length(cfg.blocks[bb_idx].preds) <= 1 ? "─" : "┄"
                printstyled(io, bb_idx_str, " ", bb_type, "─"^bb_pad, color=:light_black)
            elseif idx == last(bbrange) # print separator
                printstyled(io, "└", "─"^(1 + max_bb_idx_size), color=:light_black)
            else
                printstyled(io, "│ ", " "^max_bb_idx_size, color=:light_black)
            end
            if idx == last(bbrange)
                bb_idx += 1
            end
        end
        print(io, inlining_indent, " ")
        # convert statement index to labels, as expected by print_stmt
        if stmt isa Expr
            if stmt.head === :gotoifnot && length(stmt.args) == 2 && stmt.args[2] isa Int
                stmt = GotoIfNot(stmt.args[1], block_for_inst(cfg, stmt.args[2]::Int))
            elseif stmt.head === :enter && length(stmt.args) == 1 && stmt.args[1] isa Int
                stmt = Expr(:enter, block_for_inst(cfg, stmt.args[1]::Int))
            end
        elseif isa(stmt, GotoIfNot)
            stmt = GotoIfNot(stmt.cond, block_for_inst(cfg, stmt.dest))
        elseif stmt isa GotoNode
            stmt = GotoNode(block_for_inst(cfg, stmt.label))
        elseif stmt isa PhiNode
            e = stmt.edges
            stmt = PhiNode(Int32[block_for_inst(cfg, Int(e[i])) for i in 1:length(e)], stmt.values)
        end
        print_stmt(io, idx, stmt, used, maxlength_idx, true, show_type)
        if types isa Vector{Any} # ignore types for pre-inference code
            if !isassigned(types, idx)
                # This is an error, but can happen if passes don't update their type information
                printstyled(io, "::#UNDEF", color=:red)
            elseif show_type
                typ = types[idx]
                if typ isa Union && Base.is_expected_union(typ)
                    Base.emphasize(io, "::$typ", Base.warn_color()) # more mild user notification
                else
                    Base.emphasize(io, "::$typ")
                end
            end
        end
        println(io)
    end
    let linestart = " "^(max_bb_idx_size + 2)
        line_info_preprinter(io, linestart, 0)
    end
    nothing
end
@specialize


text:
# Autogenerated wrapper script for Qt_jll for x86_64-apple-darwin
export libqt53danimation, libqt53dcore, libqt53dextras, libqt53dinput, libqt53dlogic, libqt53dquick, libqt53dquickanimation, libqt53dquickextras, libqt53dquickinput, libqt53dquickrender, libqt53dquickscene2d, libqt53drender, libqt5bodymovin, libqt5charts, libqt5concurrent, libqt5core, libqt5dbus, libqt5gui, libqt5multimedia, libqt5multimediaquick, libqt5multimediawidgets, libqt5network, libqt5opengl, libqt5printsupport, libqt5qml, libqt5qmlmodels, libqt5qmlworkerscript, libqt5quick, libqt5quick3d, libqt5quick3dassetimport, libqt5quick3drender, libqt5quick3druntimerender, libqt5quick3dutils, libqt5quickcontrols2, libqt5quickparticles, libqt5quickshapes, libqt5quicktemplates2, libqt5quicktest, libqt5quickwidgets, libqt5sql, libqt5svg, libqt5test, libqt5widgets, libqt5xml, libqt5xmlpatterns

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt")
JLLWrappers.@declare_library_product(libqt53danimation, "@rpath/Qt3DAnimation")
JLLWrappers.@declare_library_product(libqt53dcore, "@rpath/Qt3DCore")
JLLWrappers.@declare_library_product(libqt53dextras, "@rpath/Qt3DExtras")
JLLWrappers.@declare_library_product(libqt53dinput, "@rpath/Qt3DInput")
JLLWrappers.@declare_library_product(libqt53dlogic, "@rpath/Qt3DLogic")
JLLWrappers.@declare_library_product(libqt53dquick, "@rpath/Qt3DQuick")
JLLWrappers.@declare_library_product(libqt53dquickanimation, "@rpath/Qt3DQuickAnimation")
JLLWrappers.@declare_library_product(libqt53dquickextras, "@rpath/Qt3DQuickExtras")
JLLWrappers.@declare_library_product(libqt53dquickinput, "@rpath/Qt3DQuickInput")
JLLWrappers.@declare_library_product(libqt53dquickrender, "@rpath/Qt3DQuickRender")
JLLWrappers.@declare_library_product(libqt53dquickscene2d, "@rpath/Qt3DQuickScene2D")
JLLWrappers.@declare_library_product(libqt53drender, "@rpath/Qt3DRender")
JLLWrappers.@declare_library_product(libqt5bodymovin, "@rpath/QtBodymovin")
JLLWrappers.@declare_library_product(libqt5charts, "@rpath/QtCharts")
JLLWrappers.@declare_library_product(libqt5concurrent, "@rpath/QtConcurrent")
JLLWrappers.@declare_library_product(libqt5core, "@rpath/QtCore")
JLLWrappers.@declare_library_product(libqt5dbus, "@rpath/QtDBus")
JLLWrappers.@declare_library_product(libqt5gui, "@rpath/QtGui")
JLLWrappers.@declare_library_product(libqt5multimedia, "@rpath/QtMultimedia")
JLLWrappers.@declare_library_product(libqt5multimediaquick, "@rpath/QtMultimediaQuick")
JLLWrappers.@declare_library_product(libqt5multimediawidgets, "@rpath/QtMultimediaWidgets")
JLLWrappers.@declare_library_product(libqt5network, "@rpath/QtNetwork")
JLLWrappers.@declare_library_product(libqt5opengl, "@rpath/QtOpenGL")
JLLWrappers.@declare_library_product(libqt5printsupport, "@rpath/QtPrintSupport")
JLLWrappers.@declare_library_product(libqt5qml, "@rpath/QtQml")
JLLWrappers.@declare_library_product(libqt5qmlmodels, "@rpath/QtQmlModels")
JLLWrappers.@declare_library_product(libqt5qmlworkerscript, "@rpath/QtQmlWorkerScript")
JLLWrappers.@declare_library_product(libqt5quick, "@rpath/QtQuick")
JLLWrappers.@declare_library_product(libqt5quick3d, "@rpath/QtQuick3D")
JLLWrappers.@declare_library_product(libqt5quick3dassetimport, "@rpath/QtQuick3DAssetImport")
JLLWrappers.@declare_library_product(libqt5quick3drender, "@rpath/QtQuick3DRender")
JLLWrappers.@declare_library_product(libqt5quick3druntimerender, "@rpath/QtQuick3DRuntimeRender")
JLLWrappers.@declare_library_product(libqt5quick3dutils, "@rpath/QtQuick3DUtils")
JLLWrappers.@declare_library_product(libqt5quickcontrols2, "@rpath/QtQuickControls2")
JLLWrappers.@declare_library_product(libqt5quickparticles, "@rpath/QtQuickParticles")
JLLWrappers.@declare_library_product(libqt5quickshapes, "@rpath/QtQuickShapes")
JLLWrappers.@declare_library_product(libqt5quicktemplates2, "@rpath/QtQuickTemplates2")
JLLWrappers.@declare_library_product(libqt5quicktest, "@rpath/QtQuickTest")
JLLWrappers.@declare_library_product(libqt5quickwidgets, "@rpath/QtQuickWidgets")
JLLWrappers.@declare_library_product(libqt5sql, "@rpath/QtSql")
JLLWrappers.@declare_library_product(libqt5svg, "@rpath/QtSvg")
JLLWrappers.@declare_library_product(libqt5test, "@rpath/QtTest")
JLLWrappers.@declare_library_product(libqt5widgets, "@rpath/QtWidgets")
JLLWrappers.@declare_library_product(libqt5xml, "@rpath/QtXml")
JLLWrappers.@declare_library_product(libqt5xmlpatterns, "@rpath/QtXmlPatterns")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt53danimation,
        "lib/Qt3DAnimation.framework/Versions/5/Qt3DAnimation",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dcore,
        "lib/Qt3DCore.framework/Versions/5/Qt3DCore",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dextras,
        "lib/Qt3DExtras.framework/Versions/5/Qt3DExtras",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dinput,
        "lib/Qt3DInput.framework/Versions/5/Qt3DInput",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dlogic,
        "lib/Qt3DLogic.framework/Versions/5/Qt3DLogic",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquick,
        "lib/Qt3DQuick.framework/Versions/5/Qt3DQuick",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickanimation,
        "lib/Qt3DQuickAnimation.framework/Versions/5/Qt3DQuickAnimation",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickextras,
        "lib/Qt3DQuickExtras.framework/Versions/5/Qt3DQuickExtras",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickinput,
        "lib/Qt3DQuickInput.framework/Versions/5/Qt3DQuickInput",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickrender,
        "lib/Qt3DQuickRender.framework/Versions/5/Qt3DQuickRender",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53dquickscene2d,
        "lib/Qt3DQuickScene2D.framework/Versions/5/Qt3DQuickScene2D",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt53drender,
        "lib/Qt3DRender.framework/Versions/5/Qt3DRender",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5bodymovin,
        "lib/QtBodymovin.framework/Versions/5/QtBodymovin",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5charts,
        "lib/QtCharts.framework/Versions/5/QtCharts",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5concurrent,
        "lib/QtConcurrent.framework/Versions/5/QtConcurrent",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5core,
        "lib/QtCore.framework/Versions/5/QtCore",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5dbus,
        "lib/QtDBus.framework/Versions/5/QtDBus",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5gui,
        "lib/QtGui.framework/Versions/5/QtGui",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimedia,
        "lib/QtMultimedia.framework/Versions/5/QtMultimedia",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimediaquick,
        "lib/QtMultimediaQuick.framework/Versions/5/QtMultimediaQuick",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5multimediawidgets,
        "lib/QtMultimediaWidgets.framework/Versions/5/QtMultimediaWidgets",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5network,
        "lib/QtNetwork.framework/Versions/5/QtNetwork",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5opengl,
        "lib/QtOpenGL.framework/Versions/5/QtOpenGL",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5printsupport,
        "lib/QtPrintSupport.framework/Versions/5/QtPrintSupport",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qml,
        "lib/QtQml.framework/Versions/5/QtQml",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlmodels,
        "lib/QtQmlModels.framework/Versions/5/QtQmlModels",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlworkerscript,
        "lib/QtQmlWorkerScript.framework/Versions/5/QtQmlWorkerScript",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick,
        "lib/QtQuick.framework/Versions/5/QtQuick",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3d,
        "lib/QtQuick3D.framework/Versions/5/QtQuick3D",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3dassetimport,
        "lib/QtQuick3DAssetImport.framework/Versions/5/QtQuick3DAssetImport",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3drender,
        "lib/QtQuick3DRender.framework/Versions/5/QtQuick3DRender",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3druntimerender,
        "lib/QtQuick3DRuntimeRender.framework/Versions/5/QtQuick3DRuntimeRender",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick3dutils,
        "lib/QtQuick3DUtils.framework/Versions/5/QtQuick3DUtils",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickcontrols2,
        "lib/QtQuickControls2.framework/Versions/5/QtQuickControls2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickparticles,
        "lib/QtQuickParticles.framework/Versions/5/QtQuickParticles",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickshapes,
        "lib/QtQuickShapes.framework/Versions/5/QtQuickShapes",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktemplates2,
        "lib/QtQuickTemplates2.framework/Versions/5/QtQuickTemplates2",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktest,
        "lib/QtQuickTest.framework/Versions/5/QtQuickTest",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickwidgets,
        "lib/QtQuickWidgets.framework/Versions/5/QtQuickWidgets",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5sql,
        "lib/QtSql.framework/Versions/5/QtSql",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5svg,
        "lib/QtSvg.framework/Versions/5/QtSvg",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5test,
        "lib/QtTest.framework/Versions/5/QtTest",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5widgets,
        "lib/QtWidgets.framework/Versions/5/QtWidgets",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xml,
        "lib/QtXml.framework/Versions/5/QtXml",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5xmlpatterns,
        "lib/QtXmlPatterns.framework/Versions/5/QtXmlPatterns",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    AbstractGeoSeries <: DimensionalData.AbstractDimensionalArray

Abstract supertype for high-level `DimensionalArray` that hold `GeoStacks`, `GeoArrays`,
or the paths they can be loaded from. `GeoSeries` are indexed with dimensions
as with a `AbstractGeoArray`. This is useful when you have multiple files containing
rasters or stacks of rasters spread over dimensions like time and elevation.

As much as possible, implementations should facilitate loading entire
directories and detecting the dimensions from metadata.

This allows syntax like for a series of stacks of arrays:

```julia
GeoSeries[Time(Near(DateTime(2001, 1))][:temp][Y(Between(70, 150)), X(Between(-20,20))] |> plot`
```

[`GeoSeries`](@ref) is the concrete implementation.
"""
abstract type AbstractGeoSeries{T,N,D,A} <: AbstractDimensionalArray{T,N,D,A} end

# Interface methods ####################################################

DD.metadata(A::AbstractGeoSeries) = NoMetadata()
DD.name(A::AbstractGeoSeries) = NoName()
DD.label(A::AbstractGeoSeries) = ""

"""
    modify(f, series::AbstractGeoSeries)

Apply function `f` to the data of the child object.
If the child is an `AbstractGeoStack` the function will
be passed on to its child `AbstractGeoArray`s.

`f` must return an idenically sized array.

This method triggers a complete rebuild of all objects,
and disk based objects will be transferred to memory.

This is useful for swapping out array backend for an
entire series to `CuArray` from CUDA.jl to copy data to a GPU,
and potentially other types like `DAarray` from Distributed.jl.
"""
DD.modify(f, A::AbstractGeoSeries) = map(child -> modify(f, child), values(A))

"""
    GeoSeries <: AbstractGeoSeries

    GeoSeries(arrays::AbstractArray{<:AbstractGeoArray}, dims; kw...)
    GeoSeries(stacks::AbstractArray{<:AbstractGeoStack}, dims; kw...)
    GeoSeries(filepaths::AbstractArray{<:AbstractString}, dims; child, duplicate_first, kw...)
    GeoSeries(dirpath:::AbstractString, dims; ext, child, duplicate_first, kw...)

Concrete implementation of [`AbstractGeoSeries`](@ref).
`GeoSeries` are an array of `GeoArray` or `GeoStack`, along some dimension(s).

# Arguments

- `dims`: series dimension/s.

# Keywords

- `refdims`: existing reference dimension/s.
- `child`: constructor of child objects for use with filenames are passed in,
    can be `GeoArray` or `GeoStack`. Defaults to `GeoArray`.
- `duplicate_first::Bool`: wether to duplicate the dimensions and metadata of the
    first file with all other files. This can save load time with a large
    series where dimensions are essentially identical. `true` by default to improve
    load times. If you need exact metadata, set to `false`.
- `ext`: filename extension such as ".tiff" to find when only a directory path is passed in.
- `kw`: keywords passed to the child constructor [`GeoArray`](@ref) or [`GeoStack`](@ref)
    if only file names are passed in.
"""
struct GeoSeries{T,N,D,R,A<:AbstractArray{T,N}} <: AbstractGeoSeries{T,N,D,A}
    data::A
    dims::D
    refdims::R
end
function GeoSeries(data::AbstractArray{<:Union{AbstractGeoStack,AbstractGeoArray}}, dims;
    refdims=()
)
    GeoSeries(data, DD.formatdims(data, dims), refdims)
end
function GeoSeries(filenames::NamedTuple{K}, dims; kw...) where K
    GeoSeries(map((fns...) -> NamedTuple{K}(fns), values(filenames)...), dims; kw...)
end
function GeoSeries(filenames::AbstractArray{<:Union{AbstractString,NamedTuple}}, dims;
    refdims=(), duplicate_first=true, child=nothing, resize=nothing, kw...
)
    childtype = if isnothing(child)
        eltype(filenames) <: NamedTuple ? GeoStack : GeoArray
    else
        child
    end
    data = if duplicate_first
        # We assume all dims, metadata and missingvals are the same over the series
        # We just load the first object, and swap in the filenames of the others.
        data1 = if childtype <: AbstractGeoArray
            childtype(first(filenames); kw...)
        else
            childtype(first(filenames); resize, kw...)
        end
        map(filenames) do fn
            swap_filename(data1, fn)
        end
    else
        # Load everything separately
        if childtype <: AbstractGeoArray
            [childtype(fn; kw...) for fn in filenames]
        else
            [childtype(fn; resize, kw...) for fn in filenames]
        end
    end
    return GeoSeries(data, DD.formatdims(data, dims); refdims)
end
function GeoSeries(dirpath::AbstractString, dims; ext=nothing, kw...)
    filepaths = filter_ext(dirpath, ext)
    GeoSeries(filepaths, dims; kw...)
end

@inline function DD.rebuild(
    A::GeoSeries, data, dims::Tuple, refdims=(), name=nothing, metadata=nothing,
)
    GeoSeries(data, dims, refdims)
end
@inline function DD.rebuild(
    A::GeoSeries;
    data=data(A), dims=dims(A), refdims=refdims(A), name=nothing, metadata=nothing,
)
    GeoSeries(data, dims, refdims)
end

@deprecate series(args...; kw...) GeoSeries(args...; kw...)

# Swap in the filename/s of an object for another filename, wherever it is.
# This is used to use already loaded metadata of one file with another
# file that is similar or identical besides tha actual raster data.
swap_filename(x, filename) = rebuild(x, data=swap_filename(data(x), filename))
swap_filename(x::NamedTuple, filenames::NamedTuple) = map(swap_filename, x, filenames)
swap_filename(x::FileStack, filename::AbstractString) = @set x.filename = filename
swap_filename(x::FileArray, filename::AbstractString) = @set x.filename = filename
function swap_filename(x::AbstractArray, filename::AbstractString)
    # The `FileArray` is wrapped, so use Flatten.jl to update it wherever it is
    ignore = Union{Dict,Set,Base.MultiplicativeInverses.SignedMultiplicativeInverse}
    Flatten.modify(x, FileArray, ignore) do fa
        @set fa.filename = filename
    end
end


text:
using BarFrameworks
using Test

#bf = BarFrameworks.BarFramework(V,E,edgelengths,2)

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


text:
using Condor
using Documenter

DocMeta.setdocmeta!(Condor, :DocTestSetup, :(using Condor); recursive=true)

makedocs(;
    modules=[Condor],
    authors="Yusuke Takase",
    repo="https://github.com/yusuke-takase/Condor.jl/blob/{commit}{path}#{line}",
    sitename="Condor.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://yusuke-takase.github.io/Condor.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/yusuke-takase/Condor.jl",
)


text:
abstract type Microscope end

struct MicroscopeWidefield <: Microscope end
struct MicroscopeConfocal <: Microscope end

max_stars_repo_name: JuliaMicroscopy/Microscopy

max_forks_repo_name: JuliaMicroscopy/Microscopy


text:
function test_variable_degrees()
	for num_vars = 1:10
		for min_degree = 0:3
			for max_degree = 3:5
				degrees = shuffle(collect(min_degree:max_degree))
				for d in Mayday.variable_degrees(num_vars, degrees)
					@test length(d) == num_vars
					@test sum(d) >= min_degree
					@test sum(d) <= max_degree
				end
			end
		end
	end
end


text:
using MusicManipulations
using MIDI
using Test

cd(@__DIR__)

include("midiio.jl")
include("quantizer_tests.jl")
include("timeseries.jl")
include("data_extraction.jl")
include("scale_tests.jl")
include("humanizer_tests.jl")


text:
using RobotDynamics
using Rotations
using Test
using LinearAlgebra
using StaticArrays
using BenchmarkTools
using ForwardDiff

using RobotDynamics: mass, inertia

# Temporary fix for Rotations
# function Rotations.∇rotate(q::UnitQuaternion, r::AbstractVector)
#     Rotations.check_length(r, 3)
#     rhat = UnitQuaternion(zero(eltype(r)), r[1], r[2], r[3], false)
#     R = Rotations.rmult(q)
#     2*Rotations.vmat()*Rotations.rmult(q)'*Rotations.rmult(rhat)
# end

v = @SVector rand(3)
q = rand(UnitQuaternion)
ForwardDiff.jacobian(q->UnitQuaternion(q,false)*v, Rotations.params(q)) ≈ Rotations.∇rotate(q,v)

RD.@autodiff struct Body{R} <: RobotDynamics.RigidBody{R} end
RobotDynamics.control_dim(::Body) = 6

function RobotDynamics.wrenches(model::Body, x::StaticVector, u::StaticVector, t)
    q = orientation(model, x)
    F = q * SA[u[1], u[2], u[3]]  # forces in the body frame
    M = SA[u[4], u[5], u[6]]      # moments in the body frame
    SA[F[1], F[2], F[3], M[1], M[2], M[3]]
end

function RobotDynamics.wrench_jacobian!(F, model::Body, z::RD.AbstractKnotPoint)
    inds = RobotDynamics.gen_inds(model)
    r,q,v,ω = RobotDynamics.parse_state(model, state(z))
    u = control(z)
    i3 = SA[1,2,3]
    f = u[i3]
    F[i3, inds.q] .= Rotations.∇rotate(q, f)
    F[i3, inds.u[i3]] .= RotMatrix(q)
    for i = 4:6
        F[i, inds.u[i]] = 1
    end
    return F
end

function testind(model)
    ir,iq,iv,iω,iu = gen_inds(model)
    @SVector ones(control_dim(model))
end

RobotDynamics.mass(::Body) = 2.0
RobotDynamics.inertia(::Body) = Diagonal(SA[2,3,1.])
RobotDynamics.wrench_sparsity(::Body) = SA[false true  false false true;
                                           false false false false true]

#---
model = Body{UnitQuaternion{Float64}}()
@test RD.dims(model) == (13,6,13)
x,u = rand(model)
q = orientation(model, x)
@test Rotations.params(q) == x[4:7]
@test norm(q) ≈ 1
z = RD.KnotPoint(x,u,0.0,0.01)
@test length(x) == 13
@test length(u) == 6
@test norm(x[4:7]) ≈ 1
@test RobotDynamics.LieState(model) === RobotDynamics.LieState(UnitQuaternion{Float64},(3,6))
@test RobotDynamics.rotation_type(model) == UnitQuaternion{Float64}

# Test initializers
x0,u0 = zeros(model)
@test RBState(model, x0) ≈ zero(RBState)
@test u0 ≈ zeros(RD.control_dim(model))
@test x0[4:7] ≈ [1,0,0,0]

# Test diferent rotations
for R in [UnitQuaternion{Float64}, MRP{Float64}, RodriguesParam{Float64}]
    local model = Body{R}()
	@test RD.state_dim(model) == 9 + Rotations.params(R)
	RobotDynamics.rotation_type(model) == R
	local x0,u0 = zeros(model)
	@test RBState(model, x0) ≈ zero(RBState)
end

# Test gen_inds
inds = RobotDynamics.gen_inds(model)
@test inds.q == 4:7
@test inds.u == 14:19
@test inds.r == 1:3
@test inds.v isa SVector{3,Int}

model2 = Body{MRP{Float64}}()
inds = RobotDynamics.gen_inds(model2)
@test inds.q == 4:6
@test inds.v == 7:9
@test inds.u == 13:18

# Test state building methods
r,q,v,ω = RobotDynamics.parse_state(model, x)
ir,iq,iv,iω,iu = RobotDynamics.gen_inds(model)
@test r == x[ir]
@test Rotations.params(q) == x[iq]
@test v == x[iv]
@test ω == x[iω]

x_ = RBState(r, q, v, ω)
@test RobotDynamics.build_state(model, r, q, v, ω) ≈ x
@test RobotDynamics.build_state(model, r, MRP(q), v, ω) ≈ x
@test RobotDynamics.build_state(model, Vector(r), MRP(q), v, ω) ≈ x
@test RobotDynamics.build_state(model, x_) ≈ x
@test RobotDynamics.build_state(model, SVector(x_)) ≈ x
@test RobotDynamics.build_state(model, Vector(x_)) ≈ x

model2 = Body{MRP{Float64}}()
g = Rotations.params(MRP(q))
x2 = [r; g; v; ω]
@test RobotDynamics.build_state(model2, r, q, v, ω) ≈ x2
@test RobotDynamics.build_state(model2, r, MRP(q), v, ω) ≈ x2
@test RobotDynamics.build_state(model2, r, g, v, ω) ≈ x2

@test RobotDynamics.fill_state(model, 1, 2, 3, 4) isa SVector{13,Int}
@test RobotDynamics.fill_state(model2, 1, 2, 3, 4) isa SVector{12,Int}
@test RobotDynamics.fill_state(model, 1, 2, 3, 4.) isa SVector{13,Float64}
x_ = RobotDynamics.fill_state(model, 1, 2, 3, 4)
@test position(model, x_) ≈ fill(1,3)
@test Rotations.params(orientation(model, x_)) ≈ fill(2,4)
@test linear_velocity(model, x_) ≈ fill(3,3)
@test angular_velocity(model, x_) ≈ fill(4,3)
@test RobotDynamics.fill_state(model2, 1, 2, 3, 4.)[4:6] ≈ fill(2,3)

# Test dynamics
RobotDynamics.velocity_frame(::Body) = :world
xdot = RD.dynamics(model, x, u)
@test xdot ≈ RD.dynamics(model, x, u, 1.0)
@test xdot ≈ RD.dynamics(model, z)
xdot = RBState(xdot)
x_ = RBState(x)
@test position(xdot) ≈ linear_velocity(x_)
@test Rotations.params(orientation(xdot)) ≈ Rotations.kinematics(orientation(x_), angular_velocity(x_))
ξ = RobotDynamics.wrenches(model, z)
F = ξ[SA[1,2,3]]
T = ξ[SA[4,5,6]]
@test linear_velocity(xdot) ≈ F/mass(model)
@test angular_velocity(xdot) ≈
	inertia(model) \ (T - angular_velocity(x_) × (inertia(model) * angular_velocity(x_)))

# Test body-frame velocity
RobotDynamics.velocity_frame(::Body) = :body
xdot = RD.dynamics(model, x, u)
@test xdot ≈ RD.dynamics(model, x, u, 1.0)
@test xdot ≈ RD.dynamics(model, z)
xdot = RBState(xdot)
x_ = RBState(x)
q = orientation(x_)
@test position(xdot) ≈ q * linear_velocity(x_)
@test Rotations.params(orientation(xdot)) ≈ Rotations.kinematics(orientation(x_), angular_velocity(x_))
@test linear_velocity(xdot) ≈ q \ (F / mass(model)) - angular_velocity(x_) × linear_velocity(x_)
@test angular_velocity(xdot) ≈
	inertia(model) \ (T - angular_velocity(x_) × (inertia(model) * angular_velocity(x_)))

# Test flipquat
x = rand(model)[1]
x_ = RobotDynamics.flipquat(model, x)
@test orientation(model, x) ≈ orientation(model, x_)
@test Rotations.params(orientation(model, x)) ≈ -Rotations.params(orientation(model, x_))

# Test RBState methods
x_ = RBState(model, x)
@test position(model, x) ≈ position(x_)
@test orientation(model, x) ≈ orientation(x_)
@test linear_velocity(model, x) ≈ linear_velocity(x_)
@test angular_velocity(model, x) ≈ angular_velocity(x_)
@test RBState(model, x_) isa RBState

# Not Implemented Body
struct FakeBody{R} <: RD.RigidBody{R} end
model = FakeBody{MRP{Float64}}()
@test_throws ErrorException("Not implemented") mass(model)
@test_throws ErrorException("Not implemented") inertia(model)
@test_throws ErrorException("Not implemented") RobotDynamics.forces(model, x, u)
@test_throws ErrorException("Not implemented") RobotDynamics.moments(model, x, u)
@test RobotDynamics.velocity_frame(model) == :world


text:
# Functionalities that are available when Spglib is loaded

using Bravais: reciprocalbasis
using StaticArrays
using LinearAlgebra
import ..KPaths: irrfbz_path

"""
    irrfbz_path(cell::Spglib.Cell)  -->  ::KPath{D}

Returns a **k**-path for an arbitrary 3D unit cell `cell` provided as an `Spglib.Cell`.
The **k**-path is equivalent to the standard **k**-path given by `irrfbz_path` for the
associated standard primitive lattice, but adapted to the possibly non-standard setting of
`cell`.

Note that `cell` must be a primitive unit cell and also cannot be a supercell of a smaller
primitive cell.
"""
function irrfbz_path(cell::Spglib.Cell)
    # extract a standardized primitive basis `pRs` assoc. w/ `cell` via Spglib
    dset = Spglib.get_dataset(cell)
    sgnum = dset.spacegroup_number
    pRs = SVector{3}(SVector{3,Float64}(col) for col in eachcol(dset.std_lattice))

    # if the input cell is a supercell (without any distortion), then the irrfbz algorithm
    # cannot work
    if !isapprox(det(cell.lattice), det(dset.primitive_lattice)) # check volumes agree
        error(DomainError(cell,
            "`cell` is either not primitive or a supercell and is untreatable by " *
            "`irrfbz_path(::Spglib.Cell)`"))
    end

    # Calculate kpath for standard primitive lattice
    kp = irrfbz_path(sgnum, pRs)

    # Now, we convert from the standard primitive lattice to the original lattice.
    # The conversion formula is:
    #     `cell.lattice = rotation * dset.primitive_lattice * transformation`
    # - `rotation` rotates the lattice vectors in 3D space.
    # - `transformation` is an integer-valued matrix that represents the linear combination
    #   of the lattice basis vectors.
    # `transformation` is not used, but is commented here just for information. It can be
    # obtained via
    #   `transformation = inv(primitivebasismatrix(centering(sgnum, 3))) * dset.transformation_matrix'`
    rotation = dset.std_rotation_matrix

    # Rotate k-points in Cartesian space by `rotation`
    pRs_original = SVector{3}(SVector{3,Float64}(col) for col in eachcol(cell.lattice))
    pGs_original = reciprocalbasis(pRs_original)
    cartesianize!(kp)
    for (lab, kv) in points(kp)
        points(kp)[lab] = rotation * kv
    end
    return latticize(kp, pGs_original)
end


text:
# The main workhorse
using Plots, DifferentialEquations, StaticArrays
using JLD2
pyplot()

include("ecme-dox.jl")
u0 = collect(u0_tup)

function isolated_model(u, p, t)
  (d_vm, d_dpsi, d_m_na, d_h_na, d_j_na, d_x_k, d_na_i, d_k_i, d_ca_i, d_ca_jsr, d_ca_nsr, d_ca_ss, d_ltr_ca, d_htr_ca,
     d_po1, d_po2, d_pc2, d_c0, d_c1, d_c2, d_c3, d_c4, d_o, d_cca0, d_cca1, d_cca2, d_cca3, d_y_ca, d_p0, d_p1, d_p2, d_p3, d_n1,
     d_adp_i, d_adp_ic, d_crp_i, d_crp_ic, d_ca_m, d_adp_m, d_nadh, d_isoc, d_akg, d_scoa, d_suc, d_fum, d_mal, d_oaa,
     d_sox_m, d_sox_i, d_h2o2_i, d_gsh_i, d_Q_n, d_Qdot_n, d_QH2_n, d_QH2_p, d_Qdot_p, d_Q_p, d_b1, d_b2, d_b3, d_b4,
     d_fes_ox, d_cytc1_ox, d_cytc_ox) = rhs(u, p, t)

     # Set cytosolic components 0
     d_vm = 0
   	 d_m_na = d_h_na = d_j_na = d_na_i = 0
   	 d_x_k = d_k_i = 0
   	 d_p0 = d_p1 = d_p2 = d_p3 = d_n1 = 0
   	 d_ca_i = d_ca_jsr = d_ca_nsr = d_ca_ss = d_ltr_ca = d_htr_ca = 0
   	 d_adp_i = d_adp_ic = d_crp_i = d_crp_ic = 0
     d_po1 = d_po2 = d_pc2 = d_c0 = d_c1 = d_c2 = d_c3 = d_c4 = d_o = d_cca0 = d_cca1 = d_cca2 = d_cca3 = d_y_ca = 0
     return @SVector[d_vm, d_dpsi, d_m_na, d_h_na, d_j_na, d_x_k, d_na_i, d_k_i, d_ca_i, d_ca_jsr, d_ca_nsr, d_ca_ss, d_ltr_ca, d_htr_ca,
        d_po1, d_po2, d_pc2, d_c0, d_c1, d_c2, d_c3, d_c4, d_o, d_cca0, d_cca1, d_cca2, d_cca3, d_y_ca, d_p0, d_p1, d_p2, d_p3, d_n1,
        d_adp_i, d_adp_ic, d_crp_i, d_crp_ic, d_ca_m, d_adp_m, d_nadh, d_isoc, d_akg, d_scoa, d_suc, d_fum, d_mal, d_oaa,
        d_sox_m, d_sox_i, d_h2o2_i, d_gsh_i, d_Q_n, d_Qdot_n, d_QH2_n, d_QH2_p, d_Qdot_p, d_Q_p, d_b1, d_b2, d_b3, d_b4,
        d_fes_ox, d_cytc1_ox, d_cytc_ox]
end

param = Params(
  pDOX = DOXParams(DOX=0.0),
  pCK = CKParams(V_ATPASE_CYTO= 1E-5),
  pC1 = C1Params(SCALE=1e3 * 60),
  pC3 = C3Params(SCALE=60e3 * 40),
  pC4 = C4Params(SCALE=60e3 * 30),
  pC5 = C5Params(ρF1=0.05),
  pANT = ANTParams(VMAX=5E-3),
  pSODi = SODParams(ET=1.43E-3 * 0.65),
  BCL=0)

# Run 1500 seconds
tspan = (0.0, 1.5E6)
prob = ODEProblem(isolated_model, u0, tspan, param)
baselineSol = solve(prob, Rodas5(); reltol=1e-9, abstol=1e-9, dt=0.01, progress=true, dtmax=3000.0)

paramDOX = Params(
  pDOX = DOXParams(DOX=0.25),
  pCK = CKParams(V_ATPASE_CYTO= 1E-5),
  pC1 = C1Params(SCALE=1e3 * 60),
  pC3 = C3Params(SCALE=60e3 * 40),
  pC4 = C4Params(SCALE=60e3 * 30),
  pC5 = C5Params(ρF1=0.05),
  pANT = ANTParams(VMAX=5E-3),
  pSODi = SODParams(ET=1.43E-3 * 0.65),
  BCL=0)

prob = ODEProblem(isolated_model, u0, tspan, paramDOX)
doxSol = solve(prob, Rodas5(); reltol=1e-9, abstol=1e-9, dt=0.01, progress=true, dtmax=3000.0)

JLD2.@load "mito.jld2" doxSol

pDpsi = plot(doxSol, vars=(0, nameLUT[:dpsi]), label="Mitochondrial Potential", lw=1, legend=:false, xaxis = ("time (ms)"), yaxis=("Voltage (mV)"))
savefig(pDpsi, "dpsi.png")

pATPase = plot(t-> ecme_dox(doxSol(t), paramDOX, t)[rateMap[:vATPase]], doxSol.t[1], doxSol.t[end],
     label="vATPase", ylims=(-0.003, 0.0005), xaxis = ("time (ms)"), yaxis=("ATP synthase rate (mM/ms)"), legend=:false)
savefig(pATPase, "atpase.png")

pROS = plot(doxSol, vars=(0, [nameLUT[:sox_i], nameLUT[:sox_m]]), lw=1, label=["superoxide_cyto" "superoxide_mito"], xaxis = ("time (ms)"), yaxis=("Concentration (mM)"))
savefig(pROS, "ros-burst.png")
plot!(pROS, doxSol, vars=(0, [nameLUT[:Q_n], nameLUT[:QH2_n], nameLUT[:Qdot_p], nameLUT[:Q_p], nameLUT[:QH2_p]]), lw=1, label=["Q_n" "QH2_n" "Qdot_p" "Q_p" "QH2_p"], legend=:right, xaxis = ("time (ms)"))
savefig(pROS, "ros-burst-with-Q.png")

plot()
plot(t-> ecme_dox(doxSol(t), param, t)[rateMap[:vc1]], doxSol.t[1], doxSol.t[end], label="vc1")
plot!(t-> ecme_dox(doxSol(t), param, t)[rateMap[:vROSC1]], doxSol.t[1], doxSol.t[end], label="vROSC1")
plot!(t-> ecme_dox(doxSol(t), param, t)[rateMap[:vROSC3]], doxSol.t[1], doxSol.t[end], label="vROSC3")
ylims!(-1.1e-5, 3e-5)
savefig("vROS.png")

using JLD2
@save "mito.jld2" doxSol

# Other variables you could see
plot(doxSol, vars=(0, [nameLUT[:Q_n], nameLUT[:QH2_n]]), lw=1, label=["Q_n" "QH2_n"])
plot(doxSol, vars=(0, nameLUT[:nadh]), label="[NADH]", lw=1, legend=:bottomleft)
plot(doxSol, vars=(0, nameLUT[:suc]), label="[SUC]", lw=1, legend=:topleft)
plot!(t-> ecme_dox(doxSol(t), paramDOX, t)[rateMap[:vANT]], doxSol.t[1], doxSol.t[end], label="vANT")
plot(doxSol, vars=(0, nameLUT[:adp_m]), label="[ADP]m", lw=1)
plot(doxSol, vars=(0, [nameLUT[:Q_n], nameLUT[:QH2_n]]), lw=1, label=["Q_n" "QH2_n"])
plot(doxSol, vars=(0, [nameLUT[:b1], nameLUT[:b2], nameLUT[:b3], nameLUT[:b4]]), lw=1)
plot(doxSol, vars=(0, [nameLUT[:fes_ox], nameLUT[:cytc1_ox], nameLUT[:cytc_ox]]), lw=1, label=["fes_ox" "cytc1_ox" "cytc_ox"])
plot(doxSol, vars=(0, [nameLUT[:sox_i], nameLUT[:sox_m]]), lw=1, label=["sox_i" "sox_m"])
plot(doxSol, vars=(0, nameLUT[:sox_i]), lw=1, label=["sox_i"])
plot(doxSol, vars=(0, nameLUT[:gsh_i]), lw=1)


text:
using PyPlot
println(1 + 1 + 1)

a = 1 +
1
b = 1
+1

println(a) # a => 2
println(b) # b => 1

plot(1:10)


text:
# Collections of data

struct EdgeGradient{C <: CellType,D <: CellType, NX,NY, T<: Number} <: GridData{NX,NY,T}
  dudx :: Nodes{C,NX,NY,T}
  dvdy :: Nodes{C,NX,NY,T}
  dudy :: Nodes{D,NX,NY,T}
  dvdx :: Nodes{D,NX,NY,T}
end

function EdgeGradient(T::Type{C}, dualnodedims::Tuple{Int, Int};dtype=Float64) where {C <: CellType}
    dudxdims = node_inds(T, dualnodedims)
    dudydims = node_inds(othertype(C), dualnodedims)

    EdgeGradient{T, othertype(T), dualnodedims...,dtype}(
            Nodes(T,dualnodedims,dtype=dtype),Nodes(T,dualnodedims,dtype=dtype),
            Nodes(othertype(T),dualnodedims,dtype=dtype),Nodes(othertype(T),dualnodedims,dtype=dtype)
            )
end

(::Type{EdgeGradient{R,S,NX,NY,T}})() where {R<:CellType,S<:CellType,NX,NY,T} = EdgeGradient(R, (NX, NY),dtype=T)

EdgeGradient(C, ::GridData{NX,NY,T}; dtype=T) where {NX, NY,T} = EdgeGradient(C, (NX,NY),dtype=dtype)

Base.similar(::EdgeGradient{R,S,NX,NY,T};element_type=T) where {R,S,NX,NY,T} = EdgeGradient(R, (NX, NY),dtype=element_type)

function Base.fill!(g::EdgeGradient, s::Number)
    fill!(g.dudx, s)
    fill!(g.dudy, s)
    fill!(g.dvdx, s)
    fill!(g.dvdy, s)
    g
end

Base.size(A::EdgeGradient{C,D,NX,NY,T}) where {C,D,NX,NY,T} = (2length(A.dudx)+2length(A.dudy),1)
@propagate_inbounds Base.getindex(A::EdgeGradient{C,D,NX,NY,T},i::Int) where {C,D,NX,NY,T} =
   i > length(A.dudx) ? (i > length(A.dudx)+length(A.dudy) ? (i > length(A.dudx)+2length(A.dudy) ? A.dvdy[i-length(A.dudx)-2length(A.dudy)] : A.dvdx[i-length(A.dudx)-length(A.dudy)]) : A.dudy[i-length(A.dudx)]) : A.dudx[i]
@propagate_inbounds Base.setindex!(A::EdgeGradient{C,D,NX,NY,T}, v, i::Int) where {C,D,NX,NY,T} =
   i > length(A.dudx) ? (i > length(A.dudx)+length(A.dudy) ? (i > length(A.dudx)+2length(A.dudy) ? A.dvdy[i-length(A.dudx)-2length(A.dudy)] = convert(T, v) : A.dvdx[i-length(A.dudx)-length(A.dudy)] = convert(T, v)) : A.dudy[i-length(A.dudx)]= convert(T, v)) : A.dudx[i] = convert(T, v)
Base.IndexStyle(::Type{<:EdgeGradient}) = IndexLinear()

function Base.show(io::IO, nodes::EdgeGradient{R, S, NX, NY, T}) where {R, S, NX, NY, T}
    nodedims = "(nx = $NX, ny = $NY)"
    udims = "(nx = $(size(nodes.dudx,1)), ny = $(size(nodes.dudx,2)))"
    vdims = "(nx = $(size(nodes.dudy,1)), ny = $(size(nodes.dudy,2)))"
    println(io, "Edge gradient data in a $nodedims cell grid of type $T data")
    println(io, "  du/dx and dv/dy as $R nodes : $udims")
    println(io, "  du/dy and dv/dx as $S nodes : $vdims")
end

struct NodePair{C <: CellType,D <: CellType, NX,NY,T} <: GridData{NX,NY,T}
  u :: Nodes{C,NX,NY,T}
  v :: Nodes{D,NX,NY,T}
end

function NodePair(R::Type{C},S::Type{D},dualnodedims::Tuple{Int,Int};dtype=Float64) where {C <: CellType, D <: CellType}
    udims = node_inds(R, dualnodedims)
    vdims = node_inds(S, dualnodedims)

    NodePair{R, S, dualnodedims...,dtype}(
            Nodes(R,dualnodedims,dtype=dtype),
            Nodes(S,dualnodedims,dtype=dtype)
            )
end
NodePair(R::Type{C}, U::Type{D},nodes::Nodes{S,NX,NY,T}) where
        {C <: CellType, D <: CellType, S <: CellType, NX,NY,T} = NodePair(R, U, (NX, NY),dtype=T)

NodePair(R::Type{C},dualnodedims::Tuple{Int,Int};dtype=Float64) where {C <: CellType}= NodePair(R,R,dualnodedims,dtype=dtype)

NodePair(R::Type{C}, nodes::Nodes{S,NX,NY,T}) where {C <: CellType, S <: CellType, NX,NY,T} = NodePair(R, (NX, NY),dtype=T)

(::Type{NodePair{R,S,NX,NY,T}})() where {R<:CellType,S<:CellType,NX,NY,T} = NodePair(R, S, (NX, NY),dtype=T)

NodePair(R, S, ::GridData{NX,NY,T};dtype=T) where {NX, NY,T} = NodePair(R, S, (NX,NY),dtype=dtype)

Base.similar(::NodePair{R,S,NX,NY,T};element_type=T) where {R<:CellType,S<:CellType,NX,NY,T} = NodePair(R, S, (NX, NY),dtype=element_type)

function Base.show(io::IO, nodes::NodePair{R, S, NX, NY, T}) where {R, S, NX, NY, T}
    nodedims = "(nx = $NX, ny = $NY)"
    udims = "(nx = $(size(nodes.u,1)), ny = $(size(nodes.u,2)))"
    vdims = "(nx = $(size(nodes.v,1)), ny = $(size(nodes.v,2)))"
    println(io, "($R,$S) node data pair in a $nodedims cell grid of type $T data")
    println(io, "  Number of $R nodes: $udims")
    println(io, "  Number of $S nodes: $vdims")
end

function fill!(nodes::NodePair, s::Number)
    fill!(nodes.u, s)
    fill!(nodes.v, s)
    nodes
end

Base.size(A::NodePair{C,D,NX,NY,T}) where {C,D,NX,NY,T} = (length(A.u)+length(A.v),1)
@propagate_inbounds Base.getindex(A::NodePair{C,D,NX,NY,T},i::Int) where {C,D,NX,NY,T} =
   i > length(A.u) ? A.v[i-length(A.u)] : A.u[i]
@propagate_inbounds Base.setindex!(A::NodePair{C,D,NX,NY,T}, v, i::Int) where {C,D,NX,NY,T} =
   i > length(A.u) ? A.v[i-length(A.u)] = convert(T, v) : A.u[i] = convert(T, v)
Base.IndexStyle(::Type{<:NodePair}) = IndexLinear()


text:
module ClenshawCurtisBessel

using LinearAlgebra, BandedMatrices
using SparseArrays, StaticArrays, OffsetArrays
import HypergeometricFunctions
import ArbNumerics, DoubleFloats, SpecialFunctions

include("specialfunctions.jl")  # define i.e. bessel J for high precision types
include("oliver.jl")  # transform recurrence into a boundary value problem
include("piessensbranders.jl")  # compute modified moments

export momentM

end


text:
"Hacked Parameter Space"
module Parameters

export choice, Phi

struct Phi
end

choice(::Phi, xs::AbstractVector) = rand(xs)
choice(::Phi, ::Type{T}) where T <:Number = rand(T)

end


text:
using Test, LinearAlgebra, Random, SparseArrays, ExponentialUtilities
using ExponentialUtilities: getH, getV, _baseexp!, _exp!, ExpMethodNative, ExpMethodDiagonalization, ExpMethodHigham2005,alloc_mem
using ForwardDiff, StaticArrays

@testset "exp!" begin
    n = 100
    A = randn(n, n)
    expA = exp(A);
    methodlist=[ExpMethodNative(),
                ExpMethodDiagonalization(),
                ExpMethodHigham2005()];

    for m in methodlist

        @testset "$(typeof(m))" begin
            E=_exp!(copy(A),m);
            @show E≈expA
            mem=alloc_mem(A,m)
            E=_exp!(copy(A),m);
            @show E≈expA
        end
    end
end

@testset "baseexp" begin
    n = 100
    A = randn(n, n)
    expA = exp(A)
    _baseexp!(A)
    @test A ≈ expA
    A2 = randn(n, n)
    A2 ./= opnorm(A2, 1) # test for small opnorm
    expA2 = exp(A2)
    _baseexp!(A2)
    @test A2 ≈ expA2
end

@testset "Exp" begin
    n = 100
    A = randn(n, n)
    expA = exp(A)
    _exp!(A)
    @test A ≈ expA
    A2 = randn(n, n)
    A2 ./= opnorm(A2, 1) # test for small opnorm
    expA2 = exp(A2)
    _exp!(A2)
    @test A2 ≈ expA2

    # Test all the cases for coverage of all generated code in exp
    A0=[3.0 2.0; 0.0 1.0];
    A0=A0/opnorm(A0,1);
    rhov=[0; 0.015; 0.25; 0.95; 2.1; 5.4];
    for s=1:7
        push!(rhov, rhov[end]*2);
    end
    for (i,_) in enumerate(rhov[1:end-1])
        r=(rhov[i]+rhov[i+1])/2;

        A=A0*r;
        expA=exp(A);
        _exp!(A);
        @test A ≈ expA

        A=A0*r;
        cache=[similar(A) for i=1:5];
        expA=exp(A);
        _exp!(A,caches=cache);
        @test A ≈ expA
    end

end

@testset "exp_generic" begin
    for n in [5, 10, 30, 50, 100, 500]
        M = rand(n, n)
        @test exp(M) ≈ exp_generic(M)

        M′ = M / 10opnorm(M, 1)
        @test exp(M′) ≈ exp_generic(M′)

        N = randn(n, n)
        @test exp(N) ≈ exp_generic(N)

        exp(n) ≈ exp_generic(n)
    end

    @testset "Inf" begin
        @test exp_generic(Inf) == Inf
        @test exp_generic(NaN) === NaN
        @test all(isinf, exp_generic([1 Inf; Inf 1]))
        @test all(isnan, exp_generic([1 Inf; Inf 0]))
        @test all(isnan, exp_generic([1 Inf 1 0; 1 1 1 1; 1 1 1 1; 1 1 1 1]))
    end
end

@testset "Issue 41" begin
    @test ForwardDiff.derivative(exp_generic, 0.1) ≈ exp_generic(0.1) atol=1e-15
end

@testset "Issue 42" begin
    @test exp_generic(0.0) == 1
    @test ForwardDiff.derivative(exp_generic, 0.0) == 1
    @test ForwardDiff.derivative(t -> ForwardDiff.derivative(exp_generic, t), 0.0) == 1
end

@testset "naive_matmul" begin
    A = Matrix(reshape((1.0:23.0^2) ./ 700, (23,23)));
    @test exp_generic(A) ≈ exp(A)
    # FiniteDiff.finite_difference_gradient(sum ∘ exp, A)
    @test ForwardDiff.gradient(sum ∘ exp_generic, A) ≈ [
        1057.6079767302406 1061.5362688188497 1065.4645657136762 1069.39286741472 1073.3211571002205 1077.249468413699 1081.1777677116343 1085.1060742188956 1089.0343494857439 1092.9626752178747 1096.890964903375 1100.8192569919843 1104.747561096137 1108.6758820220502 1112.6041596920072 1116.5324493775076 1120.460763094095 1124.3890623920304 1128.3173592868568 1132.245658584792 1136.173974704488 1140.1022643899885 1144.030546866163;
        1546.2686015296754 1551.8428067925008 1557.4170336833045 1562.9912269305864 1568.5654249840857 1574.1396398593458 1579.7138523314973 1585.2880503849965 1590.8622628571482 1596.4364681199734 1602.0106853983423 1607.584888258059 1613.159081505341 1618.7332915743837 1624.3075064496438 1629.881714115578 1635.455931393947 1641.0301270443374 1646.6043419195976 1652.1785327637708 1657.752740429705 1663.3269745298346 1668.9011653740079;
        2034.9292407477624 2042.1493471692606 2049.3694608000847 2056.5895816402353 2063.8097024803856 2071.02981611121 2078.2499273389253 2085.4700457759673 2092.6901642130088 2099.910270634507 2107.130393877766 2114.350493089938 2121.570616333197 2128.79073236713 2136.010845997954 2143.230966838105 2150.4510684533857 2157.6711844873184 2164.8913077305774 2172.1114261676194 2179.331554217096 2186.5516558323766 2193.771774269418;
        2523.5898631440887 2532.4558707242595 2541.3219119479522 2550.187929140558 2559.0539511393813 2567.9199659288784 2576.78599033081 2585.6520195389594 2594.518051150217 2603.384075552149 2612.2500975509724 2621.116121952904 2629.9821487579447 2638.848151531898 2647.7141975618083 2656.5802195606316 2665.446239156346 2674.3122491396257 2683.178283153992 2692.044300346598 2700.9103199423125 2709.776363569114 2718.642368746176;
        3012.2504783310887 3022.7624255196715 3033.27436309582 3043.786281447098 3054.298197395268 3064.8101421807423 3075.3220869662164 3085.8340173330384 3096.3459476998605 3106.8578804697913 3117.3698108366134 3127.881736397218 3138.393659554714 3148.905589921536 3159.4175154821405 3169.9294602676146 3180.441381022002 3190.9533210012587 3201.4652561742982 3211.9771937504465 3222.489109698616 3233.0010544840907 3243.512987254021;
        3500.911105533632 3513.0689562839966 3525.226797421926 3537.384628947421 3549.542472488459 3561.7003256419325 3573.8581523612097 3586.015995902248 3598.1738226215257 3610.331658953238 3622.489526525363 3634.647350841532 3646.805189576353 3658.9630114894126 3671.120879061538 3683.2787177963587 3695.4365589342883 3707.5943808473485 3719.7522267914956 3731.9100511076645 3744.067906664246 3756.225767027045 3768.3836033587577;
        3989.5717399455016 4003.3754942576475 4017.1792317480326 4030.9829932695047 4044.786757194085 4058.590473056492 4072.3942201593118 4086.197972068349 4100.001714364951 4113.805480692641 4127.609222989243 4141.412986913824 4155.216714791774 4169.02046910392 4182.824213803631 4196.627960906451 4210.431710412379 4224.235455112091 4238.039211827345 4251.8429445115125 4265.646701226768 4279.450453135804 4293.254188223081;
        4478.232371954263 4493.6820346344075 4509.1316828959 4524.581340769827 4540.030991434429 4555.480649308356 4570.930309585392 4586.379962653102 4601.829615720812 4617.279259176087 4632.728936274884 4648.1785749239425 4663.628242410305 4679.077907493558 4694.527553351942 4709.977206419651 4725.426871502905 4740.876519764398 4756.326184847651 4771.775830706035 4787.225483773745 4802.675148856998 4818.124801924709;
        4966.892999156807 4983.988577414276 5001.084124431332 5018.179695479476 5035.275256915185 5052.370806335351 5069.466386995929 5086.56194362542 5103.657507464238 5120.753083318599 5137.848637544982 5154.944198980691 5172.039770028834 5189.135324255217 5206.230890497143 5223.3264591421785 5240.422022980996 5257.517596432248 5274.613138643088 5291.708714497448 5308.804278336266 5325.899844578193 5342.99540841701;
        5455.553635971784 5474.29509135684 5493.036573176091 5511.778040576691 5530.519532008377 5549.26098979654 5568.002459600249 5586.743934210173 5605.485399207664 5624.226888236241 5642.968336411971 5661.70980381257 5680.451290438038 5699.1927674510725 5717.934239657889 5736.675716670922 5755.417176862195 5774.158646665904 5792.900109260285 5811.641574257776 5830.383070495679 5849.124537896278 5867.866010103095;
        5944.214258368111 5964.601636539817 5984.989024323959 6005.376400092557 6025.763761442502 6046.151166048405 6066.53852739835 6086.9259175856005 6107.313293354198 6127.70068113834 6148.088042488285 6168.475437481753 6188.862810847242 6209.250198631384 6229.637569593764 6250.024959781014 6270.412342758938 6290.7997161244275 6311.187099102352 6331.574472467841 6351.961841027113 6372.349233617471 6392.736621401613;
        6432.874883167546 6454.908169707251 6476.941465859391 6498.974752399096 6521.008041341911 6543.0413302847255 6565.074602405779 6587.10789615481 6609.141180291407 6631.17447163733 6653.20776778947 6675.241054329175 6697.2743288533375 6719.307610586826 6741.340925963836 6763.374205294215 6785.407489430812 6807.440792792278 6829.474062510223 6851.50734664682 6873.540640395851 6895.573917323122 6917.607213475262;
        6921.535507966981 6945.214717293336 6968.893909797933 6992.573109511854 7016.252299613341 7039.931506536589 7063.610687025642 7087.289881933346 7110.969084050376 7134.648286167407 7158.327473865786 7182.006666370381 7205.685868487411 7229.365068201333 7253.0442607059285 7276.723453210525 7300.402640908903 7324.081835816607 7347.761042739855 7371.440235244451 7395.119427749047 7418.798627462968 7442.477822370673;
        7410.1961447819585 7435.521238445227 7460.846356139583 7486.171457012178 7511.496562690989 7536.821680385345 7562.1467620330695 7587.471877324317 7612.796975793803 7638.122062247745 7663.447172732774 7688.772283217803 7714.097386493507 7739.422499381645 7764.747607463565 7790.072689111291 7815.397801999428 7840.722898065806 7866.048015760161 7891.37312624519 7916.698215102242 7942.023323184162 7967.348433669192;
        7898.856784000046 7925.8277980468565 7952.798792868798 7979.769828543587 8006.740820962419 8033.711822993687 8060.682839443606 8087.653858296635 8114.624865134119 8141.595881584039 8168.566883615306 8195.537909677661 8222.50890930582 8249.479918546413 8276.450947011876 8303.42195384936 8330.392951074411 8357.363974733657 8384.33498157114 8411.305986005516 8438.277016874088 8465.248026114683 8492.219037758385;
        8387.517403993263 8416.134326408073 8444.751251225991 8473.368171237693 8501.98507683074 8530.601992036225 8559.218914451034 8587.835829656518 8616.452759280653 8645.069684098573 8673.686587288512 8702.303519315758 8730.920446536784 8759.53735693605 8788.154281753968 8816.771192153235 8845.38812177737 8874.005029773529 8902.621940172796 8931.238874603148 8959.855797017957 8988.47271702966 9017.089646653794;
        8876.178026389589 8906.440864381724 8936.70368314899 8966.966497110037 8997.229342311497 9027.492173094306 9057.755011086441 9088.017827450598 9118.28065342719 9148.54347700067 9178.806319799023 9209.069126550745 9239.331950124228 9269.594790519472 9299.857633317823 9330.12043045711 9360.383270852353 9390.646101635162 9420.908927611754 9451.171768006996 9481.43458917737 9511.697424766397 9541.96024593677;
        9364.838670413894 9396.747399952266 9428.656122281312 9460.56486143212 9492.473617404688 9524.382334927519 9556.291074078326 9588.199798810481 9620.108533155071 9652.017279515205 9683.926013859795 9715.834753010602 9747.743477742757 9779.652214490456 9811.56096085059 9843.469690388963 9875.378434345988 9907.28718551234 9939.195893422733 9971.10465179841 10003.013395755435 10034.922110875155 10066.830857235289;
        9853.499285600894 9887.053928313482 9920.60857583229 9954.163220947987 9987.717863660575 10021.272501566948 10054.827149085753 10088.381801410778 10121.93643691404 10155.491101254607 10189.045729548545 10222.600367454916 10256.155007764397 10289.70966008942 10323.264307608226 10356.818945514598 10390.373585824078 10423.928238149103 10457.482880861691 10491.03752117117 10524.592161480652 10558.146825821219 10591.701454115157;
        10342.15992722209 10377.36047109335 10412.561017367721 10447.761582866962 10482.962136350658 10518.162670609485 10553.363216883856 10588.563779979988 10623.764331060574 10658.964884544272 10694.165430818643 10729.365986705447 10764.566544995361 10799.767081657297 10834.967647156536 10870.168195834016 10905.368749317713 10940.569285979649 10975.769858688214 11010.970392947042 11046.170956043174 11081.371499914434 11116.572070219892;
        10830.820552021523 10867.667009067003 10904.513456500046 10941.359930367285 10978.206382606546 11015.052854070675 11051.899306309937 11088.745756146089 11125.592222804002 11162.438677446371 11199.285141701175 11236.131596343545 11272.978065404566 11309.824522450044 11346.670989107955 11383.517436541 11420.363895989587 11457.210350631956 11494.056812483652 11530.903281544674 11567.749738590152 11604.596198038738 11641.442674309084;
        11319.48117441785 11357.973542234437 11396.465905244804 11434.958285076933 11473.450640877976 11511.94301109767 11550.435376511146 11588.927746730842 11627.420116950536 11665.91248717023 11704.404838165055 11742.897220400293 11781.389576201334 11819.881958436572 11858.374326253159 11896.866686860418 11935.359049870785 11973.851427299805 12012.343804728827 12050.836150917434 12089.328530749563 12127.820886550606 12166.313261576517;
        11808.141794411067 11848.280085014305 11888.418346780238 11928.556642189691 11968.694903955624 12008.833187349534 12048.971465937228 12089.109730106267 12129.248003887744 12169.386292087873 12209.524549047586 12249.662842053933 12289.801132657169 12329.939380004449 12370.07766339836 12410.215941986053 12450.354210961312 12490.492489549004 12530.630772942915 12570.769034708848 12610.907310893432 12651.045589481126 12691.183868068818
    ];

    for n ∈ 1:48
        A = rand(n,n);
        B = rand(n,n);
        C = similar(A);
        AB = A*B;
        @test ExponentialUtilities.naivemul!(C, A, B, axes(C,1), axes(C,2), Val(2), Val(1)) ≈ AB
        @test ExponentialUtilities.naivemul!(C, A, B, axes(C,1), axes(C,2), Val(4), Val(2)) ≈ AB
        @test ExponentialUtilities.naivemul!(C, A, B, axes(C,1), axes(C,2), Val(4), Val(3)) ≈ AB
        if n ≤ 16
            Am = MMatrix{n,n}(A)
            Bm = MMatrix{n,n}(B)
            Cm = MMatrix{n,n}(A)
            @test ExponentialUtilities.naivemul!(Cm, Am, Bm, axes(Cm,1), axes(Cm,2), Val(2), Val(2)) ≈ AB
            @test ExponentialUtilities.naivemul!(Cm, Am, Bm, axes(Cm,1), axes(Cm,2), Val(4), Val(2)) ≈ AB
            @test ExponentialUtilities.naivemul!(Cm, Am, Bm, axes(Cm,1), axes(Cm,2), Val(4), Val(3)) ≈ AB
        end
    end
    A = @SMatrix rand(7,7);
    Am = MMatrix(A);
    Aa = Matrix(A);
    @test exp(Aa) ≈ exp_generic(Aa) ≈ exp_generic(Am) ≈ exp_generic(A)
end

@testset "Phi" begin
    # Scalar phi
    K = 4
    z = 0.1
    P = fill(0., K+1); P[1] = exp(z)
    for i = 1:K
        P[i+1] = (P[i] - 1/factorial(i-1))/z
    end
    @test phi(z, K) ≈ P

    # Matrix phi (dense)
    A = [0.1 0.2; 0.3 0.4]
    P = Vector{Matrix{Float64}}(undef, K+1); P[1] = exp(A)
    for i = 1:K
        P[i+1] = (P[i] - 1/factorial(i-1)*I) / A
    end
    @test phi(A, K) ≈ P

    # Matrix phi (Diagonal)
    A = Diagonal([0.1, 0.2, 0.3, 0.4])
    Afull = Matrix(A)
    P = phi(A, K)
    Pfull = phi(Afull, K)
    for i = 1:K+1
        @test Matrix(P[i]) ≈ Pfull[i]
    end
end

@testset "Arnoldi & Krylov" begin
    Random.seed!(0)
    n = 20; m = 5; K = 4
    A = randn(n, n)
    t = 1e-2
    b = randn(n)
    direct = exp(t * A) * b
    @test direct ≈ expv(t, A, b; m=m)
    @test direct ≈ kiops(t, A, b)[1]
    P = phi(t * A, K)
    W = fill(0., n, K+1)
    for i = 1:K+1
        W[:,i] = P[i] * b
    end
    Ks = arnoldi(A, b; m=m)
    W_approx = phiv(t, Ks, K)
    @test W ≈ W_approx
    W_approx_kiops3 = kiops(t, A, hcat([b*inv(t)^i for i in 0:K-1]...))
    @test sum(W[:, 1:K], dims=2) ≈ W_approx_kiops3[1]
    @test_skip begin
        W_approx_kiops4 = kiops(t, A, hcat([b*inv(t)^i for i in 0:K]...))
        @test_broken sum(W[:, 1:K+1], dims=2) ≈ W_approx_kiops4[1]
        @test sum(W[:, 1:K+1], dims=2) ≈ W_approx_kiops4[1] atol=1e-2
    end

    # Happy-breakdown in Krylov
    v = normalize(randn(n))
    A = v * v' # A is Idempotent
    Ks = arnoldi(A, b)
    @test Ks.m == 2

    # Test Arnoldi with zero input
    z = zeros(n)
    Ksz = arnoldi(A, z)
    wz = expv(t, A, z; m=m)
    @test norm(wz) == 0.0

    # Arnoldi vs Lanczos
    A = Hermitian(randn(n, n))
    Aperm = A + 1e-10 * randn(n, n) # no longer Hermitian
    w = expv(t, A, b; m=m)
    wperm = expv(t, Aperm, b; m=m, opnorm=opnorm)
    wkiops = kiops(t, A, b; m=m)[1]
    @test w ≈ wperm
    @test w ≈ wkiops

    # Test Lanczos with zero input
    wz = expv(t, A, z; m=m)
    @test norm(wz) == 0.0
end

@testset "Complex Value" begin
    n = 20; m = 10;
    for A in [Hermitian(rand(ComplexF64, n, n)), Hermitian(rand(n, n)), rand(ComplexF64, n, n), rand(n, n)]
        for b in [rand(ComplexF64, n), rand(n)], t in [1e-2, 1e-2im, 1e-2 + 1e-2im]
            @test exp(t * A) * b ≈ expv(t, A, b; m=m)
        end
    end
end

@testset "Adaptive Krylov" begin
    # Internal time-stepping for Krylov (with adaptation)
    n = 100
    K = 4
    t = 5.0
    tol = 1e-7
    A = spdiagm(-1=>ones(n-1), 0=>-2*ones(n), 1=>ones(n-1))
    B = randn(n, K+1)
    Phi_half = phi(t/2 * A, K)
    Phi = phi(t * A, K)
    uhalf_exact = sum((t/2)^i * Phi_half[i+1] * B[:,i+1] for i = 0:K)
    u_exact = sum(t^i * Phi[i+1] * B[:,i+1] for i = 0:K)
    U = phiv_timestep([t/2, t], A, B; adaptive=true, tol=tol)
    @test norm(U[:,1] - uhalf_exact) / norm(uhalf_exact) < tol
    @test norm(U[:,2] - u_exact) / norm(u_exact) < tol
    # p = 0 special case (expv_timestep)
    u_exact = Phi[1] * B[:, 1]
    u = expv_timestep(t, A, B[:, 1]; adaptive=true, tol=tol, opnorm=opnorm)
    @test_nowarn expv_timestep(t, A, B[:, 1]; adaptive=true, tol=tol, opnorm=opnorm(A, Inf))
    @test norm(u - u_exact) / norm(u_exact) < tol
end

@testset "Krylov for Hermitian matrices" begin
    # Hermitian matrices have real spectra. Ensure that the subspace
    # matrix is representable as a real matrix.

    n = 100
    m = 15
    tol = 1e-14

    e = ones(n)
    p = -im*Tridiagonal(-e[2:end], 0e, e[2:end])

    KsA = KrylovSubspace{ComplexF64}(n, m)
    KsL = KrylovSubspace{ComplexF64, Float64}(n, m)

    v = rand(ComplexF64, n)

    arnoldi!(KsA, p, v)
    lanczos!(KsL, p, v)

    AH = view(KsA.H,1:KsA.m,1:KsA.m)
    LH = view(KsL.H,1:KsL.m,1:KsL.m)

    @test norm(AH-LH)/norm(AH) < tol
end

@testset "Alternative Lanczos expv Interface" begin
    n = 300
    m = 30

    A = Hermitian(rand(n,n))
    b = rand(ComplexF64, n)
    dt = 0.1

    atol=1e-10
    rtol=1e-10
    w = expv(-im, dt*A, b, m=m, tol=atol, rtol=rtol, mode=:error_estimate)

    function fullexp(A, v)
        w = similar(v)
        eA = exp(A)
        mul!(w, eA, v)
        w
    end

    w′ = fullexp(-im*dt*A, b)

    δw = norm(w-w′)
    @test δw < atol
    @test δw/abs(1e-16+norm(w)) < rtol

    z = zeros(ComplexF64, n)
    wz = expv(-im, dt*A, z, m=m, tol=atol, rtol=rtol, mode=:error_estimate)
    @test norm(wz) == 0
end

struct MatrixFreeOperator{T} <: AbstractMatrix{T}
    A::Matrix{T}
end
Base.eltype(A::MatrixFreeOperator{T}) where T = T
LinearAlgebra.mul!(y::AbstractVector, A::MatrixFreeOperator, x::AbstractVector) = mul!(y, A.A, x)
Base.size(A::MatrixFreeOperator, dim) = size(A.A, dim)
struct OpnormFunctor end
(::OpnormFunctor)(A::MatrixFreeOperator, p::Real) = opnorm(A.A, p)
@testset "Matrix-free Operator" begin
    Random.seed!(123)
    n = 20
    for ishermitian in (false, true)
        A = rand(ComplexF64, n, n)
        M = ishermitian ? A'A : A
        Op = MatrixFreeOperator(M)
        b = rand(ComplexF64, n)
        Ks = arnoldi(Op, b; ishermitian=ishermitian, opnorm=OpnormFunctor(), tol=1e-12)
        pv = phiv(0.01, Ks, 2)
        pv′ = hcat(map(A->A*b, phi(0.01Op.A, 2))...)

        @test pv ≈ pv′ atol=1e-12
    end
end


text:
# ---
# title: 70. Climbing Stairs
# id: problem70
# author: zhwang
# date: 2022-01-16
# difficulty: Easy
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/climbing-stairs/description/>
# hidden: true
# ---
#
# You are climbing a staircase. It takes `n` steps to reach the top.
#
# Each time you can either climb `1` or `2` steps. In how many distinct ways can
# you climb to the top?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 2
#     Output: 2
#     Explanation: There are two ways to climb to the top.
#     1. 1 step + 1 step
#     2. 2 steps
#
#
# **Example 2:**
#
#
#
#     Input: n = 3
#     Output: 3
#     Explanation: There are three ways to climb to the top.
#     1. 1 step + 1 step + 1 step
#     2. 1 step + 2 steps
#     3. 2 steps + 1 step
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 45`
#
#
## @lc code=start
using LeetCode

function climbing_stairs(n::Int)::Int
    n == 1 && return 1
    a, b = 1, 2
    for _ in 3:n
        a, b = b, a + b
    end
    return b
end

## @lc code=end


text:
using Bukdu

struct RestController <: ApplicationController
    conn::Conn
end

function init(c::RestController, region::String, site_id::Int, channel_id::Int)
    render(JSON, (:init, region, site_id, channel_id))
end

function update(c::RestController, region::String, site_id::Int, channel_id::Int)
    render(JSON, (:update, region, site_id, channel_id))
end

function build_params(c::RestController)
    region::String = c.params.region
    site_id::Int, channel_id::Int = parse.(Int, (c.params.site_id, c.params.channel_id))
    (c, region, site_id, channel_id)
end

init(c::RestController) = init(build_params(c)...)
update(c::RestController) = update(build_params(c)...)

routes() do
    get("/init/region/:region/site/:site_id/channel/:channel_id/", RestController, init)
    get("/update/region/:region/site/:site_id/channel/:channel_id/", RestController, update)
end

Bukdu.start(8080)

Base.JLOptions().isinteractive == 0 && wait()

max_stars_repo_name: kallooa-bhr/dockerfiles

max_issues_repo_name: kallooa-bhr/dockerfiles

max_forks_repo_name: kallooa-bhr/dockerfiles


text:
# """
# LocalSignaling Module
#
#     Generic methods for prepare and measure scenarios
# """

# """
# num_zero_rows(strategy):
#
#     Counts the number of empty rows in the strategy matrix.
#
# Input:
#     strategy: Matrix, a valid strategy
#
# Output:
#     num_zero_rows: Int, the number of zero rows.
# """
function num_zero_rows(strategy)
    (rows, cols) = size(strategy)

    length(filter( row -> 0 == sum(strategy[row,:]), 1:rows))
end

# """
# vertices_α_interpolable(strategy1, strategy2):
#
#     Checks if two strategies are linearly interpolable over communication protocols.
#
# Input:
#     strategy1/2: Matrix, valid strategies of the same size.
#
# Output:
#     α_interpolable: Boolean, true if a linear interpolation exists on Alice's side
# """
function vertices_α_interpolable(strategy1, strategy2)
    num_zero1 = num_zero_rows(strategy1)
    num_zero2 = num_zero_rows(strategy2)

    main_strat = (num_zero1 <= num_zero2) ? strategy1 : strategy2
    goal_strat = (num_zero1 <= num_zero2) ? strategy2 : strategy1

    (rows, cols) = size(strategy1)
    strat_perms = LocalPolytope.strategies(cols, rows)

    α_interpolable = 0 < length(filter(perm -> main_strat*perm == goal_strat, strat_perms))

    α_interpolable
end

# """
# vertices_β_interpolable(strategy1, strategy2):
#
#     Checks if two strategies are linearly interpolable over decoding strategiess.
#
# Input:
#     strategy1/2: Matrix, valid strategies of the same size.
#
# Output:
#     β_interpolable: Boolean, true if a linear interpolation exists on Bob's side
# """
function vertices_β_interpolable(strategy1, strategy2)
    num_zero1 = num_zero_rows(strategy1)
    num_zero2 = num_zero_rows(strategy2)

    main_strat = (num_zero1 <= num_zero2) ? strategy1 : strategy2
    goal_strat = (num_zero1 <= num_zero2) ? strategy2 : strategy1

    (rows, cols) = size(strategy1)
    strat_perms = LocalPolytope.strategies(cols, rows)

    β_interpolable = 0 < length(filter(perm -> perm * main_strat == goal_strat, strat_perms))

    β_interpolable
end

# """
# vertices_interpolable(strategy1, strategy2):
#
#     Checks if two strategies are linearly interpolable over encodings OR decodings.
#     If the interpolation is over encodings AND decodings then the interpolation
#     is a non-linear polynomial, the function returns false for this case.
#
# Input:
#     strategy1/2: Matrix, valid strategies of the same size.
#
# Output:
#     interpolable: Boolean, true if a linear interpolation exists for either Alice or Bob.
# """
function vertices_interpolable(strategy1, strategy2)
    (vertices_α_interpolable(strategy1, strategy2) || vertices_β_interpolable(strategy1, strategy2))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiAuthorizationV1beta1SelfSubjectAccessReviewSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus } # spec name: status

    function IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview

const _property_map_IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiAuthorizationV1beta1SelfSubjectAccessReviewSpec", Symbol("status")=>"IoK8sApiAuthorizationV1beta1SubjectAccessReviewStatus")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview }) = collect(keys(_property_map_IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview[property_name]

function check_required(o::IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview)
    (getproperty(o, Symbol("spec")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthorizationV1beta1SelfSubjectAccessReview }, name::Symbol, val)
end


text:
using Test
using JuMP
using MIPVerify
@isdefined(TestHelpers) || include("../../TestHelpers.jl")

@testset "pool.jl" begin

    @testset "Pool" begin
        strides = (1, 2, 2, 1)
        p = Pool(strides, MIPVerify.maximum)
        @test p.strides == strides
        @testset "Base.show" begin
            io = IOBuffer()
            Base.show(io, p)
            @test String(take!(io)) == "max pooling with a 2x2 filter and a stride of (2, 2)"
        end
    end

    @testset "getsliceindex" begin
        @testset "inbounds" begin
            @test MIPVerify.getsliceindex(10, 2, 3)==[5, 6]
            @test MIPVerify.getsliceindex(10, 3, 4)==[10]
        end
        @testset "out of bounds" begin
            @test MIPVerify.getsliceindex(10, 5, 4)==[]
        end
    end

    input_size = (6, 6)
    input_array = reshape(1:*(input_size...), input_size)
    @testset "getpoolview" begin
        @testset "inbounds" begin
            @test MIPVerify.getpoolview(input_array, (2, 2), (3, 3)) == [29 35; 30 36]
            @test MIPVerify.getpoolview(input_array, (1, 1), (3, 3)) == reshape([15], (1,1))
        end
        @testset "out of bounds" begin
            @test length(MIPVerify.getpoolview(input_array, (1, 1), (7, 7))) == 0
        end
    end

    @testset "maxpool" begin
        true_output = [
            8 20 32;
            10 22 34;
            12 24 36
        ]
        @testset "Numerical Input" begin
            @test MIPVerify.pool(input_array, MaxPool((2, 2))) == true_output
        end
        @testset "Variable Input" begin
            m = TestHelpers.get_new_model()
            input_array_v = map(
                i -> @variable(m, lowerbound=i-2, upperbound=i),
                input_array
            )
            pool_v = MIPVerify.pool(input_array_v, MaxPool((2, 2)))
            # elements of the input array are made to take their maximum value
            @objective(m, Max, sum(input_array_v))
            solve(m)

            solve_output = getvalue.(pool_v)
            @test solve_output≈true_output
        end
    end

end


text:
"""
    abstract type AbstractFunction <: Function end

An abstract function type.
"""
abstract type AbstractFunction <: Function end

"""
    length(f::AbstractFunction)

Returns the number of outputs of `f`.
"""
Base.length(f::AbstractFunction) = getdim(f)

"""
```
struct FunctionWrapper <: AbstractFunction
    f::Function
    dim::Int
end
```

A function wrapper type that wraps the function `f` where the function `f` is declared to return an output of dimension `dim`.
"""
@params struct FunctionWrapper <: AbstractFunction
    f::Function
    dim::Int
end

"""
    (f::FunctionWrapper)(args...; kwargs...)

Calls the wrapped function in `f` with arguments `args` and keyword arguments `kwargs` and returns the output.
"""
function (f::FunctionWrapper)(args...; kwargs...)
    out = f.f(args...; kwargs...)
    @assert length(out) == getdim(f)
    return out
end

"""
    getfunction(f::FunctionWrapper)

Returns the function wrapped in `f`.
"""
getfunction(f::FunctionWrapper) = f.f

getdim(f::FunctionWrapper) = f.dim

"""
```
struct VectorOfFunctions <: AbstractFunction
    fs
end
```

A struct for a collection of instances of `AbstractFunction`. The dimension of this function is the sum of the dimensions of the individual functions and the outputs are concatenated in a vector.
"""
@params struct VectorOfFunctions <: AbstractFunction
    fs::Union{AbstractVector{<:AbstractFunction}, Tuple{Vararg{AbstractFunction}}}
end

"""
    VectorOfFunctions(fs)

Constructs an instance of `VectorOfFunctions` made of functions `fs`. `fs` can be a vector or tuple of instances of `AbstractFunction`. If a function in `fs` is a `VectorOfFunctions`, it will be unwrapped.
"""
function VectorOfFunctions(fs) end

function VectorOfFunctions(
    fs::Union{
        AbstractVector{>:VectorOfFunctions},
        Tuple{Vararg{>:VectorOfFunctions}},
    },
)
    @assert length(fs) > 0
    new_fs = mapreduce(vcat, fs) do f
        if f isa VectorOfFunctions
            f.fs
        else
            f
        end
    end
    return VectorOfFunctions(new_fs)
end

"""
    (f::VectorOfFunctions)(args...; kwargs...)

Calls the wrapped functions in `f` with arguments `args` and keyword arguments `kwargs` and returns the concatenated output.
"""
function (f::VectorOfFunctions)(args...; kwargs...)
    ys = map(f.fs) do f
        f(args...; kwargs...)
    end
    return vcat(ys...)
end

function getdim(c::VectorOfFunctions)
    if length(c.fs) == 0
        return 0
    else
        return sum(getdim, c.fs)
    end
end

"""
    getfunction(f::VectorOfFunctions, i::Integer)

Returns the `i`^th wrapped function in `f`.
"""
getfunction(f::VectorOfFunctions, i::Integer) = getfunctions(f)[i]

"""
    getfunctions(f::VectorOfFunctions)

Returns the vector or tuple of functions wrapped in `f`.
"""
getfunctions(f::VectorOfFunctions) = f.fs

"""
    abstract type AbstractConstraint <: AbstractFunction end

An abstract constraint type.
"""
abstract type AbstractConstraint <: AbstractFunction end

"""
```
struct Objective <: AbstractFunction
    f
    multiple
end
```

The objective function to be optimized. The objective is always assumed to be of dimension 1. `multiple` is of type `Ref{<:Number}`. When an instance of `Objective` is called, its output will be scaled by a multiplier `multiple[]`. This is 1 by default.
"""
@params struct Objective <: AbstractFunction
    f
    multiple::Base.RefValue
end

"""
    Objective(f, multiple::Number = 1.0)

Constructs an instance of `Objective` wrapping `f`. When an instance of `Objective`, it calls `f` returning the output multiplied by `multiple`. `f` must return a number.
"""
Objective(f, multiple::Number = 1.0) = Objective(f, Ref(multiple))

"""
    (obj::Objective)(args...; kwargs...)

Calls the wrapped function in `obj` with arguments `args` and keyword arguments `kwargs` returning the output multiplied by `obj.multiple[]`. The output of the wrapped function must be a number.
"""
function (o::Objective)(args...; kwargs...)
    out = o.f(args...; kwargs...) * o.multiple[]
    @assert out isa Number
    return out
end

getdim(c::Objective) = 1

"""
    getfunction(obj::Objective)

Returns the function wrapped in `obj`.
"""
getfunction(obj::Objective) = obj.f

"""
```
struct IneqConstraint <: AbstractConstraint
    f::Function
    rhs
    dim::Int
end
```

A struct for an inequality constraint of the form `f(x) .- rhs .<= 0`. The dimension of the constraint is `dim`. Calling the struct will return `f(x) .- rhs`.
"""
@params struct IneqConstraint <: AbstractConstraint
    f::Function
    rhs
    dim::Int
end

"""
    IneqConstraint(f, rhs)

Constructs an instance of `IneqConstraint` with a left-hand-side function `f` and right-hand-side bound `rhs`. If `f` is an instance of `AbstractFunction`, the dimension of the constraint will be `getdim(f)`, otherwise it is assumed to be `length(rhs)` by default.
"""
IneqConstraint(f, rhs) = IneqConstraint(f, rhs, length(rhs))
IneqConstraint(f::AbstractFunction, rhs) = IneqConstraint(f, rhs, getdim(f))

"""
    (c::IneqConstraint)(args...; kwargs...)

Calls the wrapped function in the constraint `c`, `c.f`, with arguments `args` and keyword arguments `kwargs` returning the constraint violation `c.f(args...; kwargs...) .- c.rhs`.
"""
function (c::IneqConstraint)(args...; kwargs...)
    out = c.f(args...; kwargs...) .- c.rhs
    @assert length(out) == getdim(c)
    return out
end

getdim(c::IneqConstraint) = c.dim

"""
    getfunction(f::IneqConstraint)

Returns the function wrapped in `f`.
"""
getfunction(f::IneqConstraint) = f.f

"""
```
struct EqConstraint <: AbstractConstraint
    f::Function
    rhs
    dim::Int
end
```

A struct for an equality constraint of the form `f(x) .- rhs .== 0`. The dimension of the constraint is `dim`. Calling the struct will return `f(x) .- rhs`.
"""
@params struct EqConstraint <: AbstractConstraint
    f::Function
    rhs
    dim::Int
end

"""
    EqConstraint(f, rhs)

Constructs an instance of `EqConstraint` with a left-hand-side function `f` and right-hand-side bound `rhs`. If `f` is an instance of `AbstractFunction`, the dimension of the constraint will be `getdim(f)`, otherwise it is assumed to be `length(rhs)` by default.
"""
EqConstraint(f, rhs) = EqConstraint(f, rhs, length(rhs))
EqConstraint(f::AbstractFunction, rhs) = EqConstraint(f, rhs, getdim(f))

"""
    (c::EqConstraint)(args...; kwargs...)

Calls the wrapped function in the constraint `c`, `c.f`, with arguments `args` and keyword arguments `kwargs` returning the constraint violation `c.f(args...; kwargs...) .- c.rhs`.
"""
function (c::EqConstraint)(args...; kwargs...)
    out = c.f(args...; kwargs...) .- c.rhs
    @assert length(out) == getdim(c)
    return out
end

getdim(c::EqConstraint) = c.dim

"""
    getfunction(f::EqConstraint)

Returns the function wrapped in `f`.
"""
getfunction(f::EqConstraint) = f.f