text:
@testset "Persistence" begin

    # Intervals
    p = Interval(1.0 => 2.0)
    q = Interval(0.0 => 3.0)
    r = Interval(0.0 => 4.0)

    @test q < p
    @test q < r

    @test ComputationalHomology.birth(p) == -1.
    @test ComputationalHomology.death(p) == 3.
    @test diag(p)  == Interval(1.5 => 1.5)

    #######
    # from "Computational Topology - An Introduction" by Edelsbrunner & Harer, p. 184
    #######
    flt = Filtration(SimplicialComplex{Simplex{Int}}, Int)
    splxs=[ Simplex(1) => 1,
            Simplex(2) => 2,
            Simplex(3) => 3,
            Simplex(1,2) => 4,
            Simplex(2,3) => 5,
            Simplex(3,1) => 6,
            Simplex(1,2,3) => 7]
    for s in splxs
        push!(flt, s...)
    end
    @test size(complex(flt)) == (3,3,1)

    # Twist reduction
    ps, R = ComputationalHomology.pairs(TwistReduction, flt, reduced=true)
    @testset for (p, t) = zip(ps, [0=>1, 2=>4, 3=>5, 6=>7])
        @test p == t
    end

    # Standard reduction
    ps, R = ComputationalHomology.pairs(StandardReduction, flt)
    @testset for (p, t) = zip(ps, [2=>4, 3=>5, 6=>7, 1=>Inf])
        @test p == t
    end

    #######
    # from "Topology for Computing" by Zomorodian, pp.138-145
    #######
    flt = Filtration(SimplicialComplex{Simplex{Char}}, Int)
    splxs=[ Simplex('a') => 0,
            Simplex('b') => 0,
            Simplex('c') => 1,
            Simplex('d') => 1,
            Simplex('a', 'b') => 1,
            Simplex('b', 'c') => 1,
            Simplex('c', 'd') => 2,
            Simplex('a', 'd') => 2,
            Simplex('a', 'c') => 3,
            Simplex('a', 'b', 'c') => 4,
            Simplex('a', 'c', 'd') => 5]
    for s in splxs
        push!(flt, s...)
    end

    @testset "Intervals " for (d, itrs) in intervals(flt, length0=true)
        titrs = d == 0 ? intervals(0, 0=>Inf, 0=>1, 1=>1, 1=>2) : intervals(1, 3=>4, 2=>5)
        for itr in itrs
            @test itr ∈ titrs
        end
    end

    #######
    # Ex.1
    #######
    flt = Filtration(SimplicialComplex{Simplex{Int}}, Int)
    splxs=[ Simplex(1) => 0,
            Simplex(2) => 0,
            Simplex(3) => 0,
            Simplex(4) => 0,
            Simplex(5) => 1,
            Simplex(1,2) => 0,
            Simplex(2,3) => 0,
            Simplex(3,4) => 0,
            Simplex(4,1) => 0,
            Simplex(3,5) => 2,
            Simplex(4,5) => 3,
            Simplex(3,4,5) => 7]
    for s in splxs
        push!(flt, s...)
    end

    itr = intervals(flt, length0=true)
    @test Interval(0=>0) ∈ itr[0]
    @test Interval(1=>2) ∈ itr[0]
    @test Interval(1, 3=>7) ∈ itr[1]

    itr = intervals(flt, reduction = StandardReduction)
    @test Interval(1=>2) ∈ itr[0]
    @test Interval(1, 3=>7) ∈ itr[1]
    @test Interval(1, 0=>Inf) ∈ itr[1]

    #######
    # Ex.2
    #######
    flt = Filtration(SimplicialComplex{Simplex{Int}}, Float64)
    splxs=[ Simplex(1) => 1.0,
            Simplex(2) => 2.0,
            Simplex(3) => 3.0,
            Simplex(4) => 4.0,
            Simplex(5) => 5.0,
            Simplex(1,2) => 6.0,
            Simplex(2,3) => 7.0,
            Simplex(3,4) => 8.0,
            Simplex(4,1) => 9.0,
            Simplex(3,5) => 10.0,
            Simplex(4,5) => 11.0,
            Simplex(3,4,5) => 12.0]
    for s in splxs
        push!(flt, s...)
    end
    @test size(complex(flt)) == (5,6,1)
    @test length(flt) == 12

    # Betti numbers
    ∂ = boundary(flt)
    @test sparse(∂)[5,11] == 5
    @test length(∂[9]) > 0
    R = reduce(StandardReduction, ∂)
    @test ComputationalHomology.betti(flt, R, 0) == 1
    @test ComputationalHomology.betti(flt, R, 1) == 1
    @test ComputationalHomology.betti(flt, R, 2) == 0
    @test_throws AssertionError ComputationalHomology.betti(flt, R, 3)
    @test length(R[9]) == 0
    reduce!(StandardReduction, ∂)
    @test length(∂[9]) == 0

    # PH object
    ph = persistenthomology(flt)
    @test eltype(ph) == Tuple{Int64,Int64}
    @test length(ph) == 3
    @test group(ph, 0) == 1
    @test group(ph, 1) == 1
    @test group(ph, 2) == 0
    @test_throws AssertionError group(ph, 3)

    # PH iterator
    @testset "Method Comparison" for (g1, g2) in zip(homology(complex(flt), Int), persistenthomology(TwistReduction, flt))
        @test g1[1] == g2[1]
        @test g1[2] == g2[2]
    end
end


text:
@testset "Interpolate" begin
    @testset "Default" begin
        tester = ImputorTester(Interpolate)

        test_hashing(tester)
        test_equality(tester)
        test_vector(tester)
        test_matrix(tester)
        # test_cube(tester)
        test_dataframe(tester)
        test_groupby(tester)
        test_axisarray(tester)
        test_nameddimsarray(tester)
        test_keyedarray(tester)
        test_columntable(tester)
        test_rowtable(tester)

        @testset "Cube" begin
            a = allowmissing(1.0:1.0:60.0)
            a[[2, 7, 18, 23, 34, 41, 55, 59, 60]] .= missing
            C = collect(reshape(a, 5, 4, 3))

            # Cube tests are expected to fail
            @test_throws MethodError impute(C, tester.imp(; tester.kwargs...); dims=3)
        end
    end

    @testset "Floats" begin
        # Defining our missing datasets
        a = allowmissing(1.0:1.0:20.0)
        a[[2, 3, 7]] .= missing

        result = impute(a, Interpolate())
        @test result == collect(1.0:1.0:20)
        @test result == Impute.interp(a)

        # Test in-place method
        a2 = copy(a)
        Impute.interp!(a2)
        @test a2 == result

        # Test interpolation between identical points
        b = ones(Union{Float64, Missing}, 20)
        b[[2, 3, 7]] .= missing
        @test Impute.interp(b) == ones(Union{Float64, Missing}, 20)

        # Test interpolation at endpoints
        b = ones(Union{Float64, Missing}, 20)
        b[[1, 3, 20]] .= missing
        result = Impute.interp(b)
        @test ismissing(result[1])
        @test ismissing(result[20])

        # Test limiting
        c = allowmissing(1.0:1.0:20.0)
        c[13:15] .= missing

         # Limit too small for gap
        expected = copy(c)
        @test isequal(impute(c, Interpolate(; limit=2)), expected)

        # Limit matches gap size
        expected[13:15] .= [13.0, 14.0, 15.0]
        @test isequal(impute(c, Interpolate(; limit=3)), expected)
    end

    @testset "Ints" begin
        # Defining our missing datasets
        a = allowmissing(1:1:20)
        a[[2, 3, 7]] .= missing

        result = impute(a, Interpolate())
        @test result == collect(1:1:20)
        @test result == Impute.interp(a)

        # Test in-place method
        a2 = copy(a)
        Impute.interp!(a2)
        @test a2 == result

        # Test interpolation between identical points
        b = ones(Union{Float64, Missing}, 20)
        b[[2, 3, 7]] .= missing
        @test Impute.interp(b) == ones(Union{Float64, Missing}, 20)

        # Test interpolation at endpoints
        b = ones(Union{Float64, Missing}, 20)
        b[[1, 3, 20]] .= missing
        result = Impute.interp(b)
        @test ismissing(result[1])
        @test ismissing(result[20])

        # Test inexact error
        # https://github.com/invenia/Impute.jl/issues/71
        c = [1, missing, 2, 3]
        @test_throws InexactError Impute.interp(c)
    end

    # TODO Test error cases on non-numeric types
end

ext: jl


text:
list_test_str = split(
"forward 5
down 5
forward 8
up 3
down 8
forward 2",
'\n'
)

list_str = readlines("day2_input.txt")

# Part 1

function part1(list)
    depth = 0
    pos = 0
    for cmd in list
        m, k = split(cmd)
        k = parse(Int, k)
        pos += m == "forward" ? k : 0
        depth += m == "up" ? -k : m == "down" ? k : 0
    end
    return pos * depth
end

@info "Part 1"
@info "test:"
@show part1(list_test_str) == 150
@info "challenge:"
@show part1(list_str) == 1714680

# Part 2

function part2(list)
    depth = 0
    pos = 0
    aim = 0
    for cmd in list
        m, k = split(cmd)
        k = parse(Int, k)
        pos += m == "forward" ? k : 0
        depth += m == "forward" ? aim * k : 0
        aim += m == "up" ? -k : m == "down" ? k : 0
    end
    return pos * depth
end

@info "Part 2"
@info "test:"
@show part2(list_test_str) == 900
@info "challenge:"
@show part2(list_str) == 1963088820

# Alternative

# Part 1

forward(k) = [0, k] # depth and pos
up(k) = [-k, 0]
down(k) = [k, 0]

function part1_alt(list)
    commands = [split(cmd) |> u -> (Symbol(first(u)), parse(Int, last(u))) for cmd in list]
    return prod(mapreduce(u -> eval(first(u))(last(u)), +, commands))
end

@info "Part 1 alt"
@info "test:"
@show part1_alt(list_test_str) == 150
@info "challenge:"
@show part1_alt(list_str) == 1714680

# Part 2

function forward!(u, k)
    u[1] += u[3] * k
    u[2] += k
end
function up!(u, k)
    u[3] -= k
end
function down!(u, k)
    u[3] += k
end

function part2_alt2(list)
    u = [0, 0, 0]
    for l in list
        cmd = l[1]
        k = parse(Int, l[end])
        if cmd == 'f'
            forward!(u, k)
        elseif cmd == 'u'
            up!(u, k)
        else
            down!(u, k)
        end
    end
    return u[1] * u[2]
end

@info "Part 2 alt"
@info "test:"
@show part2_alt2(list_test_str) == 900
@info "challenge:"
@show part2_alt2(list_str) == 1963088820

nothing


text:
# # includes -
# include("Include.jl")
#
# # setup some global consts -
# const path_to_cobra_mat_file = "$(pwd())/cobra/config/matlab_cobra_files/modelReg.mat"
# const model_file_name = "modelReg"
# const organism_id = :eco    # we use the KEGG organism symbols
# const path_to_measurements_file = "$(pwd())/experimental_data/test_data/Glucose.json"

function moma_calculation(flux_soln_array, flux_bounds_array, ko_index)

    # load the default data_dictionary -
    default_data_dictionary = generate_default_data_dictionary(path_to_cobra_mat_file, model_file_name, organism_id);

    # pass the default dictionary to a customization method -
    updated_data_dictionary = optimize_specific_growth_rate(default_data_dictionary);

    # KO -
    flux_bounds_array[ko_index,1] = 0.0
    flux_bounds_array[ko_index,2] = 0.0
    updated_data_dictionary["flux_bounds_array"] = flux_bounds_array

    # execute MOMA calculation -
    moma_soln = moma(flux_soln_array, updated_data_dictionary)

    # return -
    return moma_soln
end

# # call an executable function -
# (flux_distribution, results_array) = maximize_specific_growth_rate(path_to_measurements_file; number_of_samples = 1000);
#
# # take the mean -
# wild_type_soln = results_array[10].flux_array
# wild_type_soln_fba = results_array[10].flux_bounds_array
# moma_soln = moma_calculation(wild_type_soln, wild_type_soln_fba, 12)


text:
module JLAssistant

    import TOML
    import LocalRegistry
    import ArgParse
    import PkgTemplates
    import LibGit2
    import Pkg: Pkg, PackageSpec
    import RegexTools: hex_escape

    using FilesTreeTools

    include("_find_unregistered_pkgs.jl")
    include("_project_toml_utils.jl")
    include("_walk_pkgs.jl")
    include("create_pkg_version.jl")
    include("commit_to_registry.jl")
    include("precompile_projects.jl")
    include("generate_pkg.jl")
    include("utils.jl")
    include("import_pkgs.jl")
    include("update_manifests.jl")
    include("find_imports.jl")
    include("check_imports.jl")

end


text:
export Times
"""
    `Times(dt, tfin; ns, Ns)`

Constructs a mesh of times.

# Arguments
- `dt`  : timestep;
- `tfin`: final time;
- `ns`  : data are stored every `ns` computations (optional, default = 1).
- `Ns`  : `Ns` data (in addition to the initial datum) are stored (optional, by default `floor( tfin/dt)).
If both `Ns` and `ns` are given, `Ns` overrules `ns`.

# Return values
`t=Times(args)` is of parametric type and offers
- `t.Nc`: number of computed times;
- `t.Ns`: number of stored times;
- `t.ns`: number of computed times between two stored times;
- `t.tfin`: the final time;
- `t.dt`: the timestep;
- `t.tc` : the vector of computed times;
- `t.ts`: the vector of stored times.

"""
struct Times

    Nc   :: Int
    Ns   :: Int
    ns   :: Vector{Int}
    tfin :: Float64
    dt   :: Float64
    tc    :: Vector{Float64}
    ts   :: Vector{Float64}

    function Times( dt, tfin; ns = 1, Ns=nothing)
        tc = range(0, stop=tfin, step = dt)
        Nc = length(tc)
        if Ns != nothing
            ind_stored = round.(Int,range(1, stop=Nc, length = minimum([Ns+1 Nc])))
        else
            ind_stored = round.(Int,range(1, stop=Nc, step = ns))
        end
        ts = tc[ind_stored]
        Ns = length(ts)
        ns = ind_stored[2:end]-ind_stored[1:end-1]
        new( Nc, Ns, ns, tfin, dt, tc, ts)
    end

end


text:
compiler_env, script, output = ARGS

if VERSION < v"0.7-"
    error("Unsupported Julia version $VERSION")
end

const Pkg =
    Base.require(Base.PkgId(Base.UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"), "Pkg"))

Pkg.activate(compiler_env)
@info "Loading PackageCompiler..."
using PackageCompiler

@info "Installing PyCall..."
Pkg.activate(".")
Pkg.add([
    Pkg.PackageSpec(
        name = "PyCall",
        uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0",
    ),
])

@info "Compiling system image..."
create_sysimage(
    [:PyCall],
    sysimage_path = output,
    project = ".",
    precompile_execution_file = script,
)

@info "System image is created at $output"

max_stars_repo_name: grahamgill/pyjulia

max_issues_repo_name: grahamgill/pyjulia

max_forks_repo_name: grahamgill/pyjulia


text:
@testset "Basic Type coercion tests" begin
    X = (x=10:10:44, y=1:4, z=collect("abcd"))

    @test_throws ScientificTypes.CoercionError coerce(X, :x=>Float64)
    @test_throws ScientificTypes.CoercionError coerce(X, :x=>Textual)

    types = Dict(:x => Continuous, :z => Multiclass)
    X_coerced = coerce(X, types)
    @test X_coerced ==  coerce(X, :x => Continuous, :z => Multiclass)
    @test scitype_union(X_coerced.x) === Continuous
    @test scitype_union(X_coerced.z) <: Multiclass
    @test !X_coerced.z.pool.ordered
    @test_throws MethodError coerce(["a", "b", "c"], Count)

    y = collect(Float64, 1:5)
    y_coerced = coerce(y, Count)
    @test scitype_union(y_coerced) === Count
    @test y_coerced == y
    y = [1//2, 3//4, 6//5]
    y_coerced = coerce(y, Continuous)
    @test scitype_union(y_coerced) === Continuous
    @test y_coerced ≈ y
    X_coerced = coerce(X, Dict(:z => OrderedFactor))
    @test X_coerced == coerce(X, :z => OrderedFactor)
    @test X_coerced.x === X.x
    @test scitype_union(X_coerced.z) <: OrderedFactor
    @test X_coerced.z.pool.ordered

    # Check no-op coercion
    y = rand(Float64, 5)
    @test coerce(y, Continuous) === y
    y = rand(Float32, 5)
    @test coerce(y, Continuous) === y
    y = rand(BigFloat, 5)
    @test coerce(y, Continuous) === y
    y = rand(Int, 5)
    @test coerce(y, Count) === y
    y = big.(y)
    @test coerce(y, Count) === y
    y = rand(UInt32, 5)
    @test coerce(y, Count) === y
    X_coerced = coerce(X, Dict{Symbol, Type}())
    @test X_coerced.x === X.x
    @test X_coerced.z === X.z
    z = categorical(X.z)
    @test coerce(z, Multiclass) === z
    z = categorical(X.z, compress=true, ordered=false)
    @test coerce(z, Multiclass) === z
    z = categorical(X.z, compress=true, ordered=true)
    @test coerce(z, OrderedFactor) === z

    # missing values
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce([4, 7, missing], Continuous))
    @test ismissing(y_coerced == [4.0, 7.0, missing])
    @test scitype_union(y_coerced) === Union{Missing,Continuous}
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Any"),
        coerce(Any[4, 7.0, missing], Continuous))
    @test ismissing(y_coerced == [4.0, 7.0, missing])
    @test scitype_union(y_coerced) === Union{Missing,Continuous}
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce([4.0, 7.0, missing], Count))
    @test ismissing(y_coerced == [4, 7, missing])
    @test scitype_union(y_coerced) === Union{Missing,Count}
    y_coerced = @test_logs(
        (:info, r"Trying to coerce from `Any"),
        coerce(Any[4, 7.0, missing], Count))
    @test ismissing(y_coerced == [4, 7, missing])
    @test scitype_union(y_coerced) === Union{Missing,Count}
    @test scitype_union(@test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce(['x', 'y', missing], Multiclass))) ===
            Union{Missing, Multiclass{2}}
    @test scitype_union(@test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce(['x', 'y', missing], OrderedFactor))) ===
            Union{Missing, OrderedFactor{2}}
    # non-missing Any vectors
    @test coerce(Any[4, 7],     Continuous) == [4.0, 7.0]
    @test coerce(Any[4.0, 7.0], Continuous) == [4, 7]

    # Finite conversions:
    @test scitype_union(coerce(['x', 'y'], Finite)) === Multiclass{2}
    @test scitype_union(@test_logs(
        (:info, r"Trying to coerce from `Union{Missing,"),
        coerce(['x', 'y', missing], Finite))) === Union{Missing, Multiclass{2}}

    # More finite conversions (to check resolution of #48):
    y = categorical([1, 2, 3, missing]) # unordered
    yc = coerce(y, Union{Missing,OrderedFactor})
    @test isordered(yc)
    @test yc[1].pool.ordered
    @test scitype(yc) == Vec{Union{Missing, OrderedFactor{3}}}
    @test scitype_union(yc) == Union{Missing, OrderedFactor{3}}
    @test scitype_union(y) == Union{Missing, Multiclass{3}}

    # tests fix for issue https://github.com/JuliaAI/ScientificTypes.jl/issues/161
    X = (x=10:10:44, y=1:4, z=collect("abcd"))
    Xc = coerce(X, :x => Continuous, "y" => Continuous)
    @test scitype_union(Xc.x) === Continuous
    @test scitype_union(Xc.y) === Continuous
end

# issue #62 (ScientficTypes)
@testset "coersion of Type=>Type" begin
    X = (x=[1,2,1,2,5,1,0,7],
         y=[0,1,0,1,0,1,0,1],
         z=['a','b','a','b','a','a',missing,missing])
    Xc = coerce(X, :y=>OrderedFactor)
    Xc = coerce(Xc, Count=>Continuous)
    @test elscitype(Xc.x) == Continuous
    @test elscitype(Xc.y) == OrderedFactor{2}
    Xc = coerce(Xc, OrderedFactor=>Count)
    @test elscitype(Xc.y) == Count
    Xc = coerce(Xc, :z=>Multiclass, verbosity=0)
    Xc = coerce(Xc, Multiclass=>OrderedFactor, verbosity=0)
    @test elscitype(Xc.z) == Union{Missing,OrderedFactor{2}}
    Xc = coerce(X, Count=>Continuous, Unknown=>Multiclass, verbosity=0)
    @test elscitype(Xc.x) == Continuous
    @test elscitype(Xc.y) == Continuous
    @test elscitype(Xc.z) == Union{Missing, Multiclass{2}}
end

# issue #13
@testset "coersion of mixture of Type=>Type and Symbol=>Type" begin
    X = (x=10:10:44, y=1:4, z=collect("abcd"), w=["a", "b", "c", missing])
    @test ScientificTypes.feature_scitype_pairs(:x => Continuous, X) ==
        [:x => Continuous, ]
    @test ScientificTypes.feature_scitype_pairs(Count => Continuous, X) ==
        [:x => Union{Continuous}, :y=> Union{Continuous}]
    X1 = coerce(X, :z => Multiclass, Count=>Continuous)
    @test schema(X1).scitypes ==
        (Continuous, Continuous, Multiclass{4}, Union{Missing,Textual})
    X2 = coerce(X, :z => Multiclass, :w=>Multiclass, Count=>Continuous,
                verbosity=0)
    @test schema(X2).scitypes ==
        (Continuous, Continuous, Multiclass{4}, Union{Missing,Multiclass{3}})
    @test_throws ArgumentError coerce(X, Count => Continuous, :x=>Multiclass)
end

@testset "coerce!" begin
    df = DataFrame((x=ones(Int,5), y=ones(5)))
     @test scitype(df) == Table{Union{AbstractArray{Continuous,1},
                                      AbstractArray{Count,1}}}
    coerce!(df, :x=>Continuous)
    @test scitype(df) == Table{AbstractArray{Continuous,1}}

    df = DataFrame((
          x=ones(Int, 50),
          y=ones(50),
          z=collect("abbabaacbcabbabaacbbcccbccbbabaaacacbcabcbccaabaaa")
          ))
     @test scitype(df) == Table{Union{AbstractArray{Continuous,1},
                                      AbstractArray{Count,1},
                                      AbstractArray{Unknown,1}}}

    coerce!(df, autotype(df, :few_to_finite))
     @test scitype(df) == Table{Union{AbstractArray{Multiclass{3},1},
                                      AbstractArray{OrderedFactor{1},1}}}

    @test_throws ScientificTypes.CoercionError coerce!(randn(5, 5))

    df2 = DataFrame(x=[1,2,3,4],
    y=["a","b","c","a"],
    z = Union{Missing,Int}[10, 20, 30, 40])
    coerce!(df2, Textual=>Finite)
    coerce!(df2, Union{Missing, Count}=>Count, tight=true) # issue #166
    @test scitype(df2) == Table{Union{AbstractArray{Count,1},
                      AbstractArray{Multiclass{3},1} }}

    # issue #9 (coerce Text => Num)
    df3 = DataFrame(x=["1","2","3"], y=[2,3,4])
    coerce!(df3, Textual=>Count)
    @test scitype(df3) == Table{AbstractArray{Count,1}}

    x = [1,2,3,4]
    @test_throws ScientificTypes.CoercionError coerce!(x, Continuous)

end


text:
module aoc_19

export SimpleRule, SeriesRule, CompoundRule, LoopingRule

abstract type AbstractRule end

struct SimpleRule <: AbstractRule
    rule::Regex
end
struct SeriesRule <: AbstractRule
    rules::Vector{Int}
end
struct CompoundRule <: AbstractRule
    rules::Tuple{SeriesRule,SeriesRule}
end
struct LoopingRule <: AbstractRule
    rule::CompoundRule
end

import Base.|
|(a::Regex, b::Regex) = Regex( "(?:" * a.pattern * ")" * "|" * "(?:" * b.pattern * ")" )

const rules_dict = Dict{Int,AbstractRule}()

Regex(r::SimpleRule) = r.rule
Regex(r::SeriesRule) = *( Regex.( [ rules_dict[i] for i in r.rules])... )
Regex(r::CompoundRule) = Regex(r.rules[1]) | Regex(r.rules[2])
function Regex(r::LoopingRule)
    inner = r.rule
    if rules_dict[ inner.rules[2].rules[end] ] == r
        # Rule 8: sub-rule can repeat
        first_rule = Regex(inner.rules[1])
        return Regex( "(?:" * first_rule.pattern * ")+" )
    elseif rules_dict[ inner.rules[2].rules[2] ] == r
        # Rule 11: sub-rules must repeat exactly n times, iterate over n
        first_rule = Regex( rules_dict[ inner.rules[1].rules[1] ])
        last_rule = Regex( rules_dict[ inner.rules[1].rules[2] ])
        # Recursive regex from Doug on Zulip
        out = Regex("((?:$(first_rule.pattern))(?1)?(?:$(last_rule.pattern)))")
        return out
    end
end

import Base.*
*(a::AbstractRule, b::AbstractRule) = SimpleRule( Regex(a) * Regex(b) )

const rule_r = r"([0-9]*):"

const rule_content_r = r"([\d\s$]+)(\|)?([\d\s$]*)"

const simple_rule_r = r"\"([a-z]+)\""

function read_input(filename, clear_rules = true)
    if clear_rules
        empty!(rules_dict)
    end
    lines = readlines(filename)
    # rules_dict = Dict{Int,AbstractRule}()
    signal = Vector{String}()
    for line in lines
        if isempty(line)
            continue
        end
        rule_m = match(rule_r , line)
        if isnothing(rule_m)
            # not a rule
            push!(signal,line)
        else
            rule_no = parse(Int,rule_m.captures[1])
            offset = length(rule_m.match) + 1
            simple_m = match(simple_rule_r, line, offset)
            if isnothing(simple_m)
                content_m = match(rule_content_r, line, offset)
                if content_m.captures[2] == "|"
                    # Compound rule
                    a = SeriesRule( parse.(Int, split( content_m.captures[1] ) ) )
                    b = SeriesRule( parse.(Int, split( content_m.captures[3] ) ) )
                    rules_dict[ rule_no ] = CompoundRule( (a,b) )
                    if rule_no in vcat(a.rules,b.rules)
                        rules_dict[ rule_no ] = LoopingRule( rules_dict[ rule_no ] )
                    end
                else
                    # Series rule
                    rules_dict[ rule_no ] = SeriesRule( parse.(Int, split( content_m.captures[1] ) ) )
                end
            else
                # Simple Rule
                rules_dict[ rule_no ] = SimpleRule( Regex(simple_m.captures[1]) )
            end
        end
    end
    signal
end

function part1()
    data = read_input("input.txt")
    r = Regex("^" * Regex( rules_dict[0] ).pattern * "\$")
    sum( endswith.(data, r) )
    # 102
end

function part2()
    data = read_input("input.txt")
    read_input("newrules.txt",false)
    r = Regex("^" * Regex( rules_dict[0] ).pattern * "\$")
    sum( endswith.(data,r) )
    # 318
end

end


text:
#Axt_mul_Bx

export Axt_mul_Bx

"""
`Axt_mul_Bx(A::AbstractOperator,B::AbstractOperator)`

Create an operator `P` such that:

`P*x == (Ax)'*(Bx)`

# Example

```julia
julia> A,B = randn(4,4),randn(4,4);

julia> P = Axt_mul_Bx(MatrixOp(A),MatrixOp(B))
▒*▒  ℝ^4 -> ℝ^1

julia> x = randn(4);

julia> P*x == [(A*x)'*(B*x)]
true

```
"""
struct Axt_mul_Bx{N,
                  L1 <: AbstractOperator,
                  L2 <: AbstractOperator,
                  C <: AbstractArray,
                  D <: AbstractArray,
                 } <: NonLinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufC::C
  bufD::D
  function Axt_mul_Bx(A::L1, B::L2, bufA::C, bufB::C, bufC::C, bufD::D) where {L1,L2,C,D}
    if ndims(A,1) == 1
      if size(A) != size(B)
        throw(DimensionMismatch("Cannot compose operators"))
      end
    elseif ndims(A,1) == 2 && ndims(B,1) == 2 && size(A,2) == size(B,2)
      if size(A,1)[1] != size(B,1)[1]
        throw(DimensionMismatch("Cannot compose operators"))
      end
    else
      throw(DimensionMismatch("Cannot compose operators"))
    end
    N = ndims(A,1)
    new{N,L1,L2,C,D}(A,B,bufA,bufB,bufC,bufD)
  end
end

struct Axt_mul_BxJac{N,
                     L1 <: AbstractOperator,
                     L2 <: AbstractOperator,
                     C <: AbstractArray,
                     D <: AbstractArray,
                    } <: LinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufC::C
  bufD::D
end

# Constructors
function Axt_mul_Bx(A::AbstractOperator,B::AbstractOperator)
  s,t = size(A,1), codomainType(A)
  bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(B,1), codomainType(B)
  bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(A,2), domainType(A)
  bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  Axt_mul_Bx(A,B,bufA,bufB,bufC,bufD)
end

# Jacobian
function Jacobian(P::Axt_mul_Bx{N,L1,L2,C,D}, x::AbstractArray) where {N,L1,L2,C,D}
  JA, JB = Jacobian(P.A, x), Jacobian(P.B, x)
  Axt_mul_BxJac{N,typeof(JA),typeof(JB),C,D}(JA,JB,P.bufA,P.bufB,P.bufC,P.bufD)
end

# Mappings
# N == 1 input is a vector
function mul!(y, P::Axt_mul_Bx{1,L1,L2,C,D}, b) where {L1,L2,C,D}
  mul!(P.bufA,P.A,b)
  mul!(P.bufB,P.B,b)
  y[1] = dot(P.bufA,P.bufB)
end

function mul!(y, J::AdjointOperator{Axt_mul_BxJac{1,L1,L2,C,D}}, b) where {L1,L2,C,D}
  #y .= conj(J.A.A'*J.A.bufB+J.A.B'*J.A.bufA).*b[1]
  mul!(y, J.A.A', J.A.bufB)
  mul!(J.A.bufD, J.A.B', J.A.bufA)
  y .= conj.( y .+  J.A.bufD ) .* b[1]
  return y
end

# N == 2 input is a matrix
function mul!(y, P::Axt_mul_Bx{2,L1,L2,C,D}, b) where {L1,L2,C,D}
  mul!(P.bufA,P.A,b)
  mul!(P.bufB,P.B,b)
  mul!(y,P.bufA',P.bufB)
  return y
end

function mul!(y, J::AdjointOperator{Axt_mul_BxJac{2,L1,L2,C,D}}, b) where {L1,L2,C,D}
  # y .= J.A.A'*((J.A.bufB)*b') + J.A.B'*((J.A.bufA)*b)
  mul!(J.A.bufC, J.A.bufB, b')
  mul!(y, J.A.A', J.A.bufC)
  mul!(J.A.bufB, J.A.bufA, b)
  mul!(J.A.bufD, J.A.B', J.A.bufB)
  y .+= J.A.bufD
  return y
end

size(P::Union{Axt_mul_Bx{1},Axt_mul_BxJac{1}}) = ((1,),size(P.A,2))
size(P::Union{Axt_mul_Bx{2},Axt_mul_BxJac{2}}) = ((size(P.A,1)[2],size(P.B,1)[2]),size(P.A,2))

fun_name(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = fun_name(L.A)*"*"*fun_name(L.B)

domainType(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = domainType(L.A)
codomainType(L::Union{Axt_mul_Bx,Axt_mul_BxJac}) = codomainType(L.A)

# utils
function permute(P::Axt_mul_Bx{N,L1,L2,C,D},
                 p::AbstractVector{Int}) where {N,L1,L2,C,D <:ArrayPartition}
  Axt_mul_Bx(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p]) )
end

remove_displacement(P::Axt_mul_Bx) =
Axt_mul_Bx(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufC, P.bufD)


text:
# ---
# title: 739. Daily Temperatures
# id: problem739
# author: Indigo
# date: 2021-01-28
# difficulty: Medium
# categories: Hash Table, Stack
# link: <https://leetcode.com/problems/daily-temperatures/description/>
# hidden: true
# ---
#
# Given a list of daily temperatures `T`, return a list such that, for each day
# in the input, tells you how many days you would have to wait until a warmer
# temperature. If there is no future day for which this is possible, put `0`
# instead.
#
# For example, given the list of temperatures `T = [73, 74, 75, 71, 69, 72, 76,
# 73]`, your output should be `[1, 1, 4, 2, 1, 1, 0, 0]`.
#
# **Note:** The length of `temperatures` will be in the range `[1, 30000]`. Each
# temperature will be an integer in the range `[30, 100]`.
#
#
## @lc code=start
using LeetCode

function daily_temperatures(T::Vector{Int})
    res = fill(0, length(T))
    stk = Int[]
    for i in 1:length(T)
        while !isempty(stk) && T[stk[end]] < T[i]
            res[stk[end]] = i - stk[end]
            pop!(stk)
        end
        push!(stk, i)
    end
    res
end
## @lc code=end


text:
module OnlineStats

import StatsBase
import StatsBase: nobs, fit, fit!, skewness, kurtosis, coef, predict
import Distributions; Ds = Distributions
using RecipesBase

export
    OnlineStat,
    # Weight
    Weight, EqualWeight, ExponentialWeight, LearningRate, LearningRate2,
    BoundedEqualWeight,
    # <: OnlineStat
    Mean, Means, Variance, Variances, Extrema, QuantileSGD, QuantileMM, Moments,
    Diff, Diffs, Sum, Sums, CovMatrix, LinReg, LogRegMM, QuantRegMM, NormalMix,
    StatLearn,
    KMeans, BiasVector, BiasMatrix, TwoWayInteractionVector, TwoWayInteractionMatrix,
    FitBeta, FitCategorical, FitCauchy, FitGamma, FitLogNormal, FitNormal,
    FitMultinomial, FitMvNormal,
    # Penalties
    Penalty, NoPenalty, RidgePenalty, LassoPenalty, ElasticNetPenalty,
    # ModelDefinition and Algorithm
    Model, LinearRegression, L1Regression, LogisticRegression,
    PoissonRegression, QuantileRegression, SVMLike, HuberRegression,
    Algorithm, SGD, AdaGrad, AdaGrad2, AdaDelta, ADAM,
    # streamstats
    BernoulliBootstrap, PoissonBootstrap, FrozenBootstrap, cached_state,
    replicates, HyperLogLog,
    # methods
    value, fit, fit!, nobs, skewness, kurtosis, sweep!, coef, predict,
    loss, cost, center, fitdistribution, classify, maprows

#-----------------------------------------------------------------------------# types
abstract Input
abstract ScalarInput    <: Input  # observation = scalar
abstract VectorInput    <: Input  # observation = vector
abstract XYInput        <: Input  # observation = (x, y) pair

abstract OnlineStat{I <: Input}

typealias VecF      Vector{Float64}
typealias MatF      Matrix{Float64}
typealias AVec{T}   AbstractVector{T}
typealias AMat{T}   AbstractMatrix{T}
typealias AVecF     AVec{Float64}
typealias AMatF     AMat{Float64}

#---------------------------------------------------------------------# printing
name(o) = replace(string(typeof(o)), "OnlineStats.", "")
printheader(io::IO, s::AbstractString) = println(io, "■ $s")
function print_item(io::IO, name::AbstractString, value)
    println(io, "  >" * @sprintf("%12s", name * ": "), value)
end
function print_value_and_nobs(io::IO, o::OnlineStat)
    print_item(io, "value", value(o))
    print_item(io, "nobs", nobs(o))
end
function Base.show(io::IO, o::OnlineStat)
    printheader(io, name(o))
    print_value_and_nobs(io, o)
end

#------------------------------------------------------------------------------# fit!
#=
There are so many fit methods because
   - Each method actually needs three implementations (ScalarInput, VectorInput, XYInput)
   - methods:
       - singleton
       - batch
       - singleton + float
       - batch + float
       - batch + vector of floats
       - batch + integer
=#
"""
Update an OnlineStat with more data.  Additional arguments after the input data
provide extra control over how the updates are done.

```
y = randn(100)
o = Mean()

fit!(o, y)      # standard usage

fit!(o, y, 10)  # update in minibatches of size 10

fit!(o, y, .1)  # update using weight .1 for each observation

wts = rand(100)
fit!(o, y, wts) # update observation i using wts[i]
```
"""
############ single observation
function fit!(o::OnlineStat{ScalarInput}, y::Real)
    updatecounter!(o)
    γ = weight(o)
    _fit!(o, y, γ)
    o
end
function fit!{T <: Real}(o::OnlineStat{VectorInput}, y::AVec{T})
    updatecounter!(o)
    γ = weight(o)
    _fit!(o, y, γ)
    o
end
function fit!{T <: Real}(o::OnlineStat{XYInput}, x::AVec{T}, y::Real)
    updatecounter!(o)
    γ = weight(o)
    _fit!(o, x, y, γ)
    o
end

############ single observation, override the weight
function fit!(o::OnlineStat{ScalarInput}, y::Real, γ::Float64)
    updatecounter!(o)
    _fit!(o, y, γ)
    o
end
function fit!{T <: Real}(o::OnlineStat{VectorInput}, y::AVec{T}, γ::Float64)
    updatecounter!(o)
    _fit!(o, y, γ)
    o
end
function fit!{T <: Real}(o::OnlineStat{XYInput}, x::AVec{T}, y::Real, γ::Float64)
    updatecounter!(o)
    _fit!(o, x, y, γ)
    o
end

############ multiple observations
function fit!(o::OnlineStat{ScalarInput}, y::AVec)
    for yi in y
        fit!(o, yi)
    end
    o
end
function fit!(o::OnlineStat{VectorInput}, y::AMat)
    for i in 1:size(y, 1)
        fit!(o, row(y, i))
    end
    o
end
function fit!(o::OnlineStat{XYInput}, x::AMat, y::AVec)
    @assert size(x, 1) == length(y)
    for i in eachindex(y)
        fit!(o, row(x, i), row(y, i))
    end
    o
end

############ multiple observations, override weight
function fit!(o::OnlineStat{ScalarInput}, y::AVec, w::AVec)
    @assert length(y) == length(w)
    for i in eachindex(y)
        fit!(o, row(y, i), w[i])
    end
    o
end
function fit!(o::OnlineStat{VectorInput}, y::AMat, w::AVec)
    n2 = nrows(y)
    @assert n2 == length(w)
    for i in 1:n2
        fit!(o, row(y, i), w[i])
    end
    o
end
function fit!(o::OnlineStat{XYInput}, x::AMat, y::AVec, w::AVec)
    @assert size(x, 1) == length(y) == length(w)
    for i in eachindex(y)
        fit!(o, row(x, i), row(y, i), w[i])
    end
    o
end

############ multiple observations, override weight, each update gets the same weight
function fit!(o::OnlineStat{ScalarInput}, y::AVec, w::Real)
    for i in eachindex(y)
        fit!(o, y[i], w)
    end
    o
end
function fit!(o::OnlineStat{VectorInput}, y::AMat, w::Real)
    n2 = nrows(y)
    for i in 1:n2
        fit!(o, row(y, i), w)
    end
    o
end
function fit!(o::OnlineStat{XYInput}, x::AMat, y::AVec, w::Real)
    @assert size(x, 1) == length(y)
    for i in eachindex(y)
        fit!(o, row(x, i), row(y, i), w)
    end
    o
end

############ fit with observations in the columns (experimental)
function fit_col!(o::OnlineStat{VectorInput}, y::AMat)
    for i in 1:size(y, 2)
        fit!(o, col(y, i))
    end
    o
end
function fit_col!(o::OnlineStat{XYInput}, x::AMat, y::AVec)
    @assert size(x, 2) == length(y)
    for i in eachindex(y)
        fit!(o, col(x, i), row(y, i))
    end
    o
end

############ multiple observations, update in batches
function fit!(o::OnlineStat{ScalarInput}, y::AVec, b::Integer)
    b = Int(b)
    n = length(y)
    0 < b <= n || warn("batch size larger than data size")
    if b == 1
        fit!(o, y)
    else
        i = 1
        while i <= n
            rng = i:min(i + b - 1, n)
            bsize = length(rng)
            updatecounter!(o, bsize)
            γ = weight(o, bsize)
            _fitbatch!(o, rows(y, rng), γ)
            i += b
        end
    end
    o
end
function fit!(o::OnlineStat{VectorInput}, y::AMat, b::Integer)
    b = Int(b)
    n = size(y, 1)
    0 < b <= n || warn("batch size larger than data size")
    if b == 1
        fit!(o, y)
    else
        i = 1
        while i <= n
            rng = i:min(i + b - 1, n)
            bsize = length(rng)
            updatecounter!(o, bsize)
            γ = weight(o, bsize)
            _fitbatch!(o, rows(y, rng), γ)
            i += b
        end
    end
    o
end
function fit!(o::OnlineStat{XYInput}, x::AMat, y::AVec, b::Integer)
    b = Int(b)
    n = length(y)
    0 < b <= n || warn("batch size larger than data size")
    if b == 1
        fit!(o, x, y)
    else
        i = 1
        while i <= n
            rng = i:min(i + b - 1, n)
            bsize = length(rng)
            updatecounter!(o, bsize)
            γ = weight(o, bsize)
            _fitbatch!(o, rows(x, rng), rows(y, rng), γ)
            i += b
        end
    end
    o
end

# error if no fitbatch! method
_fitbatch!(o, args...) = (warn("no fitbatch! method...calling fit!"); _fit!(o, args...))

#---------------------------------------------------------------------------# helpers
"""
The associated value of an OnlineStat.

```
o = Mean()
value(o)
```
"""
value(o::OnlineStat) = o.value
StatsBase.nobs(o::OnlineStat) = nobs(o.weight)
unbias(o::OnlineStat) = nobs(o) / (nobs(o) - 1)

# for updating
smooth(m::Float64, v::Real, γ::Float64) = (1.0 - γ) * m + γ * v
function smooth!(m::AVec, v::AVec, γ::Float64)
    for i in eachindex(v)
        @inbounds m[i] = smooth(m[i], v[i], γ)
    end
end
subgrad(m::Float64, γ::Float64, g::Real) = m - γ * g
function smooth!(avg::AMat, v::AMat, λ::Float64)
    n, p = size(avg)
    @assert size(avg) == size(v)
    for j in 1:p, i in 1:n
        @inbounds avg[i,j] = smooth(avg[i, j], v[i, j], λ)
    end
end
# Rank 1 update of symmetric matrix: (1 - γ) * A + γ * x * x'
function rank1_smooth!(A::AMat, x::AVec, γ::Float64)
    @assert size(A, 1) == size(A, 2)
    for j in 1:size(A, 2), i in 1:j
        @inbounds A[i, j] = (1.0 - γ) * A[i, j] + γ * x[i] * x[j]
    end
end
# # Why doesn't this work?  Tested with CovMatrix
# function rank1_smooth!(A::AMat, x::AVec, γ::Float64)
#     scale!(A, 1.0 - γ)
#     BLAS.syr!('U', γ, x, A)
# end

row(x::AMat, i::Integer) = view(x, i, :)
row(x::AVec, i::Integer) = x[i]
rows(x::AMat, rs::AVec{Int}) = view(x, rs, :)
rows(x::AVec, rs::AVec{Int}) = view(x, rs)
col(x::AMat, i::Integer) = view(x, :, i)
cols(x::AMat, rs::AVec{Int}) = view(x, :, rs)

nrows(x::AMat) = size(x, 1)
ncols(x::AMat) = size(x, 2)

Base.copy(o::OnlineStat) = deepcopy(o)

#-----------------------------------------------------------------------------# merge
function Base.merge(o::OnlineStat, o2::OnlineStat, method::Symbol = :append)
    merge!(copy(o), o2, method)
end
function Base.merge(o::OnlineStat, o2::OnlineStat, wt::Float64)
    merge!(copy(o), o2, wt)
end

function Base.merge!(o::OnlineStat, o2::OnlineStat, method::Symbol = :append)
    @assert typeof(o) == typeof(o2)
    updatecounter!(o, nobs(o2))
    if method == :append
        _merge!(o, o2, weight(o, nobs(o2)))
    elseif method == :mean
        _merge!(o, o2, 0.5 * (weight(o) + weight(o2)))
    elseif method == :singleton
        _merge!(o, o2, weight(o))
    end
    o
end

function Base.merge!(o::OnlineStat, o2::OnlineStat, wt::Float64)
    @assert typeof(o) == typeof(o2)
    updatecounter!(o, nobs(o2))
    _merge!(o, o2, wt)
    o
end

# epsilon used in special cases to avoid dividing by 0, etc.
const _ϵ = 1e-8

#---------------------------------------------------------------------------# maprows
"""
Perform operations on data in blocks.

`maprows(f::Function, b::Integer, data...)`

This function iteratively feeds `data` in blocks of `b` observations to the
function `f`.  The most common usage is with `do` blocks:

```julia
# Example 1
y = randn(50)
o = Variance()
maprows(10, y) do yi
    fit!(o, yi)
    println("Updated with another batch!")
end
```
"""
function maprows(f::Function, b::Integer, data...)
    n = size(data[1], 1)
    i = 1
    while i <= n
        rng = i:min(i + b - 1, n)
        batch_data = map(x -> rows(x, rng), data)
        f(batch_data...)
        i += b
    end
end

#----------------------------------------------------------------------# source files
include("weight.jl")
include("summary.jl")
include("distributions.jl")
include("normalmix.jl")
include("modeling/temp.jl")
include("modeling/statlearn.jl")
include("modeling/linreg.jl")
include("modeling/logregmm.jl")
include("modeling/quantregmm.jl")
include("modeling/bias.jl")
include("streamstats/bootstrap.jl")
include("streamstats/hyperloglog.jl")
include("multivariate/kmeans.jl")
include("plots.jl")

end # module


text:
if !isdefined(Main, :Groebner)
    import Groebner
end

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

function my_system_coeffs_1(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    gb = Groebner.groebner(system, reduced=false)

    aver = [(0, 0) for i in 1:length(gb)]
    for (i, poly) in enumerate(gb)
        num = sum([Base.GMP.MPZ.sizeinbase(numerator(c), 2) for c in AbstractAlgebra.coefficients(poly)])
        den = sum([Base.GMP.MPZ.sizeinbase(denominator(c), 2) for c in AbstractAlgebra.coefficients(poly)])

        aver[i] = (num, den)
    end

    Groebner.printall()

    # println(aver)

    println("Numerators: ", sum(x[1] for x in aver))
    println("Denominators: ", sum(x[2] for x in aver))

end

function my_system_coeffs_2(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    gb = Groebner.groebner(system, reduced=false)

    aver = [(0, 0) for i in 1:length(gb)]
    for (i, poly) in enumerate(gb)
        @warn "new poly" i
        println([(Base.GMP.MPZ.sizeinbase(numerator(c), 2), Base.GMP.MPZ.sizeinbase(denominator(c), 2)) for c in AbstractAlgebra.coefficients(poly)])
    end

    Groebner.printall()

    # println(aver)

    println("Numerators: ", sum(x[1] for x in aver))
    println("Denominators: ", sum(x[2] for x in aver))

end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
        #("cyclic 7", Groebner.cyclicn(7, ground=ground)),
        #("noon 8"    ,Groebner.noonn(8, ground=ground)),
        # ("eco 11"    ,Groebner.eco11(ground=ground)),
        ("katsura 8"    ,Groebner.katsuran(8, ground=ground)),
        #("root 11"    ,Groebner.rootn(11, ground=ground)),
        ]

    println()
    for (name, system) in systems
        println("-"^20)
        println("$name")
        # benchmark_system_my(system)
        my_system_coeffs_2(system)
        println("-"^20)
    end
end

ground = AbstractAlgebra.QQ
run_f4_ff_degrevlex_benchmarks(ground)

#=
cyclic 11
  50.711 ms (191971 allocations: 23.54 MiB)
noon 7
  910.447 ms (2087052 allocations: 298.79 MiB)
eco 10
  1.484 s (1994909 allocations: 417.28 MiB)
katsura 8
  4.509 s (2982135 allocations: 499.90 MiB)
=#


text:
using Lazy
import Lazy: cycle, range, drop, take
using Base.Test

@testset "Lazy" begin

@testset "Lists" begin
    @test list(1, 2, 3)[2] == 2
    @test prepend(1, list(2,3,4)) == 1:list(2, 3, 4)
    @test seq([1, 2, 3]) == list(1, 2, 3)
    @test seq(1:3) == list(1, 2, 3)
    @test constantly(1)[50] == 1
    testfn() = 1
    @test repeatedly(testfn)[50] == 1
    @test cycle([1, 2, 3])[50] == 2
    @test iterate(x->x^2, 2)[3] == 16
    @test range(1, 5)[3] == 3
    @test range(1, 5)[10] == nothing
    @test range(1, 5)[-1] == 1
    @test list(1, 2, 3) * list(4, 5, 6) == list(1, 2, 3, 4, 5, 6)
    @test first(list(1, 2, 3)) == 1
    @test tail(list(1, 2, 3)) == list(2, 3)
    @test flatten(list(1,2,list(3,4))) == list(1, 2, 3, 4)
    @test list(1,2,list(3,4))[3] == list(3, 4)
end

@testset "Fibs" begin
    fibs = @lazy 0:1:(fibs + drop(1, fibs));
    @test fibs[20] == 4181
    @test take(5, fibs) == list(0, 1, 1, 2, 3)
end

@testset "Primes" begin
    isprime(n) =
        @>> primes begin
            take_while(x -> x<=sqrt(n))
            map(x -> n % x == 0)
            any; !
        end
    primes = filter(isprime, range(2));
end

@testset "Even squares" begin
    esquares = @>> range() map(x->x^2) filter(iseven);
    @test take(5, esquares) == list(4, 16, 36, 64, 100)
end

@testset "Threading macros" begin
    temp = @> [2 3] sum
    @test temp == 5
    # Reverse from after index 2
    temp = @>> 2 reverse([1, 2, 3, 4, 5])
    @test temp == [1, 5, 4, 3, 2]
    temp = @as x 2 begin
        x^2
        x + 2
    end
    @test temp == 6
end

@testset "Listables" begin
    @test_throws MethodError sin()
end

@testset "any/all" begin
    let xs = list(true, false, false)
        @test any(identity, xs) == true
        @test any(xs) == true
        @test all(identity, xs) == false
        @test all(xs) == false
    end
    let yy = list(1, 0, 1)
        @test any(Bool, yy) == true
        @test all(Bool, yy) == false
    end
    # Base method--ensures no ambiguity with methods here
    @test all([true true; true true], 1) == [true true]
end

end


text:
"""
	grdfft(cmd0::String="", arg1=nothing, [arg2=nothing,] kwargs...)

Take the 2-D forward Fast Fourier Transform and perform one or more mathematical operations
in the frequency domain before transforming back to the space domain.

Full option list at [`grdfft`]($(GMTdoc)grdfft.html)

Parameters
----------

- **A** | **azim** :: [Type => Number]    ``Arg = azim``

    Take the directional derivative in the azimuth direction measured in degrees CW from north.
    ($(GMTdoc)grdfft.html#a)
- **C** | **upward** :: [Type => Number]    ``Arg = zlevel``

    Upward (for zlevel > 0) or downward (for zlevel < 0) continue the field zlevel meters.
    ($(GMTdoc)grdfft.html#c)
- **D** | **dfdz** :: [Type => Str or Number]		``Arg = [scale|g]``

    Differentiate the field, i.e., take d(field)/dz. This is equivalent to multiplying by kr in
    the frequency domain (kr is radial wave number).
    ($(GMTdoc)grdfft.html#d)
- **E** | **radial_power** :: [Type => Str]         ``Arg = [r|x|y][+w[k]][+n]``

    Estimate power spectrum in the radial direction [r]. Place x or y immediately after E to
    compute the spectrum in the x or y direction instead.
    ($(GMTdoc)grdfft.html#e)
- **F** | **filter** :: [Type => Str or List--        ``Arg = [r|x|y]params``

    Filter the data. Place x or y immediately after -F to filter x or y direction only; default is
    isotropic [r]. Choose between a cosine-tapered band-pass, a Gaussian band-pass filter, or a
    Butterworth band-pass filter.
    ($(GMTdoc)grdfft.html#f)
- **G** | **save** | **outgrid** | **outfile** :: [Type => Str]

    Output grid file name (or table if **radial_power** is used). Note that this is optional and to
    be used only when saving the result directly on disk. Otherwise, just use the G = grdfft(....) form.
    ($(GMTdoc)grdfft.html#g)
- **I** | **integrate** :: [Type => Str or Number]		``Arg = [scale|g]``

    Integrate the field, i.e., compute integral_over_z (field * dz). This is equivalent to divide
    by kr in the frequency domain (kr is radial wave number).
    ($(GMTdoc)grdfft.html#i)
- **N** | **inquire** :: [Type => Str]         ``Arg = [a|f|m|r|s|nx/ny][+a|[+d|h|l][+e|n|m][+twidth][+v][+w[suffix]][+z[p]]``

    Choose or inquire about suitable grid dimensions for FFT and set optional parameters. Control the FFT dimension:
    ($(GMTdoc)grdfft.html#n)
- **S** | **scale** :: [Type => Number]			``Arg = scale``

    Multiply each element by scale in the space domain (after the frequency domain operations).
    ($(GMTdoc)grdfft.html#s)
- $(GMT.opt_V)
- $(GMT.opt_f)
"""
function grdfft(cmd0::String="", arg1=nothing, arg2=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdfft", cmd0, arg1, arg2)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:G :V_params :f])
	cmd  = parse_these_opts(cmd, d, [[:A :azim], [:C :upward], [:D :dfdz], [:E :radial_power], [:F :filter],
	                                 [:I :integrate], [:N :inquire], [:S :scale]])

	cmd, got_fname, arg1, arg2 = find_data(d, cmd0, cmd, arg1, arg2)
	(isa(arg1, Array{<:Number})) && (arg1 = mat2grid(arg1))
	if (!occursin(" -E", cmd))          # Simpler case
		return common_grd(d, "grdfft " * cmd, arg1)		# Finish build cmd and run it
	else
		# Here several cases can happen: 1) arg1 only; 2) arg1 && arg2; 3) grid(s) provided via fname
		(isa(arg2, Array{<:Number})) && (arg2 = mat2grid(arg2))
		return common_grd(d, "grdfft " * cmd, arg1, arg2)
	end
end

# ---------------------------------------------------------------------------------------------------
grdfft(arg1, arg2=nothing, cmd0::String=""; kw...) = grdfft(cmd0, arg1, arg2; kw...)


text:
export enable_inline_stability_checks, inline_stability_checks_enabled
export @stable_function, stability_warn

run_inline_stability_checks = false

"""
    enable_inline_stability_checks(::Bool)

Sets whether to run inline stability checks from [`@stable_function`](@ref).

If it is set to `false` (the default value), @stable_function does not perform
any type stability checks.

The value is checked when @stable_function is evaluated, so this should useually
be set at the begining of a package definition.

See [`inline_stability_checks_enabled`](@ref).
"""
function enable_inline_stability_checks(enabled::Bool)
    global run_inline_stability_checks
    run_inline_stability_checks = enabled
end

"""
    inline_stability_checks_enabled()::Bool

Returns whether inline stability checks are enabled.

See [`enable_inline_stability_checks`](@ref).
"""
function inline_stability_checks_enabled()::Bool
    run_inline_stability_checks
end

"""
    @stable_function arg_lists function_name
    @stable_function arg_lists function_definition(s)
    @stable_function arg_lists acceptable_instability function_name
    @stable_function arg_lists acceptable_instability function_definitions(s)

Checks the type stability of the function under the given argument lists.

If the second value is a function definition, the function is defined before
checking type stability.
"""
macro stable_function(arg_lists, unstable, func)
    if run_inline_stability_checks
        if unstable == nothing || unstable == :nothing
            unstable = Dict{Symbol, Type}()
        end
        if VERSION < v"0.7.0-"
            (func_names, body) = parsebody(func, @__MODULE__)
        else
            (func_names, body) = parsebody(func, __module__)
        end
        esc(quote
            $body
            $((:(TypeStability.stability_warn($name, TypeStability.check_function($name, $arg_lists, $unstable)))
               for name in func_names)...)
        end)
    else
        esc(func)
    end
end

macro stable_function(arg_lists, func)
    esc(:(TypeStability.@stable_function $arg_lists nothing $func))
end

"""
    parsebody(body, mod::Any)

Internal method to parse the last argument of @stable_function
`mod` should be a Module in Julia 0.7+, but don't matter on Julia 0.6 due to changes in `macroexpand`
"""
function parsebody(body::Expr, mod; require_function=true)
    if body.head == :function || (body.head == :(=) && isa(body.args[1], Expr))
        if body.args[1] isa Symbol
            func_names = [body.args[1]]
        else
            subExpr = body.args[1]
            if subExpr.head == :where
                subExpr = subExpr.args[1]
            end
            if subExpr.head == :(::)
                subExpr = subExpr.args[1]
            end
            if subExpr.head == :call
                func_names = [subExpr.args[1]]
            else
                error("Cannot find function name in $(dump(body))")
            end
        end
    elseif body.head == :macrocall
        if VERSION < v"0.7.0-"
            expanded_body = macroexpand(body)
        else
            expanded_body = macroexpand(mod, body)
        end
        if isa(expanded_body, Expr)
            (func_names, _) = parsebody(expanded_body, mod; require_function=false)
        elseif isa(expanded_body, Symbol)
            func_names = [expanded_body]
        elseif require_function
            error("Cannot find a function name in macro expansion of $body")
        else
            func_names = Symbol[]
        end
    elseif body.head == :block
        func_names = Symbol[]
        for expr in body.args
            if isa(expr, Expr)
                (expr_func_names, _) = parsebody(expr, mod; require_function=false)
                append!(func_names, expr_func_names)
            end
        end
        func_names = unique(func_names)
        if require_function && length(func_names) == 0
            error("Cannot find any function names in $body")
        end
    elseif require_function
        error("Don't know how to find function names in $body")
    else
        func_names = Symbol[]
    end
    (func_names, body)
end

function parsebody(func::Symbol, mod)
    ([func], quote end)
end


text:
module Covid19

using HTTP, Gumbo, CSV, Cascadia, DataFrames, JSON
using Cascadia:matchFirst

include("DatesInGerman.jl")

const URL = "https://www.hamburg.de/corona-zahlen/"
const CSV_INFECTED = joinpath(@__DIR__, "infected.csv")
const CSV_BOROUGHS = joinpath(@__DIR__, "boroughs.csv")
const CSV_AGEGROUPS = joinpath(@__DIR__, "agegroups.csv")
const JSON_INFECTED = joinpath(@__DIR__, "infected.json")

function fetch()
    response = HTTP.get(URL)
    parsehtml(String(response))
end

function parseall(html)
    infected = parseinfected(html.root)
    vaccinations = parsevaccinated(html.root)
    recordedat = parsedateinfected(html.root)
    deaths = parsedeaths(html.root)
    hospitalizations = parsehospitalizations(html.root)
    trend = parsetrend(html.root)
    boroughs = parseboroughs(html.root)
    agegroups = parseagegroups(html.root)

    Dict(:infected => infected,
       :deaths => Dict(:new => deaths[1], :total => deaths[2]),
       :hospitalizations => Dict(:total => hospitalizations[1], :intensivecare => hospitalizations[2]),
       :vaccinations => Dict(:first_vaccination => vaccinations[1], :second_vaccination => vaccinations[2]),
       :trend => trend,
       :boroughs => boroughs,
       :agegroups => agegroups)
end

function build(current)
    Dict(:infected => buildinfected(current),
         :boroughs => buildboroughs(current),
         :agegroups => buildagegroups(current))
end

function parseinfected(root)
    parsed = map(parsenumbers, eachmatch(sel".nav-main__wrapper .dashboar_number", root)[1:3])
    recordedat = parsedateinfected(root)
    Dict(:new => parsed[1], :total => parsed[2], :recovered => parsed[3], :recordedat => recordedat)
end

function parsedateinfected(root)
    daterecorded = matchFirst(sel".chart_publication", root)[1].text
    DatesInGerman.parsefrom(daterecorded, inwords=false)
end

function parsevaccinated(root)
    map(parsenumbers, eachmatch(sel".nav-main__wrapper .dashboar_number", root)[5:6])
end

function parsedateboroughs(root)
    daterecorded = eachmatch(sel".table-article + p", root)[end][1].text
    DatesInGerman.parsefrom(daterecorded, inwords=false)
end

function parsedeaths(root)
    map(parsenumbers, eachmatch(sel".nav-main__wrapper .dashboar_number", root)[9:10])
end

function parsehospitalizations(root)
    map(parsenumbers, eachmatch(sel".nav-main__wrapper .dashboar_number", root)[7:8])
end

function parsetrend(root)
    map(el -> parse(Int, el[1].text), eachmatch(sel".cv_chart_container .value_show", root))
end

function parseboroughs(root)
    rows = eachmatch(sel"table tr", root)[18:end]
    mapped = Dict{String,Any}()
    foreach(rows) do row
        name = matchFirst(sel"td:first-child", row)[1].text
        num = parse(Int, matchFirst(sel"td:last-child", row)[1].text)
        mapped[name] = num
    end
    mapped["recordedat"] = parsedateboroughs(root)
    mapped
end

function parseagegroups(root)
    rows = eachmatch(sel"table tr", root)[6:16]
    daterecorded = eachmatch(sel"table+p", root)[1][1].text |> DatesInGerman.parsefrom

    map(rows) do row
        age = matchFirst(sel"[data-label=\"Alter\"]", row)[1].text |> parseage |> string
        male = parse(Int, row[2][1].text)
        female = parse(Int, row[3][1].text)
        Dict(:male => male, :female => female, :age => age, :recordedat => daterecorded)
    end
end

function parseage(age::String)::Union{Int,UnitRange{Int}}
    agerange = match(r"^(\d+).+?(\d+)", age)
    if ! isnothing(agerange)
        return parse(Int, agerange.captures[1]):parse(Int, agerange.captures[2])
    end

    young = match(r"^bis.(\d+).+$", age)
    if ! isnothing(young)
        return 0:parse(Int, young.captures[1])
    end

    old = match(r"über.(\d+).*$", age)
    if ! isnothing(old)
        return parse(Int, old.captures[1])
    end
    throw(ArgumentError("Age $age could not be parsed."))
end

function parsenumbers(el)
    text = el[1].text
    parse(Int, match(r"\d+", text).match)
end

function buildinfected(current)::DataFrame
    infected = current[:infected]
    infected[:deaths] = current[:deaths][:total]
    infected[:hospitalizations] = current[:hospitalizations][:total]
    infected[:intensivecare] = current[:hospitalizations][:intensivecare]
    infected[:first_vaccination] = current[:vaccinations][:first_vaccination]
    infected[:second_vaccination] = current[:vaccinations][:second_vaccination]
    df = DataFrame(infected)
    persisted = CSV.read(CSV_INFECTED, DataFrame)
    unique(vcat(df, persisted), :recordedat)
end

function buildboroughs(current)::DataFrame
    df = DataFrame(current[:boroughs])
    persisted = CSV.read(CSV_BOROUGHS, DataFrame)
    unique(vcat(df, persisted), :recordedat)
end

function buildagegroups(current)::DataFrame
    df = DataFrame(current[:agegroups])
    persisted = CSV.read(CSV_AGEGROUPS, DataFrame)
    unique(vcat(df, persisted), [:recordedat, :age])
end

function record()::Nothing
    datasets = fetch() |> parseall |> build
    recordinfected(datasets[:infected])
    recordboroughs(datasets[:boroughs])
    recordagegroups(datasets[:agegroups])
    nothing
end

function recordinfected(df)
    open(JSON_INFECTED, "w") do f
        write(f, JSON.json(df))
    end
    df |> CSV.write(CSV_INFECTED)
end

function recordboroughs(df)
    CSV.write(CSV_BOROUGHS, df)
end

function recordagegroups(df)
    CSV.write(CSV_AGEGROUPS, df)
end

end # module


text:
module TensorGraphicalModels

using SparseArrays
using LinearAlgebra
using Kronecker
using SpecialMatrices
using BandedMatrices
using Distributions
using Random
using DataFrames
using MatrixEquations
using PDMats
using TensorToolbox
using LightGraphs
using GLMNet
using Arpack
using Statistics
using Debugger
using Printf
using Plots

include("utils.jl")
include("glasso.jl")
include("kglasso.jl")
include("sg_palm.jl")
include("kron_pca.jl")
include("utils_teralasso.jl")
include("teralasso.jl")
include("utils_sim.jl")
include("enkf.jl")

export syglasso_palm
export kglasso
export teralasso
export kron_pca
export robust_kron_pca
export glasso
export enkf

end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function acoustic_one_step(w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny)
    acoustic_one_step_ = load_op_and_grad("../build/libADSeismic.so","acoustic_one_step_cpu", multiple=true)
    w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny = convert_to_tensor(Any[w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny], [Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Int64,Int64])
    acoustic_one_step_(w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny)
end

# TODO: specify your input parameters
nx = 10
ny = 10
mask = ones(nx+2, ny+2)
mask = mask[:]

using Random; Random.seed!(1234)
hx = constant(0.1)
hy = constant(0.1)
dt = constant(0.1)
c = rand((nx+2)*(ny+2)) .* mask
phi = rand((nx+2)*(ny+2)) .* mask
psi = rand((nx+2)*(ny+2)) .* mask
wold = rand((nx+2)*(ny+2)) .* mask
w = rand((nx+2)*(ny+2)) .* mask
sigma = rand((nx+2)*(ny+2)) .* mask
tau = rand((nx+2)*(ny+2)) .* mask
u = acoustic_one_step(w,wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny)
sess = Session(); init(sess)
@show run(sess, u)

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    return sum(sum(acoustic_one_step(w, wold,phi,psi,sigma,tau,c,dt,hx,hy,nx,ny).^2))
end
# w, psi

# psi, w

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand((nx+2)*(ny+2)) .* mask)
v_ = rand((nx+2)*(ny+2)) .* mask
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 0.01 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")
savefig("gradtest")


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Types and structures for the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export TextFormat, Highlighter, CustomTextCell

"""
    TextFormat

# Fields

- `up_right_corner::Char`: Character in the up right corner.
- `up_left_corner::Char`: Character in the up left corner.
- `bottom_left_corner::Char`: Character in the bottom left corner.
- `bottom_right_corner::Char`: Character in the bottom right corner.
- `up_intersection::Char`: Character in the intersection of lines in the up
    part.
- `left_intersection::Char`: Character in the intersection of lines in the left
    part.
- `right_intersection::Char`: Character in the intersection of lines in the
    right part.
- `middle_intersection::Char`: Character in the intersection of lines in the
    middle of the table.
- `bottom_intersection::Char`: Character in the intersection of the lines in the
    bottom part.
- `column::Char`: Character in a vertical line inside the table.
- `left_border::Char`: Character used as the left border.
- `right_border::Char`: Character used as the right border.
- `row::Char`: Character in a horizontal line inside the table.
- `hlines::Vector{Symbol}`: Horizontal lines that must be drawn by default.
- `vlines::Union{Symbol, Vector{Symbol}}`: Vertical lines that must be drawn by
    default.

# Pre-defined formats

The following pre-defined formats are available: `unicode` (**default**),
`mysql`, `compact`, `markdown`, `simple`, `ascii_rounded`, and `ascii_dots`.
"""
@kwdef struct TextFormat
    up_right_corner::Char                  = '┐'
    up_left_corner::Char                   = '┌'
    bottom_left_corner::Char               = '└'
    bottom_right_corner::Char              = '┘'
    up_intersection::Char                  = '┬'
    left_intersection::Char                = '├'
    right_intersection::Char               = '┤'
    middle_intersection::Char              = '┼'
    bottom_intersection::Char              = '┴'
    column::Char                           = '│'
    row::Char                              = '─'
    hlines::Vector{Symbol}                 = [:begin, :header, :end]
    vlines::Union{Symbol, Vector{Symbol}}  = :all
end

"""
    Highlighter

Defines the default highlighter of a table when using the text backend.

# Fields

- `f::Function`: Function with the signature `f(data, i, j)` in which should
    return `true` if the element `(i,j)` in `data` must be highlighter, or
    `false` otherwise.
- `fd::Function`: Function with the signature `f(h, data, i, j)` in which `h` is
    the highlighter. This function must return the `Crayon` to be applied to the
    cell that must be highlighted.
- `crayon::Crayon`: The `Crayon` to be applied to the highlighted cell if the
    default `fd` is used.

# Remarks

This structure can be constructed using three helpers:

    Highlighter(f::Function; kwargs...)

where it will construct a `Crayon` using the keywords in `kwargs` and apply it
to the highlighted cell,

    Highlighter(f::Function, crayon::Crayon)

where it will apply the `crayon` to the highlighted cell, and

    Highlighter(f::Function, fd::Function)

where it will apply the `Crayon` returned by the function `fd` to the
highlighted cell.
"""
@kwdef struct Highlighter
    f::Function
    fd::Function = (h, data, i, j) -> h.crayon

    # Private
    crayon::Crayon = Crayon()
end

Highlighter(f; kwargs...) = Highlighter(f = f, crayon = Crayon(;kwargs...))
Highlighter(f, crayon::Crayon) = Highlighter(f = f, crayon = crayon)
Highlighter(f::Function, fd::Function) = Highlighter(f = f, fd = fd)

"""
    Display

Store the information of the display and the current cursor position.

!!! note
    This is not the real cursor position with respect to the display, but with
    respect to the point in which the table is printed.

# Fields

- `size::Tuple{Int, Int}`: Display size.
- `row::Int`: Current row.
- `col::Int`: Current column.
- `has_color::Bool`: Indicates if the display has color support.
- `cont_char::Char`: The character that indicates the line is cropped.
- `cont_space_char::Char`: Space character to be printed before `cont_char`.
"""
@kwdef mutable struct Display
    size::Tuple{Int,Int}  = (-1, -1)
    row::Int              = 1
    col::Int              = 0
    has_color::Bool       = false
    cont_char::Char       = '⋯'
    cont_space_char::Char = ' '

    # Buffer that will store the current line.
    buf_line::IOBuffer = IOBuffer()
end

"""
    CustomTextCell

Abstract type of custom cells in the text backend.

Each type must implement the following API:

- `get_printable_cell_text`: A function that must return a vector of strings
    with the printable text, *i.e.* without any non-printable character.
- `get_rendered_line`: A function that must return the rendered line that will
    be printed to the display.
- `apply_line_padding!`: Apply a certain number of spaces to the left and right
    of a specific line.
- `crop_line!`: A function that must crop a certain number of printable
    characters from the end of the line.

- `append_suffix_to_line!`: Append a string suffix to a line of the custom cell.
- `apply_line_padding!`: Apply left and right padding to a line of the custom
    cell.
- `crop_line!`: Crop a certain number of characters from a line of the custom
    cell.
- `get_printable_cell_line`: Get a printable line of the custom cell.
- `get_rendered_line`: Get a rendered line of the custom cell.
- `parse_cell_text`: Parse the cell text and return a `Vector{String}` with the
    printable lines.
- `reset!`: Reset all the temporary fields. This function is not required.
"""
abstract type CustomTextCell end

################################################################################
#                                  Constants
################################################################################

# Crayon used to reset all the styling.
const _default_crayon = Crayon()
const _reset_crayon   = Crayon(reset = true)


text:
export computeGradMisfit

function computeGradMisfit(sig,Dc::Array,pMis::MisfitParam)
	#
	#	gc = computeGradMisfit(mc,model,Dc,Dobs,Wd,misfit,pFor)
	#
	#	computes gradient of misfit for single forward problem
	#
	#	Note: all variables have to be in memory of the worker executing this method
	#
	try
		sigma,dsigma = pMis.modelfun(sig)
		sigmaloc = interpGlobalToLocal(sigma,pMis.gloc.PForInv,pMis.gloc.sigmaBackground)
		F,dF,d2F = pMis.misfit(Dc,pMis.dobs,pMis.Wd)
		return dsigma'*interpLocalToGlobal(getSensTMatVec(dF,sigmaloc,pMis.pFor),pMis.gloc.PForInv)
	catch err
		if isa(err,InterruptException)
			return -1
		else
			throw(err)
		end
	end
end

function computeGradMisfit(sigmaRef::RemoteChannel,
                           DcRef::Future,
                           pMisRef::RemoteChannel,
                           dFiRef::RemoteChannel)
	#
	#	gc = computeGradMisfit(mc,model,Dc,Dobs,Wd,misfit,pFor)
	#
	#	compute gradient of misfit for single forward problem
	#
	#	Note: model and forward problem are represented as RemoteReferences.
	#
	#	!! make sure that everything is stored on this worker to avoid communication !!
	#
	rrlocs = [sigmaRef.where pMisRef.where DcRef.where]
	if !all(rrlocs .== myid())
		warn("computeGradMisfit: Problem on worker $(myid()) not all remote refs are stored here, but rrlocs=$rrlocs")
	end

	t = time_ns();
	pMis  = take!(pMisRef) # this is a no-op if pFor is stored on this worker
	sigma = fetch(sigmaRef)
	Dc    = fetch(DcRef)
    finalize(DcRef)  # to prevent memory leak
	commTime = (time_ns() - t)/1e+9;

	compTime = @elapsed dFt  = computeGradMisfit(sigma,Dc,pMis)

	t = time_ns();
	dFi  = take!(dFiRef)
	put!(dFiRef,dFi+dFt)
	put!(pMisRef,pMis)    # does not require communication if PF lives on this worker
	commTime += (time_ns() - t)/1e+9;

#	if commTime/compTime > 1.0
#		warn("computeGradMisfit: Communication time larger than computation time! commTime/compTime = $(commTime/compTime)")
#	end
	return true,commTime,compTime
end

function computeGradMisfit(sigma,
	                   DcRef::Array{Future,1},
	                   pMisRefs::Array{RemoteChannel,1},
	                   indFors=1:length(pMisRefs))
	#
	#	gc = computeGradMisfit(mc,model,Dc,Dobs,Wd,misfit,pFor)
	#
	#	compute gradient of misfit for multiple forward problems
	#
	#	Note: models and forward problems are represented as RemoteReferences.
	#
	#
	# find out which workers are involved
	workerList = []
	for k=1:length(pMisRefs)
		push!(workerList,pMisRefs[k].where)
	end
	workerList = unique(workerList)
	# send sigma to all workers
	sigmaRef = Array{RemoteChannel}(undef,maximum(workers()))
	dFiRef   = Array{RemoteChannel}(undef,maximum(workers()))
	dF = zeros(length(sigma))

	commTime = 0.0
	compTime = 0.0
	updateTimes(c1,c2) = (commTime+=c1; compTime+=c2)
	updateDF(x) = (dF+=x)

	@sync begin
		for p = workers()
				@async begin
					# send model to worker and get a remote ref
					 t = time_ns();
						sigmaRef[p] = initRemoteChannel(identity,p,sigma)   # send sigma to workers
						dFiRef[p] = initRemoteChannel(identity,p,zeros(length(sigma))) # get remote Ref to part of gradient
					c1 = (time_ns() - t)/1e+9;
					updateTimes(c1,0.0)

					# the actual computation
					for idx=indFors
						if pMisRefs[idx].where==p
							isDone,c1,c2 = remotecall_fetch(computeGradMisfit,p, sigmaRef[p],DcRef[idx],pMisRefs[idx],dFiRef[p])
							updateTimes(c1,c2)
						end
					end

					# fetch result and add
					c1 = @elapsed updateDF(fetch(dFiRef[p]))
					updateTimes(c1,0.0)
				end
		end
	end

	return dF
end

###############################################################################
#                                                                             #
# Method below is old, currently broken and possibly obsolete. However,       #
# it seems to include functionality for interrupting timed-out workers        #
# that is not present anywhere else in jInv. For that reason I keep it around #
# but commented out in case anyone wants to resurrect it in the future.       #
#                                                                             #
###############################################################################

# function computeGradMisfit(mc,Dc::Array,pMis::Array{Future},M2M::Array{Future},indFors=1:length(PF))
# 	#
# 	#	gc = computeGradMisfit(mc,model,Dc,Dobs,Wd,misfit,pFor)
# 	#
# 	#	compute gradient of misfit for multiple forward problems
# 	#
# 	#	Note: forward problems and interpolation matrices are represented as RemoteReferences.
# 	#
# 	#	!! there is only one(!) model. Interpolation matrix has to be entered for each problem !!
# 	#
# 	#	NEW: First worker that finishes takes the time and interrupt hanging workers.
# 	#
#
# 	dF     = zeros(length(mc))
# 	times  = fill(-1.0,length(PF))
#
# 	updateRes(dFi,tm,idx) = (dF+=dFi; times[idx]=tm)
#
# 	fwid = 0; tmStop = 0.0
# 	stopTime() = tmStop
# 	function isFirstWorker(p)
# 		if fwid ==0 # p is first worker that finished
# 			fwid = p # store p's id
# 			tmStop = time() + 3 * mean(times[times.>-1])  # set stop time
# 			return true
# 		else
# 			return false
# 		end
# 	end
#
# 	sigma,dsigmadm = mfun(mc)
#
# 	@sync begin
# 		for p = workers()
# 			@async begin
# 				for idx=indFors
# 					if PF[idx].where==p
# 						tic()
# 						dFi = remotecall_fetch(p, computeGradMisfit,sigma,gloc,Dc[idx],Dobs[idx],Wd[idx],misfit,PF[idx],M2M[idx])
# 						tm = toq()
# 						if dFi != -1; updateRes(dFi,tm,idx); else; break; end
# 					end
# 				end
#
# 				if isFirstWorker(p)
# 					while true
# 					  wt = remotecall_fetch(p,sleep,.01)
# 					  if minimum(times[indFors])>0.0 # all workers are done
# 						  break
# 					  elseif time()  > stopTime() # time to stop
# 							for pw=workers();
# 								if p!=pw; interrupt(pw) end;
# 							end
# 						  break
# 					  end
# 				  	end
# 				end
# 			end
# 		end
# 	end
# 	dF = (dsigmadm' * dF) / length(times[times.>-1])
# 	if length(times[times.>-1]) != length(indFors)
# 		println("computeGradMisfit was interrupted after computing ", length(times[times.>-1]) ," of ", length(indFors), " gradients")
# 	end
# 	return dF
# end

function computeGradMisfit(sigma,Dc::Array,pMis::Array{MisfitParam},indFors=1:length(pMis))
	#
	#	gc = computeGradMisfit(mc,model,Dc,Dobs,Wd,misfit,pFor)
	#
	#	compute gradient of misfit for multiple forward problems
	#
	#	Note: models, interpolations and forward problems are stored on main process and then sent to workers
	#
	#	!! this method may lead to more communication than the ones above !!
	#
	#
	numFor = length(pMis);

	# get process map
	i      = 1; nextidx = p -> i;
	procMap = zeros(Int64,numFor)

	###
	## Eran Treister:
	## The code below is a hack in case we're using a MUMPSsolver, who is no longer an integral part of jInv, due to dependency issues.
	## TODO: define the problem and fix it in a more general way for all solver with external data structures.
	## DO NOT DELETE THE CODE BELOW:
	## This appears also in HessMatVec.jl
	###
	# if isa(pMis[1].pFor.Ainv,MUMPSsolver) &&  (pMis[1].pFor.Ainv.Ainv.ptr !=-1)
		# for ii=1:numFor
			# if any(ii.==indFors)
				# procMap[ii] = pMis[ii].pFor.Ainv.Ainv.worker
			# end
		# end
		# nextidx = p -> ( ind = find(procMap.==p);
			# if !isempty(ind);
				# ind = ind[1];
				# procMap[ind] = -1;
			# end;
			# return ind)
	# else
		nextidx = p -> (idx=i; i+=1; idx)
	# end
	dF = zeros(length(sigma))
	updateRes(dFi) = (dF+=dFi)

	@sync begin
		for p = workers()
				@async begin
					while true
						idx = nextidx(p)
						if isempty(idx) || idx > numFor
							break
						end
						if any(idx.==indFors)
							dFi = remotecall_fetch(computeGradMisfit,p,sigma,Dc[idx],pMis[idx])
							updateRes(dFi)
						end
					end
				end
		end
	end
	return dF
end


text:
function solve()
    n = parse(Int, readline())
    as = [parse(Int, x) for x in split(readline())]

    ret = 0
    m = 0
    for x=2:1000
        s = sum([a%x==0 for a in as])
        if s > m
            m = s
            ret = x
        end
    end
    ret
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
cumtrapz(y, x, x0=0) = cumsum(vcat(0, diff(x).*(y[1:end-1] .+ y[2:end])/2)) .+ x0
invcumtrapz(y, x, x0=0) = cumsum(vcat(0, 2*diff(x)./(y[1:end-1] .+ y[2:end]))) .+ x0

@inline function distance2(line::Tuple{<:AbstractVector,<:AbstractVector}, x)
    v = line[2] - line[1]
    λ = clamp(v⋅(x - line[1])/(v⋅v), 0, 1)
    p = (1 - λ)*line[1] + λ*line[2]
    (p - x)⋅(p - x)
end

max_stars_repo_name: Michelle-NYX/CollisionAvoidance

max_issues_repo_name: Michelle-NYX/CollisionAvoidance

max_forks_repo_name: Michelle-NYX/CollisionAvoidance


text:
"""
    mutable struct Choi{dA, dB} <: QuantumVariable

An `AbstractVariable` which represents the Choi matrix of a quantum channel.

If `T` is a quantum channel from system A to system B, then its Choi matrix
is a `dA*dB` by `dA*dB` matrix given by `J` where

    J = T ⊗ id (Ω)

and `Ω =|Ω⟩⟨Ω|` where `|Ω⟩ = ∑ᵢ |i⟩|i⟩` is the unnormalized maximally
entangled state.
"""
mutable struct Choi{T, dA, dB} <: QuantumVariable{T}
    head::Symbol
    id_hash::UInt64
    size::Tuple{Int, Int}
    value::Convex.ValueOrNothing
    vexity::Convex.Vexity
    function Choi{T}(dA::Int, dB::Int) where {T <: Number}
        this = new{T, dA, dB}(:Choi, 0, (dA*dB,dA*dB), nothing, Convex.AffineVexity())
        this.id_hash = objectid(this)
        Convex.id_to_variables[this.id_hash] = this
        this
    end

    Choi{T}(d::Int) where {T <: Number} = Choi{T}(d, d)

    Choi(args...) = Choi{ComplexF64}(args...)
end

Convex.constraints(J::Choi{T, dA, dB}) where {T, dA, dB} = [ J ⪰ 0, Convex.partialtrace(J, 1, [dB, dA]) == I(dA) ]
Convex.sign(::Choi) = Convex.ComplexSign()
Convex.vartype(::Choi) = Convex.ContVar


text:
module BinTests

using Test

using PicoQuant
include("test_utils.jl")

@testset "Test executable files in PicoQuant/bin" begin
    gen_qft_circuit_executed = false
    qasm2tn_executed = false
    tn2plan_executed = false
    contract_tn_executed = false
    contract_qasm_executed = false

    n = 4 # number of qubits for qft circuit
    test_dir = "tmptest"
    if isdir(test_dir) rm(test_dir, recursive=true) end
    mkdir("tmptest")

    # Try executing the scripts
    try
        test_ARGS = ["--number-qubits", "$n", "-o", joinpath(test_dir, "qft_$n.qasm")]
        include("../bin/gen_qft_circuit.jl")
        main(test_ARGS)
        gen_qft_circuit_executed = isfile(joinpath(test_dir, "qft_$n.qasm"))

        test_ARGS = ["--qasm", joinpath(test_dir, "qft_$n.qasm")]
        include("../bin/qasm2tn.jl")
        main(test_ARGS)
        qasm2tn_executed = isfile(joinpath(test_dir, "qft_$n.json"))  &&
                           isfile(joinpath(test_dir, "qft_$n.tl")) &&
                           isfile(joinpath(test_dir, "qft_$n.h5"))

        test_ARGS = ["--tng", joinpath(test_dir, "qft_$n.json"),
                     "--tng_data", joinpath(test_dir, "qft_$n.h5"),
                     "--dsl", joinpath(test_dir, "qft_$(n)_contract"),
                     "--output_tng", joinpath(test_dir, "qft_$(n)_contract.json")]
        include("../bin/tn2dsl.jl")
        main(test_ARGS)
        tn2plan_executed = isfile(joinpath(test_dir, "qft_$(n)_contract.tl")) &&
                           isfile(joinpath(test_dir, "qft_$(n)_contract.h5")) &&
                           isfile(joinpath(test_dir, "qft_$(n)_contract.json"))

        test_ARGS = ["--dsl", joinpath(test_dir, "qft_$(n)_contract.tl"),
                     "--dsl_data", joinpath(test_dir, "qft_$(n)_contract.h5"),
                     "--output", joinpath(test_dir, "qft_$(n)_contract_out.h5")]
        include("../bin/contract_tn.jl")
        main(test_ARGS)
        contract_tn_executed = isfile(joinpath(test_dir, "qft_$(n)_contract_out.h5"))

        test_ARGS = ["--qasm", joinpath(test_dir, "qft_$n.qasm"),
                     "--output", joinpath(test_dir, "output.json"),
                     "--output_data", joinpath(test_dir, "output.h5")]
        include("../bin/contract_qasm.jl")
        main(test_ARGS)
        contract_qasm_executed = isfile(joinpath(test_dir, "output.json")) &&
                                 isfile(joinpath(test_dir, "output.h5"))

    finally
        # clean up any created files
        for fn in ["qft_$n.qasm",
                   "qft_$n.json",
                   "qft_$n.tl",
                   "qft_$n.h5",
                   "qft_$(n)_contract.tl",
                   "qft_$(n)_contract.h5",
                   "qft_$(n)_contract.json",
                   "qft_$(n)_contract_out.h5",
                   "output.json",
                   "output.h5"]
            rm(joinpath(test_dir, fn), force=true)
        end
        rm(test_dir, recursive=true)
    end

    # Test if scripts executed
    @test gen_qft_circuit_executed == true
    @test qasm2tn_executed == true
    @test tn2plan_executed == true
    @test contract_tn_executed == true
    @test contract_qasm_executed == true
end
end


text:
module Utils

using Recyclers

function recycler_types()
    ts = Any[Recyclers.ShardedRecycler]
    t = try
        Recyclers.CentralizedRecycler
    catch
        nothing
    end
    t === nothing || push!(ts, t)
    return ts
end

function cache_types()
    ts = Any[Recyclers.ShardedCache]
    t = try
        Recyclers.CentralizedCache
    catch
        nothing
    end
    t === nothing || push!(ts, t)
    return ts
end

end  # module


text:
using Test
using Printf
using FileIO

# SortedOffsets depends on this
include("../../src/image.jl")
# include("../src/synthResult.jl")

# the tested code
include("../../src/sortedOffsets.jl")

include("tinyImage.jl")

@testset "SortedOffsets" begin

    testImage = createTestMaskedImage()

    testSortedOffsets = SortedOffsets(testImage.image)

    @debug "Sorted offsets" testSortedOffsets

    # offset 0,0 is not in
    @test ! any( offset == CartesianIndex(0,0) for offset in testSortedOffsets.offsets )

    # TODO add .offsets below

    # lower right of image is last element of sortedOffsets
    # Is  this a robust test?  Other large distance offset could be in last element?
    dimensions = size(testImage.image)
    @test last(testSortedOffsets) == CartesianIndex(dimensions)

    # sortedOffsets are unique
    @test length(unique(testSortedOffsets)) == length(testSortedOffsets)

    # sortedOffsets spans
    # approximately four times
    # TODO exact is (2x+1)(2y+1) == 4xy + 2x + 2y + 1
    @test length(testSortedOffsets) >=  4 * length(testImage.image)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
function PlainChain( dims...; σ = identity )
    Flux.Chain( collect( Flux.Dense(i,f,σ) for (i,f) in zip( dims[1:(end-1)], dims[2:end] ) )...  )
end

function TChain( dims...; σ = identity )
    Flux.Chain( collect( TenseLayer(i,f,σ) for (i,f) in zip( dims[1:(end-1)], dims[2:end] ) )...  )
end


text:
"""
`prepare_interpolation(nodes::Matrix{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare the interpolating normal spline by constructing and factoring Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space and
           `n_1` is the number of function value nodes. It means that each column in the matrix defines one node.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_interpolation(nodes::Matrix{T},
                               kernel::RK = RK_H0()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_interpolation(nodes, kernel)
     return spline
end

"""
`construct_interpolation(spline::NormalSpline{T, RK}, values::Vector{T}) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Construct the interpolating normal spline by calculating its coefficients and
completely initializing the `NormalSpline` object.
# Arguments
- `spline`: the prepared `NormalSpline` object.
- `values`: function values at interpolation nodes.

Return: constructed `NormalSpline` object.
"""
function construct_interpolation(spline::NormalSpline{T, RK},
                                 values::Vector{T}
                                ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
    spline = _construct_interpolation(spline, values)
    return spline
end

"""
`interpolate(nodes::Matrix{T}, values::Vector{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare and construct the interpolating normal spline.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space
           and `n_1` is the number of function value nodes.
           It means that each column in the matrix defines one node.
- `values`: function values at interpolation nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.
Return: constructed `NormalSpline` object.
"""
function interpolate(nodes::Matrix{T},
                     values::Vector{T},
                     kernel::RK = RK_H0()
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_interpolation(nodes, kernel)
     spline = _construct_interpolation(spline, values)
     return spline
end

########

"""
`prepare_interpolation(nodes::Matrix{T}, d_nodes::Matrix{T}, es::Matrix{T}, kernel::RK = RK_H1())
                       where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Prepare the spline by constructing and factoring a Gram matrix of the interpolation problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space and
           `n_1` is the number of function value nodes.
            It means that each column in the matrix defines one node.
- `d_nodes`: function directional derivatives nodes.
             This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
             `n_2` is the number of function directional derivative nodes.
- `es`: Directions of the function directional derivatives.
        This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
        `n_2` is the number of function directional derivative nodes.
        It means that each column in the matrix defines one direction of the function directional derivative.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_interpolation(nodes::Matrix{T},
                               d_nodes::Matrix{T},
                               es::Matrix{T},
                               kernel::RK = RK_H1()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
     spline = _prepare_interpolation(nodes, d_nodes, es, kernel)
     return spline
end

"""
`construct_interpolation(spline::NormalSpline{T, RK}, values::Vector{T}, d_values::Vector{T})
                         where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Construct the interpolating normal spline by calculating its coefficients and
completely initializing the `NormalSpline` object.
# Arguments
- `spline`: the prepared `NormalSpline` object.
- `values`: function values at interpolation nodes.
- `d_values`: function directional derivative values at function derivative nodes.

Return: constructed `NormalSpline` object.
"""
function construct_interpolation(spline::NormalSpline{T, RK},
                                 values::Vector{T},
                                 d_values::Vector{T}
                                ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
     spline = _construct_interpolation(spline, values, d_values)
     return spline
end

"""
`interpolate(nodes::Matrix{T}, values::Vector{T}, d_nodes::Matrix{T}, es::Matrix{T}, d_values::Vector{T}, kernel::RK = RK_H1())
             where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Prepare and construct the interpolating normal spline.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n×n_1` matrix, where `n` is dimension of the sampled space
           and `n_1` is the number of function value nodes.
           It means that each column in the matrix defines one node.
- `values`: function values at interpolation nodes.
- `d_nodes`: function directional derivative nodes.
            This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
            `n_2` is the number of function directional derivative nodes.
- `es`: Directions of the function directional derivatives.
       This should be an `n×n_2` matrix, where `n` is dimension of the sampled space and
       `n_2` is the number of function directional derivative nodes.
       It means that each column in the matrix defines one direction of the function directional derivative.
- `d_values`: function directional derivative values at function derivative nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: constructed `NormalSpline` object.
"""
function interpolate(nodes::Matrix{T},
                     values::Vector{T},
                     d_nodes::Matrix{T},
                     es::Matrix{T},
                     d_values::Vector{T},
                     kernel::RK = RK_H1()
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
     spline = _prepare_interpolation(nodes, d_nodes, es, kernel)
     spline = _construct_interpolation(spline, values, d_values)
     return spline
end

############################## One-dimensional case

"""
`prepare_interpolation(nodes::Vector{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare the 1D interpolating normal spline by constructing and factoring a Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n_1` vector where `n_1` is the number of function value nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_interpolation(nodes::Vector{T},
                               kernel::RK = RK_H0()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_interpolation(Matrix(nodes'), kernel)
     return spline
end

"""
`interpolate(nodes::Vector{T}, values::Vector{T}, kernel::RK = RK_H0()) where {T <: AbstractFloat, RK <: ReproducingKernel_0}`

Prepare and construct the 1D interpolating normal spline.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n_1` vector where `n_1` is the number of function value nodes.
- `values`: function values at `n_1` interpolation nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H0` if the spline is constructing as a continuous function,
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: constructed `NormalSpline` object.
"""
function interpolate(nodes::Vector{T},
                     values::Vector{T},
                     kernel::RK = RK_H0()
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_0}
     spline = _prepare_interpolation(Matrix(nodes'), kernel)
     spline = _construct_interpolation(spline, values)
     return spline
end

"""
`prepare_interpolation(nodes::Vector{T}, d_nodes::Vector{T}, kernel::RK = RK_H1()) where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

Prepare the 1D interpolating normal spline by constructing and factoring a Gram matrix of the problem.
Initialize the `NormalSpline` object.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n_1` vector where `n_1` is the number of function value nodes.
- `d_nodes`: The function derivatives nodes.
             This should be an `n_2` vector where `n_2` is the number of function derivatives nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: prepared `NormalSpline` object.
"""
function prepare_interpolation(nodes::Vector{T},
                               d_nodes::Vector{T},
                               kernel::RK = RK_H1()
                              ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}

     es = ones(T, length(d_nodes))
     spline = _prepare_interpolation(Matrix(nodes'), Matrix(d_nodes'), Matrix(es'), kernel)
     return spline
end

"""
`interpolate(nodes::Vector{T}, values::Vector{T}, d_nodes::Vector{T}, d_values::Vector{T},
             kernel::RK = RK_H1())
             where {T <: AbstractFloat, RK <: ReproducingKernel_1}`

prepare_interpolation and construct the 1D interpolating normal spline.
# Arguments
- `nodes`: function value interpolation nodes.
           This should be an `n_1` vector where `n_1` is the number of function value nodes.
- `values`: function values at `nodes` nodes.
- `d_nodes`: The function derivatives nodes.
             This should be an `n_2` vector where `n_2` is the number of function derivatives nodes.
- `d_values`: function derivative values at function derivative nodes.
- `kernel`: reproducing kernel of Bessel potential space the normal spline is constructed in.
            It must be a struct object of the following type:
              `RK_H1` if the spline is constructing as a differentiable function,
              `RK_H2` if the spline is constructing as a twice differentiable function.

Return: constructed `NormalSpline` object.
"""
function interpolate(nodes::Vector{T},
                     values::Vector{T},
                     d_nodes::Vector{T},
                     d_values::Vector{T},
                     kernel::RK = RK_H1()
                    ) where {T <: AbstractFloat, RK <: ReproducingKernel_1}
     es = ones(T, length(d_nodes))
     spline = _prepare_interpolation(Matrix(nodes'), Matrix(d_nodes'), Matrix(es'), kernel)
     spline = _construct_interpolation(spline, values, d_values)
     return spline
end