text:
"""
    tf!(dtm::AbstractMatrix{Real}, tf::AbstractMatrix{AbstractFloat})

Overwrite `tf` with the term frequency of the `dtm`.

Works correctly if `dtm` and `tf` are same matrix.

See also: [`tf`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
function tf!(dtm::AbstractMatrix{T1}, tf::AbstractMatrix{T2}) where {T1 <: Real, T2 <: AbstractFloat}
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    for i in 1:n
        words_in_document = 0
        for j in 1:p
            words_in_document += dtm[i, j]
        end
        tf[i, :] = dtm[i, :] ./ max(words_in_document, one(T1))
    end

    return tf
end

"""
    tf!(dtm::SparseMatrixCSC{Real}, tf::SparseMatrixCSC{AbstractFloat})

Overwrite `tf` with the term frequency of the `dtm`.

`tf` should have the has same nonzeros as `dtm`.

See also: [`tf`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
function tf!(dtm::SparseMatrixCSC{T}, tf::SparseMatrixCSC{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    tfvals = nonzeros(tf)
    @assert size(dtmvals) == size(tfvals)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = sum(dtm,dims=2)

    n, p = size(dtm)
    for i = 1:p
       for j in nzrange(dtm, i)
          row = rows[j]
          tfvals[j] = dtmvals[j] / max(words_in_documents[row], one(T))
       end
    end
    return tf
end

tf!(dtm::AbstractMatrix{T}) where {T <: Real} = tf!(dtm, dtm)

tf!(dtm::SparseMatrixCSC{T}) where {T <: Real} = tf!(dtm, dtm)

"""
    tf(dtm::DocumentTermMatrix)
    tf(dtm::SparseMatrixCSC{Real})
    tf(dtm::Matrix{Real})

Compute the `term-frequency` of the input.

# Example

```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
              StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> m = DocumentTermMatrix(crps)

julia> tf(m)
2×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 10 stored entries:
  [1, 1]  =  0.166667
  [2, 1]  =  0.166667
  [1, 2]  =  0.333333
  [2, 3]  =  0.333333
  [1, 4]  =  0.166667
  [2, 4]  =  0.166667
  [1, 5]  =  0.166667
  [2, 5]  =  0.166667
  [1, 6]  =  0.166667
  [2, 6]  =  0.166667
```

See also: [`tf!`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
tf(dtm::DocumentTermMatrix) = tf(dtm.dtm)

tf(dtm::Matrix{T}) where {T <: Real} = tf!(dtm, Array{Float64}(undef, size(dtm)...))

tf(dtm::SparseMatrixCSC{T}) where {T <: Real} =  tf!(dtm, similar(dtm, Float64))

"""
    tf_idf!(dtm::AbstractMatrix{Real}, tf_idf::AbstractMatrix{AbstractFloat})

Overwrite `tf_idf` with the tf-idf (Term Frequency - Inverse Doc Frequency) of the `dtm`.

`dtm` and `tf-idf` must be matrices of same dimensions.

See also: [`tf`](@ref), [`tf!`](@ref) , [`tf_idf`](@ref)
"""
function tf_idf!(dtm::AbstractMatrix{T1}, tfidf::AbstractMatrix{T2}) where {T1 <: Real, T2 <: AbstractFloat}
    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    tf!(dtm, tfidf)

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=1))
    idf = log.(n ./ documents_containing_term)

    # TF-IDF is the product of TF and IDF
    for i in 1:n
        for j in 1:p
           tfidf[i, j] = tfidf[i, j] * idf[j]
        end
    end

    return tfidf
end

"""
    tf_idf!(dtm::SparseMatrixCSC{Real}, tfidf::SparseMatrixCSC{AbstractFloat})

Overwrite `tfidf` with the tf-idf (Term Frequency - Inverse Doc Frequency) of the `dtm`.

The arguments must have same number of nonzeros.

See also: [`tf`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
function tf_idf!(dtm::SparseMatrixCSC{T}, tfidf::SparseMatrixCSC{F}) where {T <: Real, F <: AbstractFloat}
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    tfidfvals = nonzeros(tfidf)
    @assert size(dtmvals) == size(tfidfvals)

    n, p = size(dtm)

    # TF tells us what proportion of a document is defined by a term
    words_in_documents = F.(sum(dtm, dims=2))
    oneval = one(F)

    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=1))
    idf = log.(n ./ documents_containing_term)

    for i = 1:p
       for j in nzrange(dtm, i)
          row = rows[j]
          tfidfvals[j] = dtmvals[j] / max(words_in_documents[row], oneval) * idf[i]
       end
    end

    return tfidf
end

"""
    tf_idf!(dtm)

Compute tf-idf for `dtm`
"""
tf_idf!(dtm::AbstractMatrix{T}) where {T <: Real} = tf_idf!(dtm, dtm)

tf_idf!(dtm::SparseMatrixCSC{T}) where {T <: Real} = tf_idf!(dtm, dtm)

# This does not make sense, since DocumentTermMatrix is based on an array of integers
#tf_idf!(dtm::DocumentTermMatrix) = tf_idf!(dtm.dtm)

"""
    tf(dtm::DocumentTermMatrix)
    tf(dtm::SparseMatrixCSC{Real})
    tf(dtm::Matrix{Real})

Compute `tf-idf` value (Term Frequency - Inverse Document Frequency) for the input.

In many cases, raw word counts are not appropriate for use because:

- Some documents are longer than other documents
- Some words are more frequent than other words

A simple workaround this can be done by performing `TF-IDF` on a `DocumentTermMatrix`

# Example

```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
              StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> m = DocumentTermMatrix(crps)

julia> tf_idf(m)
2×6 SparseArrays.SparseMatrixCSC{Float64,Int64} with 10 stored entries:
  [1, 1]  =  0.0
  [2, 1]  =  0.0
  [1, 2]  =  0.231049
  [2, 3]  =  0.231049
  [1, 4]  =  0.0
  [2, 4]  =  0.0
  [1, 5]  =  0.0
  [2, 5]  =  0.0
  [1, 6]  =  0.0
  [2, 6]  =  0.0
```

See also: [`tf!`](@ref), [`tf_idf`](@ref), [`tf_idf!`](@ref)
"""
tf_idf(dtm::DocumentTermMatrix) = tf_idf(dtm.dtm)

tf_idf(dtm::SparseMatrixCSC{T}) where {T <: Real} =  tf_idf!(dtm, similar(dtm, Float64))

tf_idf(dtm::Matrix{T}) where {T <: Real} = tf_idf!(dtm, Array{Float64}(undef, size(dtm)...))

function bm_25!(dtm::AbstractMatrix{T},
                bm25::AbstractMatrix{F};
                κ::Int=2,
                β::Float64=0.75
               ) where {T<:Real, F<:AbstractFloat}
    @assert size(dtm) == size(bm25)
    # Initializations
    k = F(κ)
    b = F(β)
    p, n = size(dtm)
    oneval = one(F)
    # TF tells us what proportion of a document is defined by a term
    words_in_documents = F.(sum(dtm, dims=1))
    ln = words_in_documents./mean(words_in_documents)
    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=2)) .+ one(T)
    idf = log.(n ./ documents_containing_term) .+ oneval
    # BM25 is the product of IDF and a fudged TF
    tf_bm25!(dtm, bm25)
    @inbounds @simd for i in 1:n
        for j in 1:p
            bm25[j, i] = idf[j] *
                ((k + 1) * bm25[j, i]) /
                (k * (oneval - b + b * ln[i]) + bm25[j, i])
        end
    end
    return bm25
end

function bm_25!(dtm::SparseMatrixCSC{T},
                bm25::SparseMatrixCSC{F};
                κ::Int=2,
                β::Float64=0.75
               ) where {T<:Real, F<:AbstractFloat}
    @assert size(dtm) == size(bm25)
    # Initializations
    k = F(κ)
    b = F(β)
    rows = rowvals(dtm)
    dtmvals = nonzeros(dtm)
    bm25vals = nonzeros(bm25)
    @assert size(dtmvals) == size(bm25vals)
    p, n = size(dtm)
    # TF tells us what proportion of a document is defined by a term
    words_in_documents = F.(sum(dtm, dims=1))
    ln = words_in_documents./mean(words_in_documents)
    oneval = one(F)
    # IDF tells us how rare a term is in the corpus
    documents_containing_term = vec(sum(dtm .> 0, dims=2)) .+ one(T)
    idf = log.(n ./ documents_containing_term) .+ oneval
    for i = 1:n
       for j in nzrange(dtm, i)
          row = rows[j]
          tf = sqrt.(dtmvals[j] / max(words_in_documents[i], oneval))
          bm25vals[j] = idf[row] * ((k + 1) * tf) /
                        (k * (oneval - b + b * ln[i]) + tf)
       end
    end
    return bm25
end

bm_25(dtm::AbstractMatrix{T}; κ::Int=2, β::Float64=0.75) where T<:Integer =
    bm_25!(dtm, similar(dtm, Float64), κ=κ, β=β)

bm_25(dtm::AbstractMatrix{T}; κ::Int=2, β::Float64=0.75) where T<:AbstractFloat =
    bm_25!(dtm, similar(dtm, T), κ=κ, β=β)

bm_25(dtm::DocumentTermMatrix; κ::Int=2, β::Float64=0.75) =
    bm_25(dtm.dtm, κ=κ, β=β)

bm_25!(dtm::DocumentTermMatrix; κ::Int=2, β::Float64=0.75) =
    bm_25!(dtm.dtm, κ=κ, β=β)

# The score was modified according to for bm25:
#   https://opensourceconnections.com/blog/2015/10/16/bm25-the-next-generation-of-lucene-relevation/
function tf_bm25!(dtm::AbstractMatrix{T}, tf::AbstractMatrix{F}
            ) where {T<:Real, F<:AbstractFloat}
    @assert size(dtm) == size(tf)
    p, n = size(dtm)
    # TF tells us what proportion of a document is defined by a term
    words_in_documents = sum(dtm, dims=1)
    @inbounds for i in 1:n
        tf[:, i] = sqrt.(dtm[:, i] ./ max(words_in_documents[i], one(T)))
    end
    return tf
end


text:
using Fermi.Options
import TBLIS
import TensorOperations: contract!

@testset "Tensor Contraction" begin
    Fermi.tblis_set_num_threads(1)
    @test TBLIS.get_num_threads() == 1

    # Using Aᵢₘₙ⋅Bₘₙⱼ = Cᵢⱼ
    oindA = (1,)
    cindA = (2,3)
    oindB = (3,)
    cindB = (1,2)
    tindC = (1,2)
    x = Fermi.oind2eins(oindA, cindA, oindB, cindB, tindC)
    @test x == ("Bbc", "bcC", "BC")

    A = FermiMDrand(2,5,5)
    B = FermiMDrand(5,5,2)
    C1 = FermiMDzeros(2,2)
    C2 = FermiMDzeros(2,2)

    @set tblis true
    α = rand()
    β = rand()
    contract!(α, A, :N, B, :N, β, C1, oindA, cindA, oindB, cindB, tindC)
    contract!(α, A.data, :N, B.data, :N, β, C2.data, oindA, cindA, oindB, cindB, tindC)
    @test C1 ≈ C2

    @test_throws FermiException contract!(α, A, :C, B, :N, β, C1, oindA, cindA, oindB, cindB, tindC)

    @set tblis false
    C2 = FermiMDzeros(2,2)
    contract!(α, A, :N, B, :N, β, C2, oindA, cindA, oindB, cindB, tindC)
    @test C1 ≈ C2
end

ext: jl


text:
function test_connection()
  println("Using libpq")
    libpq = PostgreSQL.libpq_interface

    println("Checking basic connect")
    conn = connect(Postgres, "localhost", "postgres", "", "julia_test")
    @test isa(conn, PostgreSQL.PostgresDatabaseHandle)
    @test conn.status == PostgreSQL.CONNECTION_OK
    @test errcode(conn) == PostgreSQL.CONNECTION_OK
    @test !conn.closed
    @test bytestring(libpq.PQdb(conn.ptr)) == "julia_test"
    @test bytestring(libpq.PQuser(conn.ptr)) == "postgres"
    @test bytestring(libpq.PQport(conn.ptr)) == "5432"

    disconnect(conn)
    @test conn.closed
    println("Basic connection passed")

    println("Checking doblock")
    conn = connect(Postgres, "localhost", "postgres", "", "julia_test") do conn
        @test isa(conn, PostgreSQL.PostgresDatabaseHandle)
        @test conn.status == PostgreSQL.CONNECTION_OK
        @test errcode(conn) == PostgreSQL.CONNECTION_OK
        @test !conn.closed
        return conn
    end
    @test conn.closed
    println("Doblock passed")

    println("Testing connection with DSN string")
    conn = connect(Postgres; dsn="postgresql://postgres@localhost:5432/julia_test")
    @test isa(conn, PostgreSQL.PostgresDatabaseHandle)
    @test conn.status == PostgreSQL.CONNECTION_OK
    @test errcode(conn) == PostgreSQL.CONNECTION_OK
    @test !conn.closed
    @test bytestring(libpq.PQdb(conn.ptr)) == "julia_test"
    @test bytestring(libpq.PQuser(conn.ptr)) == "postgres"
    @test bytestring(libpq.PQport(conn.ptr)) == "5432"

    disconnect(conn)
    @test conn.closed
    println("DSN connection passed")

#=    println("Testing connection with DSN string and doblock")
    conn = connect(Postgres; dsn="postgresql://postgres@localhost:5432/postgres") do conn
        @test isa(conn, PostgreSQL.PostgresDatabaseHandle)
        @test conn.status == PostgreSQL.CONNECTION_OK
        @test errcode(conn) == PostgreSQL.CONNECTION_OK
        @test !conn.closed
        return conn
    end
    @test conn.closed
    println("DSN connection passed in a do block")=#
end

test_connection()

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
using LightGraphs, SimpleWeightedGraphs
export Day7, Day7_1, Day7_2
abstract type Day7 <: Day end
abstract type Day7_1 <: Day7 end
abstract type Day7_2 <: Day7 end

struct Bag
    hue
    color
    Bag(hue, color) = new(Symbol(hue), Symbol(color))
    Bag(str::AbstractString) = str |> split |> p -> new(Symbol(p[1]), Symbol(p[2]))
end

struct BagGraph
    bag_graph::Dict{Bag,Int}
    graph_bag::Dict{Int,Bag}
    graph::SimpleWeightedDiGraph
    BagGraph(lines::Vector) = begin
        g = SimpleWeightedDiGraph()
        bag_graph = Dict{Bag,Int}()
        graph_bag = Dict{Int,Bag}()

        node_count = 0
        for (b, bag_map) in lines
            if b ∉ keys(bag_graph)
                node_count += 1
                bag_graph[b] = node_count
                graph_bag[node_count] = b
                add_vertex!(g)
            end
            for (b_to, v) in bag_map
                if b_to ∉ keys(bag_graph)
                    node_count += 1
                    bag_graph[b_to] = node_count
                    graph_bag[node_count] = b_to
                    add_vertex!(g)
                end
                add_edge!(g, bag_graph[b], bag_graph[b_to], v)
            end
        end
        if !is_directed(g)
            throw("Not Directed")
        end

        new(bag_graph, graph_bag, g)
    end
end

line_parser(::Type{T}, line) where T<: Day7 = begin
    line = line |> l->replace(l,Pair(".","")) |> strip
    bag_first, rest = split(line, " contain ")
    node = Bag(bag_first)

    can_contain = split(rest, ", ")
    contained = Dict{Bag,Int}()
    for b in can_contain
        if startswith(b, "no")
            continue
        else
            count, bag_next = split(b; limit=2)
            contained[Bag(bag_next)] = parse(Int, count)
        end
    end
    [node, contained]
end

function parents(s,g)
    r = neighbors(g,s)
    if length(r) == 0
        return Set([s])
    else
        res = [parents(n,g) for n in r]
        return union(Set([s]), res...)
    end
end

function bags_contained(s,g, w)
    r = neighbors(g,s)
    if length(r) == 0
        return 0
    else
        sub_bags = [(x,bags_contained(x,g, w), max(w[x,s],w[s,x])) for x in r]
        tot = sum([(1 + n)*v for (x,n,v) in sub_bags])
        return tot
    end
end

solve(::Type{T}, lines::Vector) where T <: Day7 = lines |> BagGraph |> bg -> solve(T, bg)

solve(::Type{T}, bg::BagGraph) where T <: Day7_1 = bg |> bg->parents(bg.bag_graph[Bag("shiny","gold")],DiGraph(bg.graph)) |> p->length(p) - 1
solve(::Type{T}, bg::BagGraph) where T <: Day7_2 = bg |>  bg->bags_contained(bg.bag_graph[Bag("shiny","gold")],reverse(DiGraph(bg.graph)),weights(bg.graph)) |> Int

tree_root(hue,color,dir, bg::BagGraph) = dfs_tree(bg.graph, bg.bag_graph[Bag(hue, color)];dir=dir)


text:
include("../../__lib__/std/util/src/util_module.jl")
include("../../__lib__/std/collections/gridded/src/gridded_module.jl")

module compgraph
using ..util

function contour_level(xs, ys, z, level)
    #   vertex    edges
	##  1-----2   .-1-.
    ##  |     |   |   |
    ##  |     |   4   2
    ##  |     |   |   |
    ##  4-----3   .-3-.

    Z = z .- level
    mask = Z .>= 0
    mark = mark_cells(mask)

    # show(stdout, MIME("text/plain"), mark)

    accum_array = zeros(2, 0)

    (n_row, n_col) = size(mark)

    # top search
    step_dir   = down::StepDirections
    first           = true
    start_edge      = 1
    rr              = 1
    for cc = 1:n_col
        id = mark[rr, cc]
        if id > 0 && (  id == 2 || id == 6 || id == 14 ||  id == 1 || id == 9 || id == 13 )
            (mark, accum_array) = __search_path!(xs, ys, Z, mark, accum_array, rr, cc, start_edge, first, step_dir, level)
        end
    end

    # bottom search
    step_dir   = up::StepDirections
    first           = true
    start_edge      = 3
    rr              = n_row
    for cc = 1:n_col
        id = mark[rr, cc]
        if id > 0 && ( id == 8 || id == 9 || id == 11 || id == 4 || id == 6 || id == 7)
            (mark, accum_array) = __search_path!(xs, ys, Z, mark, accum_array, rr, cc, start_edge, first, step_dir, level)
        end
    end

    # left search
    step_dir   = right::StepDirections
    first           = true
    start_edge      = 4
    cc              = 1
    for rr = 1:n_row
        id = mark[rr, cc]
        if id > 0 && ( id == 1 || id == 3 || id == 7 || id == 8 || id == 12 || id == 14)
            (mark, accum_array) = __search_path!(xs, ys, Z, mark, accum_array, rr, cc, start_edge, first, step_dir, level)
        end
    end

    # right search
    step_dir   = left::StepDirections
    first           = true
    start_edge      = 2
    cc              = n_col
    for rr = 1:n_row
        id = mark[rr, cc]
        if id > 0 &&  (id == 4 || id == 12 || id == 13 || id == 2 || id == 3 || id == 11)
            (mark, accum_array) = __search_path!(xs, ys, Z, mark, accum_array, rr, cc, start_edge, first, step_dir, level)
        end
    end

    # interior search
    step_dir = none::StepDirections
    first = true
    start_edge = 1
    for cc = 1:n_col-1
        for rr = 1:n_row-1
            id = mark[rr, cc]
            if id > 0
                (mark, accum_array) = __search_path!(xs, ys, Z, mark, accum_array, rr, cc, start_edge, first, step_dir, level)
            end
        end
    end

    return (accum_array, mark)
end

@enum StepDirections begin
    none
    up
    right
    down
    left
end

function __search_path!(xs, ys, Z, mark, accum_array, r, c, init_start_edge, first, step_dir::StepDirections, level)

    (n_row, n_col) = size(mark)
    start_idx = size(accum_array, 2) + 1

    # step_dir = none::StepDirections
    prev_step_dir = none::StepDirections # overrite for now

    stop_edge = Int8(127)
    while r >= 1 && c >= 1 && r <= n_row && c <= n_col && mark[r, c] > 0

        px = ( xs[c], xs[c+1], xs[c+1], xs[c] )
        py = ( ys[r], ys[r],   ys[r+1], ys[r+1] )
        pz = ( Z[r, c], Z[r, c+1], Z[r+1, c+1], Z[r+1, c] )

        # find start edge
        (start_edge, stop_edge) = find_start_and_stop_edge(mark, r, c)
        # if first
        #     start_edge = init_start_edge
        # end

        #   vertex    edges
        ##  1-----2   .-1-.
        ##  |     |   |   |
        ##  |     |   4   2
        ##  |     |   |   |
        ##  4-----3   .-3-.

        # switch the start and stop edges
        if step_dir == none::StepDirections
            # do nothing
        elseif step_dir == up::StepDirections
            if start_edge != 3
                start_edge, stop_edge = stop_edge, start_edge
            end
        elseif step_dir == down::StepDirections
            if start_edge != 1
                start_edge, stop_edge = stop_edge, start_edge
            end
        elseif step_dir == left::StepDirections
            if start_edge != 2
                start_edge, stop_edge = stop_edge, start_edge
            end
        elseif step_dir == right::StepDirections
            if start_edge != 4
                start_edge, stop_edge = stop_edge, start_edge
            end
        end

        if start_edge == 127
            # unsuccessful search, break up
            break
        end

        if first

            pt = (start_edge, util.modulo_index(start_edge+1,4))
            (ct_x, ct_y) = lerp_interp(px, py, pz, pt)

            accum_array = add_point!(accum_array, level, 0)
            accum_array = add_point!(accum_array, ct_x, ct_y)
            first = false
        end

        pt = (stop_edge, util.modulo_index(stop_edge+1,4))
        (ct_x, ct_y) = lerp_interp(px, py, pz, pt)
        accum_array = add_point!(accum_array, ct_x, ct_y)

        mark[r, c] -= mark[r, c]

        prev_step_dir = step_dir
        if (stop_edge == 1)
            r -= 1;
            step_dir = up::StepDirections;
        elseif (stop_edge == 2)
            c += 1;
            step_dir = right::StepDirections;
        elseif (stop_edge == 3)
            r += 1;
            step_dir = down::StepDirections;
        elseif (stop_edge == 4)
            c -= 1;
            step_dir = left::StepDirections;
        end

    end # while

    end_idx = size(accum_array, 2)
    if start_idx >= 1
        accum_array[2,start_idx] = end_idx - start_idx
    end

    return (mark, accum_array)
end

function add_point!(accum_array, x, y)
    accum_array = hcat(accum_array, [x, y])
    return accum_array
end

function mark_cells(mask)
    (n_row, n_col) = size(mask)

    mark = zeros(Int8, n_row-1, n_col-1)

    for jj = 1:n_col-1
        for ii = 1:n_row-1

            f1 = mask[ii, jj]
            f2 = mask[ii, jj+1]
            f3 = mask[ii+1, jj+1]
            f4 = mask[ii+1, jj]

            mark[ii, jj] = f1 + f2 * 2 + f3 * 4 + f4 * 8

        end
    end

    mark[ mark .>= 15] .= 0

    return mark

end

function lerp_interp(px, py, pz, pt)

    tmp = pz[pt[2]] / pz[pt[1]];
    tmp = abs(tmp);

    if !isnan(tmp)
        lamda = 1 / (1 + tmp)
        ct_x = px[pt[1]] * (1-lamda) + px[pt[2]] * lamda
        ct_y = py[pt[1]] * (1-lamda) + py[pt[2]] * lamda
    else
        lamda = 0.5;
        ct_x = px[pt[1]] * (1-lamda) + px[pt[2]] * lamda
        ct_y = py[pt[1]] * (1-lamda) + py[pt[2]] * lamda
    end

    return (ct_x, ct_y)
end

function find_start_and_stop_edge(mark, r, c)
    id = mark[r, c]

    start_edge = Int8(127)
    stop_edge = Int8(127)

    if id == 0 || id == 15
        # impossible condition
        return (start_edge, stop_edge)
    end

    if id == 1 || id == 3 || id == 7
        start_edge = 4
    elseif id == 2 || id == 6 || id == 14
        start_edge = 1
    elseif id == 4 || id == 12 || id == 13
        start_edge = 2
    elseif id == 8 || id == 9 || id == 11
        start_edge = 3
    end

    if id == 1 || id == 9 || id == 13
        stop_edge = 1
    elseif id == 2 || id == 3 || id == 11
        stop_edge = 2
    elseif id == 4 || id == 6 || id == 7
        stop_edge = 3
    elseif id == 8 || id == 12 || id == 14
        stop_edge = 4
    end

    if id == 5 || id == 10
        start_edge = -1
        stop_edge = -1
    end

    start_edge = Int8(start_edge)
    stop_edge = Int8(stop_edge)

    return (start_edge, stop_edge)

end

end

module ctest
using ..util
using ..compgraph
using ..gridded
using PyPlot
PyPlot.pygui(true)

Nx = 10 # 10;
Ny = 8 # 8;

x_min = -1.0;
x_max = 2.0;

y_min = -2.0;
y_max = 3.0;

xs = collect( range(x_min, x_max, Nx) );
ys = collect( range(y_min, y_max, Ny) );

(xm, ym) = util.mesh_grid(xs, ys);

xm = xm'
ym = ym'
# zm = xm .+ ym;

# zm = xm.^2 .+ ym.^2

gmap = gridded.peaks()
zm = gmap.z'
xm = gmap.x
ym = gmap.y

xs = xm
ys = ym

level = 1.5;

# show(stdout, MIME("text/plain"), zm)
# println("")

(accum_array, mark) = compgraph.contour_level(xs, ys, zm, level)

PyPlot.figure()
PyPlot.grid()
PyPlot.contour(xm, ym, zm, [level], linewidths = 5, colors = :red)

PyPlot.plot( accum_array[1,2:end], accum_array[2,2:end], marker = :., markersize = 10 )

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
### INTERNAL FUNCTIONS ###

function padleft(str::String,n)
    out=str
    while length(out)<n
        out=" "*out
    end
    out
end

function padright(str::String,n)
    out=str
    while length(out)<n
        out*=" "
    end
    out
end

function showfloat(x::Float64,dec,n)
    y=x*10^dec
    y=Int(round(y))
    s=string(y)
    if dec<=0
        return padleft(s,n)
    end
    while length(s)<dec+1
        s="0"*s
    end
    r=s[1:end-dec]*"."*s[end-dec+1:end]
    padleft(r,n)
end

function initsim(sim::FastSim,Tdur::Float64,Tout::Float64)
    Tend=sim.t+Tdur
    sim.Tend=Tend
    leaddigits=Int(ceil(log10(Tend)))+1
    decimals=Int(ceil(-log10(Tout)))
    if decimals<1
        decimals=1
    end
    if leaddigits<1
        leaddigits=1
    end
    sim.timedec=decimals
    sim.timedigits=max(leaddigits+decimals+1,4)
    nothing
end

function printlinktype(io::IO,lt::LinkType)
    from=lt.from
    to=lt.to
    if isa(from,Array)
        fromtext=reduce((x,y)->"$x/$y",from)
    else
        if from===0
            fromtext="any"
        else
            fromtext=string(from)
        end
    end
    if isa(to,Array)
        totext=reduce((x,y)->"$x/$y",to)
    else
        if to===0
            totext="any"
        else
            totext=string(to)
        end
    end
    if (lt.dir==2)
        print(io,"Links of the form:  ($fromtext) --- ($totext)" )
    else
        print(io,"Links of the form:  ($fromtext) --> ($totext)" )
    end
end

function linkstate_f!(net::FastNet,kid::Int,cls::Int)
    @inbounds begin
        while net.kstate[kid]<cls
            upstate_link!(net,kid)
        end
        while net.kstate[kid]>cls
            destate_link!(net,kid)
        end
    end
    nothing
end

function destate_node!(net::FastNet,nid) #This function won't check links
    @inbounds begin
        pos=net.npos[nid]
        cls=net.nstate[nid]
        firstpos=net.ncstart[cls]+1
        firstid=net.nid[firstpos]
        net.nid[firstpos]=nid
        net.nid[pos]=firstid
        net.npos[nid]=firstpos
        net.npos[firstid]=pos
        net.ncstart[cls]+=1
        net.nclen[cls]-=1
        net.nclen[cls-1]+=1
        net.nstate[nid]=cls-1
    end
    nothing
end

function upstate_node!(net::FastNet,nid) #This function won't check links
    @inbounds begin
        pos=net.npos[nid]
        cls=net.nstate[nid]
        lastpos=net.ncstart[cls]+net.nclen[cls]
        lastid=net.nid[lastpos]
        net.nid[lastpos]=nid
        net.nid[pos]=lastid
        net.npos[nid]=lastpos
        net.npos[lastid]=pos
        net.ncstart[cls+1]-=1
        net.nclen[cls]-=1
        net.nclen[cls+1]+=1
        net.nstate[nid]=cls+1
    end
    nothing
end

function destate_link!(net::FastNet,kid)
    @inbounds begin
        pos=net.kpos[kid]
        cls=net.kstate[kid]
        firstpos=net.kcstart[cls]+1
        firstid=net.kid[firstpos]
        net.kid[firstpos]=kid
        net.kid[pos]=firstid
        net.kpos[kid]=firstpos
        net.kpos[firstid]=pos
        net.kcstart[cls]+=1
        net.kclen[cls]-=1
        net.kclen[cls-1]+=1
        net.kstate[kid]=cls-1
    end
    nothing
end

function upstate_link!(net::FastNet,kid)
    @inbounds begin
        pos=net.kpos[kid]
        cls=net.kstate[kid]
        lastpos=net.kcstart[cls]+net.kclen[cls]
        lastid=net.kid[lastpos]
        net.kid[lastpos]=kid
        net.kid[pos]=lastid
        net.kpos[kid]=lastpos
        net.kpos[lastid]=pos
        net.kcstart[cls+1]-=1
        net.kclen[cls]-=1
        net.kclen[cls+1]+=1
        net.kstate[kid]=cls+1
    end
    nothing
end

function fastprint(sim::FastSim,x...)
    if isa(sim.printresults,Bool)
        if sim.printresults
            print(x...)
        end
    else
        print(sim.printresults,x...)
    end
end

function _configmodelsetup!(net,degreedist,N,S)
    nullgraph!(net)
    n=0
    try
        n=convert(Int,N)
    catch e
        throw(ArgumentError("Configuration model expects the number of nodes N to be an integer."))
    end
    if n==0
        n=net.N
    end
    if n<1
        throw(ArgumentError("The number of nodes N specified for configuration model is too low"))
    end
    if n<0
        throw(ArgumentError("Configuration model expects the number of nodes, N, to be positive"))
    end
    if n>net.N
        throw(ArgumentError("The requested configuration model network exceeds the max nodecount for the underlying FastNet"))
    end
    checknodestate(net::FastNet,S,"Trying to create configuration model network")
    counts=_countsfromdd(degreedist,n)
    isolated=N-sum(counts)
    totallinks=0
    mx=length(counts)
    for i=1:mx
        totallinks+=i*counts[i]
    end
    if totallinks÷2>net.K
        throw(ArgumentError("The requested configuration model exceeds the max linkcount allowed by the underlying FastNet"))
    end
    makenodes_f!(net,n,S)
    (counts,totallinks)
end

function _randindexdd(rng,list::Union{Array,Tuple})
    L=length(list)
    r=rand(rng)-list[1]
    ret=1
    while r>0.0 && ret<L
        ret+=1
        r-=list[ret]
    end
    if r>0.0
        0
    else
        ret
    end
end

function _randindex(rng,list::Union{Array,Tuple})
    L=length(list)
    t=sum(list)
    r=rand(rng)*t-list[1]
    ret=1
    while r>0.0 && ret<L
        ret+=1
        r-=list[ret]
    end
    ret
end

function _countsfromdd(degreedist,n)
    if !(isa(degreedist,AbstractVector) || isa(degreedist,Tuple))
         throw(ArgumentError("The passed degree distribution needs to be a Vector of Floats"))
    end
    mx=length(degreedist)
    if mx<1
        throw(ArgumentError("The passed degree distribution needs to contain at least one element"))
    end
    dd=Array{Float64,1}(undef,mx)
    for i=1:mx
        try
            dd[i]=convert(Float64,degreedist[i])
        catch e
            throw(ArgumentError("Cannot interpret elment number $i of the degree distribution as a Float variable"))
        end
        if dd[i]<0.0
            throw(ArgumentError("The degree distribution should contain non-negative elements, but element $i seems to be negative"))
        end
        if dd[i]>1.0
            throw(ArgumentError("The degree distribution should be probabilities, but element $i seems to be greater than 1"))
        end
    end
    if sum(dd)>1.0
        throw(ArgumentError("The specified degree distribution sums to more than 1.0"))
    end
    zerop=1.0-sum(dd)
    counts=[Int(floor(dd[i]*n)) for i=1:mx]
    diffs=[dd[i]*n-Float64(counts[i]) for i=1:mx]
    nodessofar=sum(counts)
    nonisolated=Int(round(sum(dd)*n))
    while nodessofar<nonisolated
        (val,idx)=findmax(diffs)
        diffs[idx]=0.0
        nodessofar+=1
        counts[idx]+=1
    end
    stubs=0
    for i=1:mx
        stubs+=i*counts[i]
    end
    if isodd(stubs)               # All this stuff is just for fixing odd numbers of stubs
        peak=(1,2)
        peakdel=-1.0
        diff=1.0
        for i=2:mx
            if (counts[i]==0)
                continue
            end
            for j=1:i-1
                if (counts[j]==0)
                    continue
                end
                del=(counts[i]-dd[i]*n)-(counts[j]-dd[j]*n)
                if abs(del)>peakdel
                    if del<0
                        peakdel=-del
                        diff=-1
                    else
                        peakdel=del
                        diff=1
                    end
                    peakdel=del
                    peak=(i,j)
                end
            end
        end
        counts[peak[1]]-=diff
        counts[peak[2]]+=diff
    end
    stubs=0
    for i=1:mx
        stubs+=i*counts[i]
    end
    if isodd(stubs)
        idx=1
        while idx<mx && count[idx]==0
            idx+=2
        end
        count[idx]-=1
    end
    counts
end

function hhtest(counts)
    c=copy(counts)
    last=length(c)
    adds=zeros(last)
    while last>1
        while c[last]>0
            c[last]-=1
            dist=last
            cur=last
            while dist>0
                if cur<1
                    return false
                end
                sub=dist
                if sub>c[cur]
                    sub=c[cur]
                end
                dist-=sub
                c[cur]-=sub
                if (cur>1)
                    adds[cur-1]=sub
                end
                cur-=1
            end
            for i=1:last
                c[i]+=adds[i]
                adds[i]=0
            end
        end
        last-=1;
    end
    if iseven(c[1])
        true
    else
        false
    end
end


text:
# SimHistory
# maintained by @zsunberg

abstract type SimHistory end
abstract type AbstractMDPHistory{S,A} <: SimHistory end
abstract type AbstractPOMDPHistory{S,A,O,B} <: SimHistory end

# Ordering of a complete step tuple:
# (The general structure is (mdp core, t, mdp info, pomdp core, pomdp info))
const COMPLETE_POMDP_STEP = (:s,:a,:r,:sp,:t,:i,:ai,:b,:o,:bp,:ui)
const COMPLETE_MDP_STEP = COMPLETE_POMDP_STEP[1:7]

"""
An object that contains a MDP simulation history

Returned by simulate when called with a HistoryRecorder. Iterate through the (s, a, r, s') tuples in MDPHistory h like this:

    for (s, a, r, sp) in eachstep(h)
        # do something
    end
"""
struct MDPHistory{S,A} <: AbstractMDPHistory{S,A}
    state_hist::Vector{S}
    action_hist::Vector{A}
    reward_hist::Vector{Float64}
    info_hist::Vector{Any}
    ainfo_hist::Vector{Any}

    discount::Float64

    # if an exception is captured, it will be stored here
    exception::Union{Nothing, Exception}
    backtrace::Union{Nothing, Any}
end

"""
An object that contains a POMDP simulation history

Returned by simulate when called with a HistoryRecorder. Iterate through the (s, b, a, r, s', o) tuples in POMDPHistory h like this:

    for (s, b, a, r, sp, o) in eachstep(h, "s,b,a,r,sp,o")
        # do something
    end
"""
struct POMDPHistory{S,A,O,B} <: AbstractPOMDPHistory{S,A,O,B}
    state_hist::Vector{S}
    action_hist::Vector{A}
    observation_hist::Vector{O}
    belief_hist::Vector{B}
    reward_hist::Vector{Float64}
    info_hist::Vector{Any}
    ainfo_hist::Vector{Any}
    uinfo_hist::Vector{Any}

    discount::Float64

    # if an exception is captured, it will be stored here
    exception::Union{Nothing, Exception}
    backtrace::Union{Nothing, Any}
end

# accessors: use these to access the members - in case the implementation changes
n_steps(h::SimHistory) = length(h.state_hist)-1

state_hist(h::SimHistory) = h.state_hist
action_hist(h::SimHistory) = h.action_hist
observation_hist(h::SimHistory) = h.observation_hist
belief_hist(h::SimHistory) = h.belief_hist
reward_hist(h::SimHistory) = h.reward_hist
info_hist(h::SimHistory) = h.info_hist
ainfo_hist(h::SimHistory) = h.ainfo_hist
uinfo_hist(h::SimHistory) = h.uinfo_hist

exception(h::SimHistory) = h.exception
Base.backtrace(h::SimHistory) = h.backtrace
POMDPs.discount(h::SimHistory) = h.discount

undiscounted_reward(h::SimHistory) = sum(reward_hist(h))
function discounted_reward(h::SimHistory)
    disc = 1.0
    r_total = 0.0
    for i in 1:length(reward_hist(h))
        r_total += disc*reward_hist(h)[i]
        disc *= discount(h)
    end
    return r_total
end

# iteration
Base.length(h::SimHistory) = n_steps(h)
function Base.iterate(h::SimHistory, i::Int=1)
    if i > n_steps(h)
        return nothing
    else
        return (step_tuple(h, i), i+1)
    end
end

function step_tuple(h::MDPHistory, i::Int)
    return (s=state_hist(h)[i],
            a=action_hist(h)[i],
            r=reward_hist(h)[i],
            sp=state_hist(h)[i+1]
           )
end
function step_tuple(h::POMDPHistory, i::Int)
    return (s=state_hist(h)[i],
            b=belief_hist(h)[i],
            a=action_hist(h)[i],
            r=reward_hist(h)[i],
            sp=state_hist(h)[i+1],
            o=observation_hist(h)[i]
           )
end

const Inds = Union{AbstractRange,Colon,Real}

Base.view(h::AbstractMDPHistory, inds::Inds) = SubMDPHistory(h, inds)
Base.view(h::AbstractPOMDPHistory, inds::Inds) = SubPOMDPHistory(h, inds)

struct SubMDPHistory{S,A,H<:AbstractMDPHistory,I<:Inds} <: AbstractMDPHistory{S,A}
    parent::H
    inds::I
end
SubMDPHistory(h::AbstractMDPHistory{S,A}, inds::I) where {S,A,I<:Inds} = SubMDPHistory{S,A,typeof(h),I}(h, inds)

struct SubPOMDPHistory{S,A,O,B,H<:AbstractPOMDPHistory,I<:Inds} <: AbstractPOMDPHistory{S,A,O,B}
    parent::H
    inds::I
end
SubPOMDPHistory(h::AbstractPOMDPHistory{S,A,O,B}, inds::I) where {S,A,O,B,I<:Inds} = SubPOMDPHistory{S,A,O,B,typeof(h),I}(h, inds)

const SubHistory = Union{SubMDPHistory, SubPOMDPHistory}

n_steps(h::SubHistory) = length(h.inds)

state_hist(h::SubHistory) = state_hist(h.parent)[minimum(h.inds):maximum(h.inds)+1]
action_hist(h::SubHistory) = action_hist(h.parent)[h.inds]
observation_hist(h::SubHistory) = observation_hist(h.parent)[h.inds]
belief_hist(h::SubHistory) = belief_hist(h.parent)[h.inds]
reward_hist(h::SubHistory) = reward_hist(h.parent)[h.inds]
info_hist(h::SubHistory) = info_hist(h.parent)[h.inds]
ainfo_hist(h::SubHistory) = ainfo_hist(h.parent)[h.inds]
uinfo_hist(h::SubHistory) = uinfo_hist(h.parent)[h.inds]

step_tuple(h::SubHistory, i::Int) = step_tuple(h.parent, h.inds[i])

exception(h::SubHistory) = exception(h.parent)
Base.backtrace(h::SubHistory) = backtrace(h.parent)
POMDPs.discount(h::SubHistory) = discount(h.parent)

# iterators
struct HistoryIterator{H<:SimHistory, SPEC}
    history::H
end

# Note this particular function is not type-stable
function HistoryIterator(history::SimHistory, spec::String)
    # XXX should throw warnings for unrecognized specification characters
    syms = [Symbol(m.match) for m in eachmatch(r"(sp|bp|ai|ui|s|a|r|b|o|i|t)", spec)]
    if length(syms) == 1
        return HistoryIterator{typeof(history), first(syms)}(history)
    else
        return HistoryIterator{typeof(history), tuple(syms...)}(history)
    end
end

function HistoryIterator(history::SimHistory, spec::Tuple)
    @assert all(isa(s, Symbol) for s in spec)
    return HistoryIterator{typeof(history), spec}(history)
end
HistoryIterator(h::SimHistory, spec::Symbol) = HistoryIterator{typeof(h), spec}(h)

"""
    for t in eachstep(hist, [spec])
        ...
    end

Iterate through the steps in `SimHistory` `hist`. `spec` is a tuple of symbols or string that controls what is returned for each step.

For example,
```julia
for (s, a, r, sp) in eachstep(h, "(s, a, r, sp)")
    println("reward \$r received when state \$sp was reached after action \$a was taken in state \$s")
end
```
returns the start state, action, reward and destination state for each step of the simulation.

Alternatively, instead of expanding the steps implicitly, the elements of the step can be accessed as fields (since each step is a `NamedTuple`):
```julia
for step in eachstep(h, "(s, a, r, sp)")
    println("reward \$(step.r) received when state \$(step.sp) was reached after action \$(step.a) was taken in state \$(step.s)")
end
```

The possible valid elements in the iteration specification are
- `s` - the initial state in a step
- `b` - the initial belief in the step (for POMDPs only)
- `a` - the action taken in the step
- `r` - the reward received for the step
- `sp` - the final state at the end of the step (s')
- `o` - the observation received during the step (note that this is usually based on `sp` instead of `s`)
- `bp` - the belief after being updated based on `o` (for POMDPs only)
- `i` - info from the state transition (from `generate_sri` for MDPs or `generate_sori` for POMDPs)
- `ai` - info from the policy decision (from `action_info`)
- `ui` - info from the belief update (from `update_info`)
- `t` - the timestep index
"""
eachstep(hist::SimHistory, spec) = HistoryIterator(hist, spec)

eachstep(mh::AbstractMDPHistory) = eachstep(mh, COMPLETE_MDP_STEP)
eachstep(mh::AbstractPOMDPHistory) = eachstep(mh, COMPLETE_POMDP_STEP)

function sym_to_call(sym::Symbol)
    if sym == :s
        return :(state_hist(it.history)[i])
    elseif sym == :a
        return :(action_hist(it.history)[i])
    elseif sym == :r
        return :(reward_hist(it.history)[i])
    elseif sym == :sp
        return :(state_hist(it.history)[i+1])
    elseif sym == :b
        return :(belief_hist(it.history)[i])
    elseif sym == :o
        return :(observation_hist(it.history)[i])
    elseif sym == :bp
        return :(belief_hist(it.history)[i+1])
    elseif sym == :i
        return :(info_hist(it.history)[i])
    elseif sym == :ai
        return :(ainfo_hist(it.history)[i])
    elseif sym == :ui
        return :(uinfo_hist(it.history)[i])
    elseif sym == :t
        return :(i)
    end
end

@generated function step_tuple(it::HistoryIterator, i::Int)
    spec = it.parameters[2]
    if isa(spec, Tuple)
        calls = []
        for sym in spec
            push!(calls, :($sym = $(sym_to_call(sym))))
        end

        return quote
            return ($(calls...),)
        end
    else
        @assert isa(spec, Symbol)
        return quote
            return $(sym_to_call(spec))
        end
    end
end

Base.length(it::HistoryIterator) = n_steps(it.history)
Base.getindex(it::HistoryIterator, i) = step_tuple(it, i)
function Base.iterate(it::HistoryIterator, i::Int = 1)
    if i > length(it)
        return nothing
    else
        return (step_tuple(it, i), i+1)
    end
end


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# int tbx_name2id(tbx_t *tbx, const char *ss);
@defun tbx_name2id(tbx::Ptr{tbx_t}, ss::Cstring)::Cint

# BGZF *hts_get_bgzfp(htsFile *fp);
@defun hts_get_bgzfp(fp::Ptr{htsFile})::Ptr{BGZF}

# int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_pos_t *end);
@defun tbx_readrec(fp::Ptr{BGZF}, tbxv::Ptr{Cvoid}, sv::Ptr{Cvoid}, tid::Ptr{Cint}, beg::Ptr{hts_pos_t}, end_::Ptr{hts_pos_t})::Cint

# tbx_t *tbx_index(BGZF *fp, int min_shift, const tbx_conf_t *conf);
@defun tbx_index(fp::Ptr{BGZF}, min_shift::Cint, conf::Ptr{tbx_conf_t})::Ptr{tbx_t}

# int tbx_index_build(const char *fn, int min_shift, const tbx_conf_t *conf);
@defun tbx_index_build(fn::Cstring, min_shift::Cint, conf::Ptr{tbx_conf_t})::Cint

# int tbx_index_build2(const char *fn, const char *fnidx, int min_shift, const tbx_conf_t *conf);
@defun tbx_index_build2(fn::Cstring, fnidx::Cstring, min_shift::Cint, conf::Ptr{tbx_conf_t})::Cint

# int tbx_index_build3(const char *fn, const char *fnidx, int min_shift, int n_threads, const tbx_conf_t *conf);
@defun tbx_index_build3(fn::Cstring, fnidx::Cstring, min_shift::Cint, n_threads::Cint, conf::Ptr{tbx_conf_t})::Cint

# tbx_t *tbx_index_load(const char *fn);
@defun tbx_index_load(fn::Cstring)::Ptr{tbx_t}

# tbx_t *tbx_index_load2(const char *fn, const char *fnidx);
@defun tbx_index_load2(fn::Cstring, fnidx::Cstring)::Ptr{tbx_t}

# tbx_t *tbx_index_load3(const char *fn, const char *fnidx, int flags);
@defun tbx_index_load3(fn::Cstring, fnidx::Cstring, flags::Cint)::Ptr{tbx_t}

# const char **tbx_seqnames(tbx_t *tbx, int *n);
@defun tbx_seqnames(tbx::Ptr{tbx_t}, n::Ptr{Cint})::Ptr{Cstring}

# void tbx_destroy(tbx_t *tbx);
@defun tbx_destroy(tbx::Ptr{tbx_t})::Cvoid


text:
"""
    split_model(model::EoSModel)

Takes in a model for a multi-component system and returns a vector of model for each pure system.

## Example:
```julia-repl
julia> gerg2 = GERG2008(["propane","pentane"])
GERG008 model with 2 components:
"propane"
"pentane"

julia> split_model(gerg2)
2-element Vector{GERG2008}:
 GERG2008("propane")
 GERG2008("pentane")
```
"""
function split_model end

"""
    is_splittable(model)::Bool

Trait to determine if a `EoSModel` should be splitted by itself or can be simply filled into a vector.
This is useful in the case of models without any parameters, as those models are impossible by definition to split, because they don't have any underlying data.

The Default is `is_splittable(model) = true`.
"""
is_splittable(model) = true
is_splittable(null::Union{Nothing,Missing}) = false
is_splittable(::Number) = false
is_splittable(::String) = false

function split_model(param::SingleParam{T},
    splitter =split_model(1:length(param.components))) where T
    function generator(I)
    return SingleParam(
    param.name,
    param.components[I],
    deepcopy(param.values[I]),#for SingleParam{Vector{X}}
    param.ismissingvalues[I],
    param.sourcecsvs,
    param.sources
    )
    end
    return [generator(i) for i ∈ splitter]
end

function split_model(param::ClapeyronParam,groups::GroupParam)
    return split_model(param,groups.i_groups)
end
#this conversion is lossy, as interaction between two or more components are lost.

function split_model(param::PairParam{T},
    splitter = split_model(1:length(param.components))) where T
    function generator(I)
        value = param.values[I,I]
        diagvalue = view(value,diagind(value))
        ismissingvalues = param.ismissingvalues[I,I]
        components = param.components[I]
        return PairParam{T}(
                param.name,
                components,
                value,
                diagvalue,
                ismissingvalues,
                param.sourcecsvs,
                param.sources
                )
    end
    return [generator(I) for I ∈ splitter]
end

function split_model(param::AbstractVector,splitter = ([i] for i ∈ 1:length(param)))
    return [param[i] for i ∈ splitter]
end

function split_model(param::UnitRange{Int},splitter = ([i] for i ∈ 1:length(param)))
    return [1:length(i) for i ∈ splitter]
end

function _split_model(assoc::CompressedAssocMatrix{T},I) where T
    len = length(assoc.values)
    iszero(len) && return CompressedAssocMatrix{T}()
    old_idx = assoc.outer_indices
    idx_bool = findall(x -> (first(x) ∈ I)&(last(x) ∈ I),old_idx)
    iszero(length(idx_bool)) && return CompressedAssocMatrix{T}()
    values = assoc.values[idx_bool]
    outer_indices = assoc.outer_indices[idx_bool]
    inner_indices = assoc.inner_indices[idx_bool]
    out_val = length(I)
    outer_size = (out_val,out_val)
    inner_size = assoc.inner_size
    len2 = length(outer_indices)
    for i ∈ 1:len2
        i1,j1 = outer_indices[i]
        i2,j2 = findfirst(==(i1),I),findfirst(==(j1),I)
        outer_indices[i] = (i2,j2)
    end
    return CompressedAssocMatrix(values,outer_indices,inner_indices,outer_size,inner_size)
end

#this conversion is lossy, as interaction between two or more components are lost.
#also, this conversion stores the site values for other components. (those are not used)
function split_model(param::AssocParam{T},
    splitter = split_model(1:length(param.components))) where T
    function generator(I)
        _value  = _split_model(param.values,I)

        return AssocParam{T}(
                param.name,
                param.components[I],
                _value,
                param.sites[I],
                param.sourcecsvs,
                param.sources
                )
        end
    return [generator(I) for I ∈ splitter]
end

#this param has a defined split form
function split_model(groups::GroupParam)
    len = length(groups.components)
    function generator(i)
        return GroupParam(
        [groups.components[i]],
        [groups.groups[i]],
        [groups.n_groups[i]],
        [collect(1:length(groups.n_groups[i]))],
        groups.flattenedgroups[groups.i_groups[i]],
        [groups.n_groups[i]],
        1:length(groups.n_groups[i]),
        groups.sourcecsvs)
    end
    [generator(i) for i ∈ 1:len]
end

function split_model(param::SiteParam,
    splitter = split_model(param.components))
    function generator(I)
        return SiteParam(
            param.components[I],
            param.sites[I],
            param.n_sites[I],
            param.i_sites[I],
            param.flattenedsites,
            param.n_flattenedsites[I],
            param.i_flattenedsites[I],
        param.sourcecsvs)
    end
    return [generator(i) for i ∈ splitter]
end

function split_model(Base.@nospecialize(params::EoSParam),splitter)
    T = typeof(params)
    split_paramsvals = [split_model(getfield(params,i),splitter) for i  ∈ fieldnames(T)]
    return T.(split_paramsvals...)
end

export SingleParam, SiteParam, PairParam, AssocParam, GroupParam
#

split_model(model::EoSModel,subset=nothing) = auto_split_model(model,subset)

function auto_split_model(Base.@nospecialize(model::EoSModel),subset=nothing)
    try
        allfields = Dict{Symbol,Any}()
        if has_groups(typeof(model))
            raw_splitter = model.groups.i_groups
            subset !== nothing && throw("using subsets is not supported with Group Contribution models")
        else
            raw_splitter = split_model(Vector(1:length(model.components)))
        end
        if subset === nothing
            splitter = raw_splitter
        elseif eltype(subset) <: Integer
            splitter = raw_splitter[subset]
        elseif eltype(subset) <: AbstractVector
            splitter = subset
        else
            throw("invalid type of subset.")
        end

        len = length(splitter)
        M = typeof(model)
        allfieldnames = fieldnames(M)
        if hasfield(typeof(model),:groups) #TODO implement a splitter that accepts a subset
            allfields[:groups] = split_model(model.groups)
            allfields[:components] = getfield.(allfields[:groups],:components)
            for modelkey in allfieldnames
                if !haskey(allfields,modelkey)
                    modelx = getproperty(model,modelkey)
                    if modelx isa Vector{<:EoSModel}
                        allfields[modelkey] = split_model(modelx) #this is an api problem
                    end
                end
            end
        end

        #add here any special keys, that behave as non_splittable values
        for modelkey in [:references]
            if modelkey in allfieldnames
                if !haskey(allfields,modelkey)
                    allfields[modelkey] = fill(getproperty(model,modelkey),len)
                end
            end
        end

        for modelkey ∈ allfieldnames
            if !haskey(allfields,modelkey)
                modelx = getproperty(model,modelkey)
                if is_splittable(modelx)
                    if modelx isa EoSModel
                        allfields[modelkey]= split_model(modelx,subset)
                    else
                        allfields[modelkey]= split_model(modelx,splitter)
                    end
                else
                    allfields[modelkey] = fill(modelx,len)
                end
            end
        end

        return [M((allfields[k][i] for k ∈ fieldnames(M))...) for i ∈ 1:len]
    catch e
        M = typeof(model)
        @error "$M cannot be splitted"
        rethrow(e)
    end
end

##fallback,around 50 times slower if there is any need to read csvs

function simple_split_model(Base.@nospecialize(model::EoSModel),subset = nothing)
    MODEL = typeof(model)
    pure = Vector{MODEL}(undef,0)
    if subset === nothing
        comps = model.components
    else
        comps = model.components[subset]
    end
    for comp ∈ comps
        push!(pure,MODEL([comp]))
    end
    return pure
end

export split_model


text:
"""
# TerminalButtons.jl

Terminal-mode Push Buttons for Linux Touch Screens.

Create simple touch screen interfaces without X11.

Uses [LinuxTouchEvents.jl](https://github.com/notinaboat/LinuxTouchEvents.jl)
to get touch screen input (e.g from a
[Raspberry Pi Touch Display](https://www.raspberrypi.org/documentation/hardware/display/README.md)).

Uses [Terming.jl](https://github.com/foldfelis/Terming.jl) and
[Crayons.jl](https://github.com/KristofferC/Crayons.jl) to draw buttons.

![](buttons.png)
"""
module TerminalButtons

using LinuxTouchEvents
using Crayons
using UnixIO
using UnixIO.Debug
using Terming

struct Rect
    x::Int
    y::Int
    width::Int
    height::Int
end

function Base.in(p, r::Rect)
    x, y = p
    x >= r.x && x <= r.x + r.width &&
    y >= r.y && y <= r.y + r.height
end

mutable struct Button
    id::Any
    text::String
    style::Crayon
    active_style::Crayon
    shape::Union{Nothing,Rect}
    state::Union{Nothing,Symbol}
    function Button(id, text,
                    style = crayon"fg:white bg:blue",
                    active_style = crayon"fg:white bg:yellow") where T
        new(id, text, style, active_style, nothing, nothing)
    end
end

"""
    divide(x, n) -> [a, b, c ...]

Divide `x` into a `n`-element vector such that `sum(divide(x, n)) == x`
"""
function divide(x, n)
    d = x ÷ n
    err = x - (n * d)
    result = fill(d, n)
    for i in 1:err
        result[1 + (i-1) % n] += 1
    end
    result
end

function render_horizontal(io, buttons::Vector{Button}, shape::Rect)

    # Fit buttons evenly into available width...
    widths = divide(shape.width + 1, length(buttons))

    # Update shape of each button for touch detection...
    x = shape.x
    for (w, b) in zip(widths, buttons)
        b.shape = Rect(x, shape.y, w-1, shape.height)
        x += w
    end

    # Draw buttons...
    y = shape.y
    for i in 1:shape.height
        cmove(io, y, shape.x)
        for (j, b) in enumerate(buttons)
            j == 1 || print(io, " ")
            render_line(io, b, b.style, i)
        end
        y += 1
    end
end

function render_line(io, b::Button, style, row)
    print(io, style)
    if row == (b.shape.height + 1) ÷ 2
        # Draw Text on middle row,
        pad = b.shape.width - length(b.text)
        print(io, " "^(pad÷2), b.text, " "^(pad-(pad÷2)))
    else
        #...otherwise pad with spaces.
        print(io, " "^b.shape.width)
    end
    print(io, inv(style))
end

function render_vertical(io, buttons::Vector{Button}, shape::Rect)

    # Fit buttons evenly into available height...
    heights = divide(shape.height + 1, length(buttons))

    # Update shape of each button for touch detection...
    y = shape.y
    for (h, b) in zip(heights, buttons)
        b.shape = Rect(shape.x, y, shape.width, h-1)
        y += h
    end

    # Draw buttons...
    for (i, b) in enumerate(buttons)
        if i > 1
            cmove_line_down(io)
            clear_line(io)
        end
        draw_button(io, b)
    end
end

function draw_button(io, b::Button, style=b.style)
    y = b.shape.y
    for i in 1:b.shape.height
        cmove(io, y, b.shape.x)
        render_line(io, b, style, i)
        y += 1
    end
end

function select_button(io, buttons::Vector{Button}, touch)
    for b in buttons
        if touch in b.shape
            draw_button(io, b, b.active_style)
            return b
        end
    end
    return nothing
end

"""
    choose_button(io, [id => button, ...]; vertical=false) -> selected_button

Draw a selection of buttons and wait for one of them to be pressed.

e.g.

    x = TerminalButtons.choose_button([1=>"Abort", 2=>"Retry", 3=>"Fail"])
    if x == 1
        ...
    elseif x == 2
        ...
    elseif x == 3
        ...
    end
"""
@db function choose_button(io, buttons;
                           vertical=false, rect=nothing, size=nothing)

    if size == nothing
        size = UnixIO.tiocgwinsz(io)
        screen_h, screen_w = size.ws_row, size.ws_col
    else
        screen_h, screen_w = size
    end
                                                        @db 3 screen_h screen_w
    if rect == nothing
        if vertical
            rect = Rect(1, 1, screen_w, screen_h)
        else
            # Draw buttons across the bottom of the screen by default.
            button_h = 3
            rect = Rect(1, screen_h  + 1 - button_h, screen_w, button_h)
        end
    end

    # Draw buttons.
    buttons = [Button(id, text) for (id, text) in buttons]
    if vertical
        render_vertical(io, buttons, rect)
    else
        render_horizontal(io, buttons, rect)
    end

    # Wait for touch event.
    c = TouchEventChannel()                                            ;@db 3 c
    while true
        x, y = take!(c)
        x = round(Int, x * screen_w)
        y = round(Int, y * screen_h)

        # Select button at touch location.
        result = select_button(io, buttons, (x, y))
        if result != nothing
            close(c)
            @db return result.id
        end
    end
end

choose_button(io, v::Vector{String}; kw...) =
    choose_button(io, [x => x for x in v]; kw...)

choose_button(io, v::Vector{Pair{String,Function}}; kw...) =
    choose_button(io, [f => s for (s, f) in v]; kw...)()

# Documentation.

readme() = join([
    Docs.doc(@__MODULE__),
    "## Interface\n",
    Docs.@doc choose_button
   ], "\n\n")

end # module


text:
# module Bukdu.Plug

import ..ApplicationRouter, ..Routing, ..RouterRoute

"""
plug `Router` to give the routes into the Endpoint.

```julia
Endpoint() do
    plug(Router)
end
```
"""
function plug{AR<:ApplicationRouter}(::Type{AR})
    if haskey(Routing.routing_map, AR)
        append!(RouterRoute.routes, Routing.routing_map[AR])
    end
end


text:
function assembly_Theta(X_domain, X_boundary, set_sigma)
    N_domain = size(X_domain)[1]
    N_boundary = size(X_boundary)[1]
    X_all = [X_domain;X_boundary]
    Theta = zeros((2*N_domain + N_boundary, 2*N_domain + N_boundary))

    for iter_i = 1:N_domain
        for iter_j = 1:N_domain
            Theta[iter_i,iter_j] = Delta_x_y_kappa(X_domain[iter_i,1], X_domain[iter_i,2], X_domain[iter_j,1], X_domain[iter_j,2], set_sigma)
        end
        # for loop is good!
        for iter_j = N_domain+1:2*N_domain+N_boundary
            Theta[iter_i,iter_j] = Delta_x_kappa(X_domain[iter_i,1], X_domain[iter_i,2], X_all[iter_j-N_domain,1], X_all[iter_j-N_domain,2], set_sigma)
            Theta[iter_j,iter_i] = Theta[iter_i,iter_j]
        end
    end
    for iter_i = N_domain+1:2*N_domain+N_boundary
        for iter_j = N_domain+1:2*N_domain+N_boundary
            Theta[iter_i,iter_j] = kappa(X_all[iter_i-N_domain,1], X_all[iter_i-N_domain,2], X_all[iter_j-N_domain,1], X_all[iter_j-N_domain,2], set_sigma)
        end
    end
    return Theta
end


text:
using Test, YaoBase, YaoBlocks, BitBasis, YaoArrayRegister, Random
using StatsBase: mean

@testset "measure ghz" begin
    # GHZ state
    st = normalize!(arrayreg(bit"0000") + arrayreg(bit"1111"))
    Random.seed!(1234)

    # measure it at 1, 2
    # should collapse to 0000 or 1111 since entangled
    g = Measure(4; locs = (1, 2))
    @test occupied_locs(g) == (1, 2)
    st |> g
    @test (copy(st) |> g) isa ArrayReg

    @test g.results[1] == 0 ? st.state[end] == 0 : st.state[1] == 0
    g = Measure(4; locs = (1, 2), resetto = 2)
    @test g.postprocess isa ResetTo{BitStr64{2}}

    m = Measure(4)
    @test occupied_locs(m) == (1, 2, 3, 4)
    @test m.postprocess isa NoPostProcess
end

@testset "resetto" begin
    Random.seed!(1234)

    st = rand_state(5; nbatch = 3)
    g = Measure(5; locs = (1, 2), resetto = bit"00011")
    @test num_measured(g) == 2
    st |> g
    for k = 1:32
        if !(st.state[k] ≈ 0.0)
            @test all(BitStr64{5}(k - 1)[1:2] .== 1)
        end
    end
    @test Measure(5; locs = (1, 2), resetto = 0b0011).postprocess isa ResetTo{BitStr64{2}}
end

@testset "error handling" begin
    @test_throws ErrorException Measure(
        5;
        locs = (1, 2),
        resetto = bit"00011",
        remove = true,
    )
    @test_throws ErrorException mat(Measure(5; locs = (1, 2), resetto = bit"00011"))
end

@testset "op-measures" begin
    Random.seed!(1234)

    reg = rand_state(6, nbatch = 10)
    op = repeat(3, X)

    # measure!
    reg2 = reg |> copy
    res = measure!(op, reg2, 2:4)
    res2 = measure!(op, reg2, 2:4)
    @test size(res) == (10,)
    @test res2 == res

    # measure_resetto! and measure_remove! for operators are no-longer supported due to its ill property.
    reg = repeat(ArrayReg(ComplexF64[1, -1] / sqrt(2.0)), 10)
    @test measure!(X, reg) |> mean ≈ -1
    reg = repeat(ArrayReg(ComplexF64[1.0, 0]), 1000)
    @test abs(measure!(X, reg) |> mean) < 0.1
end

@testset "op-measures" begin
    reg = rand_state(8; nbatch = 32)
    op = repeat(5, X, 1:5)

    # measure!
    reg2 = reg |> copy
    res = measure!(op, reg2, 2:6)
    res2 = measure!(op, reg2, 2:6)
    @test size(res) == (32,)
    @test res2 == res

    # measure
    reg2 = reg |> copy
    res = measure(op, reg2, 2:6; nshots = 100)
    @test size(res) == (100, 32)
    @test reg ≈ reg2

    reg = repeat(ArrayReg([1, -1 + 0im] / sqrt(2.0)), 10)
    @test measure!(X, reg) |> mean ≈ -1
    reg = repeat(ArrayReg([1.0, 0 + 0im]), 1000)
    @test abs(measure!(X, reg) |> mean) < 0.1

    m = Measure(5)
    @test chmeasureoperator(m, X) == Measure(5, operator = X)
end

@testset "measure an operator correctly" begin
    op = kron(Z, Z)
    reg = ArrayReg(ComplexF64[1/sqrt(2), 0, 0, 1/sqrt(2)])
    res = measure!(op, reg)
    @test res == 1
    @test reg ≈ ArrayReg(ComplexF64[1/sqrt(2), 0, 0, 1/sqrt(2)])
    reg = uniform_state(2)
    res = measure!(op, reg)
    @test count(!iszero, reg.state) == 2

    # batched
    reg =
        BatchedArrayReg(reshape(ComplexF64[1/sqrt(2), 0, 0, 1/sqrt(2), 0.5, 0.5, 0.5, 0.5], 4, 2))
    res = measure!(op, reg)
    @test length(reg) == 2 && res[1] == 1
    @test reg.state[:, 1] ≈ ComplexF64[1/sqrt(2), 0, 0, 1/sqrt(2)]
    @test count(!iszero, reg.state) == 4
    @test isnormalized(reg)

    # with virtual dimension
    reg =
        ArrayReg(
            reshape(ComplexF64[1/sqrt(2), 0, 0, 1/sqrt(2), 0.5, 0.5, 0.5, 0.5], 4, 2),
        ) / sqrt(2)
    res = measure!(op, reg)
    @test reg isa ArrayReg
    c = count(!iszero, reg.state)
    @test (c == 4 && res ≈ 1) || (c == 2 && res == -1)
    @test isnormalized(reg)

    # measure zero space
    reg = ArrayReg(ComplexF64[1.0])
    res = measure!(matblock(fill(3.0 + 0im, 1, 1)), reg)
    @test res ≈ 3.0
    @test reg == zero_state(0)
end


text:
include("processors/assets.jl")
include("processors/callback.jl")
include("processors/dependecies.jl")
include("processors/index.jl")
include("processors/layout.jl")
include("processors/reload_hash.jl")
include("processors/resource.jl")

function start_reload_poll(state::HandlerState)
    folders = Set{String}()
    push!(folders, get_assets_path(state.app))
    push!(folders, state.registry.dash_renderer.path)
    for pkg in values(state.registry.components)
        push!(folders, pkg.path)
    end
    state.reload.task = @async poll_folders(folders; interval = get_devsetting(state.app, :hot_reload_watch_interval)) do file, ts, deleted
        state.reload.hard = true
        state.reload.hash = generate_hash()
        assets_path = get_assets_path(state.app)
        if startswith(file, assets_path)
            state.cache.need_recache = true
            rel_path = lstrip(
                    replace(relpath(file, assets_path), '\\'=>'/'),
                    '/'
                )
            push!(state.reload.changed_assets,
                ChangedAsset(
                    asset_path(state.app, rel_path),
                    trunc(Int, ts),
                    endswith(file, "css")
                    )
            )
        end

    end
end

validate_layout(layout::Component) = validate(layout)

validate_layout(layout::Function) = validate_layout(layout())

validate_layout(layout) = error("The layout must be a component, tree of components, or a function which returns a component.")

function stack_from_callback(full_stack; prune_errors)
    !prune_errors && return full_stack

    result = Base.StackTraces.StackFrame[]
    it = iterate(full_stack)
    while !isnothing(it) && it[1].func != :process_callback_call
        push!(result, it[1])
        it = iterate(full_stack, it[2])
    end
    return result
end

function exception_handling(ex; prune_errors)
    st = stack_from_callback(stacktrace(catch_backtrace()), prune_errors = prune_errors)
    @error "error handling request" exception = (ex, st)
    return HTTP.Response(500)
end

function debug_exception_handling(ex; prune_errors)
    response = HTTP.Response(500, ["Content-Type" => "text/html"])
    io = IOBuffer()
    write(io,
         "<!DOCTYPE HTML>",
         "",
         "",
         "",
         "",
         "
",
         "
",
         "
")
    showerror(io, ex)
    write(io,
        "
",
        "
",
        "
",
        "
",
        "<textarea cols=\"50\" rows=\"10\" name=\"code\" readonly>"
    )
    st = stack_from_callback(stacktrace(catch_backtrace()), prune_errors = prune_errors)
    showerror(io, ex)
    write(io, "\n")
    Base.show_backtrace(io, st)

    write(io, "</textarea>")

    write(io, "")

    response.body = take!(io)

    @error "error handling request"  exception = (ex, st)

    return response
end

#For test purposes, with the ability to pass a custom registry
function make_handler(app::DashApp, registry::ResourcesRegistry; check_layout = false)

    state = HandlerState(app, registry)
    prefix = get_setting(app, :routes_pathname_prefix)
    assets_url_path = get_setting(app, :assets_url_path)
    check_layout && validate_layout(get_layout(app))

    router = Router()
    add_route!(process_layout, router, "$(prefix)_dash-layout")
    add_route!(process_dependencies, router, "$(prefix)_dash-dependencies")
    add_route!(process_reload_hash, router, "$(prefix)_reload-hash")
    add_route!(process_resource, router, "$(prefix)_dash-component-suites/<namespace>/<path>")
    add_route!(process_assets, router, "$(prefix)$(assets_url_path)/<file_path>")
    add_route!(process_callback, router, "POST", "$(prefix)_dash-update-component")
    add_route!(process_index, router, "$prefix/*")
    add_route!(process_index, router, "$prefix")

    handler = state_handler(router, state)
    if get_devsetting(app, :ui)
        handler = exception_handling_handler(handler) do ex
            debug_exception_handling(ex, prune_errors = get_devsetting(app, :prune_errors))
        end
    else
        handler = exception_handling_handler(handler) do ex
            exception_handling(ex, prune_errors = get_devsetting(app, :prune_errors))
        end
    end
    get_setting(app, :compress) && (handler = compress_handler(handler))

    compile_request = HTTP.Request("GET", prefix)
    HTTP.setheader(compile_request, "Accept-Encoding" => "gzip")
    HTTP.handle(handler, compile_request) #For handler precompilation

    get_devsetting(app, :hot_reload) && start_reload_poll(state)

    return handler
end

function make_handler(app::DashApp)
    return make_handler(app, main_registry(), check_layout = true)
end


text:
export monovec, monovectype, emptymonovec, sortmonovec, mergemonovec

monomials(v::AbstractVariable, args...) = monomials((v,), args...)

"""
    emptymonovec(p::AbstractPolynomialLike)

Returns an empty collection of the type of `monomials(p)`.
"""
emptymonovec(p) = monomialtype(p)[]

"""
    monovec(X::AbstractVector{MT}) where {MT<:AbstractMonomialLike}

Returns the vector of monomials `X` in decreasing order and without any duplicates.

### Examples

Calling `monovec` on ``[xy, x, xy, x^2y, x]`` should return ``[x^2y, xy, x]``.
"""
function monovec(X::AbstractVector{MT}) where {MT<:AbstractMonomial}
    Y = sort(X, rev=true)
    dups = findall(i -> Y[i] == Y[i-1], 2:length(Y))
    deleteat!(Y, dups)
    Y
end
monovec(X::AbstractVector{TT}) where {TT<:AbstractTermLike} = monovec(AbstractVector{monomialtype(TT)}(X))

"""
    monovec(a, X::AbstractVector{MT}) where {MT<:AbstractMonomialLike}

Returns `b, Y` where `Y` is the vector of monomials of `X` in decreasing order
and without any duplicates and `b` is the vector of corresponding coefficients
in `a`, where coefficients of duplicate entries are summed together.

### Examples

Calling `monovec` on ``[2, 1, 4, 3, -1], [xy, x, xy, x^2y, x]`` should return
``[3, 6, 0], [x^2y, xy, x]``.
"""
function monovec(a, x)
    if length(a) != length(x)
        throw(ArgumentError("There should be as many coefficient than monomials"))
    end
    σ, X = sortmonovec(x)
    b = a[σ]
    if length(x) > length(X)
        rev = Dict(X[j] => j for j in eachindex(σ))
        for i in eachindex(x)
            j = rev[x[i]]
            if i != σ[j]
                b[j] += a[i]
            end
        end
    end
    return b, X
end

"""
    monovectype(X::AbstractVector{MT}) where {MT<:AbstractMonomialLike}

Returns the return type of `monovec`.
"""
monovectype(X::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where {PT<:APL} = monovectype(PT)
monovectype(::Union{PT, Type{PT}}) where {PT <: APL} = Vector{monomialtype(PT)}

# If there are duplicates in X, the coefficients should be summed for a polynomial and they should be equal for a measure.
"""
    sortmonovec(X::AbstractVector{MT}) where {MT<:AbstractMonomialLike}

Returns `σ`, the orders in which one must take the monomials in `X` to make them sorted and without any duplicate and the sorted vector of monomials, i.e. it returns `(σ, X[σ])`.

### Examples

Calling `sortmonovec` on ``[xy, x, xy, x^2y, x]`` should return ``([4, 1, 2], [x^2y, xy, x])``.
"""
function sortmonovec(X::AbstractVector{MT}) where {MT<:AbstractMonomial}
    σ = sortperm(X, rev=true)
    dups = findall(i -> X[σ[i]] == X[σ[i-1]], 2:length(σ))
    deleteat!(σ, dups)
    σ, X[σ]
end
sortmonovec(X::AbstractVector{TT}) where {TT<:AbstractTermLike} = sortmonovec(AbstractVector{monomialtype(TT)}(X))
sortmonovec(X::Tuple) = sortmonovec(vec(X))

"""
    mergemonovec{MT<:AbstractMonomialLike, MVT<:AbstractVector{MT}}(X::AbstractVector{MVT}}

Returns the vector of monomials in the entries of `X` in decreasing order and without any duplicates, i.e. `monovec(vcat(X...))`

### Examples

Calling `mergemonovec` on ``[[xy, x, xy], [x^2y, x]]`` should return ``[x^2y, xy, x]``.
"""
mergemonovec(X) = monovec(vcat(X...))


text:
export TinkerforgeSilentStepperParams, TinkerforgeSilentStepper

Base.@kwdef struct TinkerforgeSilentStepperParams <: DeviceParams
  host::IPAddr = ip"127.0.0.1"
  port::Integer = 4223
  uid::String
  motorCurrent::typeof(1.0u"A")
  stepResolution::Integer
  microStepResolution::Integer
  gearRatio::Rational = 1//1
  allowedVoltageDrop
  stallThreshold::Integer = 20000
end
TinkerforgeSilentStepperParams(dict::Dict) = params_from_dict(TinkerforgeSilentStepperParams, dict)

Base.@kwdef mutable struct TinkerforgeSilentStepper <: StepperMotor
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimulatedGaussMeterParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  brick = missing
  stallSum = 0
end

function init(motor::TinkerforgeSilentStepper)
  @debug "Initializing Tinkerforge silent stepper unit with ID `$(motor.deviceID)`."

  tinkerforge = PyCall.pyimport("tinkerforge")
  ipConnection = tinkerforge.ip_connection.IPConnection()

  try
    ipConnection.connect(host, port)
  catch
    throw(ScannerConfigurationError("Could not establish a connection to the motor with the given parameters: "))
  end

  motor.brick = tinkerforge.brick_silent_stepper.BrickSilentStepper(uid(motor), ipConnection)

  # connect the all_data callback to catch relevant motor information and stop on an error
  # motor.brick.set_all_data_period(100)  # /ms, return values every x milliseconds
  # motor.brick.register_callback(tinkerforge.brick_silent_stepper.BrickSilentStepper.CALLBACK_ALL_DATA, allDataCallback)

  motor.brick.set_max_velocity(maxVelocity)  # /steps per s, max velocity in steps per second, depends on step_resolution
  motor.brick.set_speed_ramping(acceleration, deacceleration)  # steps per s^2, acceleration and deacceleration of the motor, 8000 steps per s in 10 s equals 800 steps per s^2
  motor.brick.set_motor_current(motorCurrent * 1000)  # /mA, sets the current to drive the motor.
  motor.brick.set_step_configuration(stepResolution, true)  # sets the defines stepResolution and activates Interpolation

  motor.present = true
end

checkDependencies(motor::TinkerforgeSilentStepper) = true

Base.close(motor::TinkerforgeSilentStepper) = nothing

# Note: This is commented out because I am not sure how to handle the callback. Using a global variable would only work with one device.
# function allDataCallback(current_velocity, current_position, remaining_step, stack_voltage, external_voltage,
#                     current_consumption)
#     if ssBrick.get_driver_status().stallguard_result > 300:
#         stallSum = stallSum + ssBrick.get_driver_status().stallguard_result
#         print("stallSum: "+str(stallSum))
#     if (stallSum > stallThreashold):
#         stopMotorInSec(0)

#         print("Motor is stalled, so it was stopped\n")
#         infoString = "The motor stopped at:\n\tcurrent velocity = \t" + str(
#             current_velocity) + " steps per second\n\tcurrent velocity = \t" + str(
#             sps2rpm(current_velocity, microStepResolution, 1)
#         ) + "in rpm\n\tcurrent position = \t" + str(
#             current_position) + " in no of steps:\n\tremaining steps = \t" + str(
#             remaining_step) + " in no of steps\n\tcurrent consumption = \t" + str(
#             current_consumption) + " in milli amperes."
#         print(infoString)
#         print(getMotorStatus())
#         stalledVelocity = current_velocity
#         stalledBool = True
#     if getMotorDirection() == "still":
#         stallSum = 0

host(motor::TinkerforgeSilentStepper) = motor.params.host
port(motor::TinkerforgeSilentStepper) = motor.params.port
uid(motor::TinkerforgeSilentStepper) = motor.params.uid
motorCurrent(motor::TinkerforgeSilentStepper) = motor.params.motorCurrent
stepResolution(motor::TinkerforgeSilentStepper) = motor.params.stepResolution
microStepResolution(motor::TinkerforgeSilentStepper) = motor.params.microStepResolution
gearRatio(motor::TinkerforgeSilentStepper) = motor.params.gearRatio
allowedVoltageDrop(motor::TinkerforgeSilentStepper) = motor.params.allowedVoltageDrop
stallThreshold(motor::TinkerforgeSilentStepper) = motor.params.stallThreshold

# Function that calculates the velocity in steps per second to get x rpm at the motor
rpm2sps(motor::TinkerforgeSilentStepper, rpm::Real) = rpm * 200 / 60 * microStepResolution(motor) * gearRatio(motor)
sps2rpm(motor::TinkerforgeSilentStepper, sps::Real) = sps * 60 / 200 * 1 / microStepResolution(motor) * 1 / gearRatio(motor)

function direction(motor::TinkerforgeSilentStepper)
  if motor.brick.get_current_velocity() > 0
    return MOTOR_FORWARD
  elseif motor.brick.get_current_velocity() < 0
    return MOTOR_BACKWARD
  else
    return MOTOR_STILL
  end
end

function emergencyBreak(motor::TinkerforgeSilentStepper)
  motor.brick.full_brake()
  motor.brick.disable()
end

function driveSteps(motor::TinkerforgeSilentStepper, numSteps::Integer, direction::MotorDirection)
  if direction == MOTOR_BACKWARD
    numSteps = -numSteps
  end

  # If the motor is not still: brake and drive afterwards
  if direction(motor) != MOTOR_STILL
    stop(motor, 0u"s")
  end

  # Enable the motor if not done already
  if !motor.brick.is_enabled()
    motor.brick.enable()
  end

  # Set the numper of steps the motor should drive. The sign defines the direction.
  motor.brick.set_steps(numSteps)
end

function drive(motor::TinkerforgeSilentStepper, direction::MotorDirection)
  if direction(motor) != MOTOR_STILL
    stop(motor, 0u"s")
  end

  if !motor.brick.is_enabled()
    motor.brick.enable()
  end

  if direction(motor) == MOTOR_BACKWARD
    motor.brick.drive_backward() # drives motor with set parameters till drive_forward or stop is called
  else
    motor.brick.drive_forward() # drives motor with set parameters till drive_backward or stop is called
  end
end

# Try to reach velocity while decreasing the acceleration and increasing the velocity
# Commented out because with it CI failed
#=function driveToVelocity(motor::TinkerforgeSilentStepper, rpm::typeof(u"1/s"), direction::MotorDirection)
  global stallThreashold
  stallThreashold = 10000
  if not (getMotorDirection() == "still"):
    stopMotorInSec(0)
  if not (ssBrick.is_enabled()):
    ssBrick.enable()  # enable the motor if not done already
    vSPS = rpm2sps(rpm, microStepResolution, gearRatio)
  while round(sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio)) < rpm:
  newAcceleration = 1000
  setupMotorParameters(vSPS, newAcceleration, vSPS / 2, motorCurrent, stepResolution)

  if sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio) > 500:
  newAcceleration = 50
  elif sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio) > 400:
  newAcceleration = 100
  elif sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio) > 300:
  newAcceleration = 200
  elif sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio) > 200:
  newAcceleration = 500
  elif sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio) > 100:
  newAcceleration = 1000
  else:
  newAcceleration = 1

  if (direction == "backward"):
  ssBrick.drive_backward()  # drives motor with set parameters till drive_forward or stop is called
  else:
  ssBrick.drive_forward()  # drives motor with set parameters till drive_backward or stop is called
  time.sleep(0.1)
  print(sps2rpm(ssBrick.get_current_velocity(),microStepResolution,gearRatio))
  if getMotorDirection() == "still":
  break

  if (direction == "backward"):
  ssBrick.drive_backward()  # drives motor with set parameters till drive_forward or stop is called
  else:
  ssBrick.drive_forward()  # drives motor with set parameters till drive_backward or stop is called
  stallThreashold = 20000
end=#

function stop(motor::TinkerforgeSilentStepper, delay::typeof(1.0u"s") = 0u"s")
  sleep(ustrip(u"s", delay))  # Delay the stop maneuver
  motor.brick.stop()  # Stop the motor with the deacceleration set in the motor parameters

  # Wait till the motor has actually stopped before disabling the driver
  sleep((motor.brick.get_max_velocity() / motor.brick.get_speed_ramping().deacceleration) + 0.5)  # steps per second / steps per second^2 = second
  motor.brick.disable()

  # global stallSum
  # stallSum=0
end

function progress(motor::TinkerforgeSilentStepper)
  setSteps = motor.brick.get_steps()
  remainingSteps = motor.brick.get_remaining_steps()
  donePercentage = 100 * remainingSteps / setSteps
  return donePercentage
end


text:
using Dojo

# Parameters
origin = Origin{Float64}()
pbody = Sphere(0.5, 10.0)
cbody = Box(1.0, 1.0, 2.0, 1.0)
joint1 = JointConstraint(Floating(origin, pbody))
joint2 = JointConstraint(Floating(origin, cbody))

bodies = [pbody, cbody]
joints = [joint1, joint2]

side = 1.0
corners = [
            [[ side / 2.0;  side / 2.0; -side]]
            [[ side / 2.0; -side / 2.0; -side]]
            [[-side / 2.0;  side / 2.0; -side]]
            [[-side / 2.0; -side / 2.0; -side]]
            [[ side / 2.0;  side / 2.0;  side]]
            [[ side / 2.0; -side / 2.0;  side]]
            [[-side / 2.0;  side / 2.0;  side]]
            [[-side / 2.0; -side / 2.0;  side]]
        ]

normal = [[0.0, 0.0, 1.0] for i = 1:8]
contact_radius = [0.0 for i = 1:8]
friction_coefficient = 0.5 * ones(8)

sphere_contact = contact_constraint(pbody, [0.0; 0.0; 1.0],
    friction_coefficient=0.5,
    contact_origin=SA[0.0; 0.0; 0.0],
    contact_radius=0.5)

box_contacts = contact_constraint(cbody, normal,
    friction_coefficient=friction_coefficient,
    contact_origins=corners,
    contact_radius=contact_radius,
    contact_type=:nonlinear)

collision = SphereBoxCollision{Float64,2,3,6}(
    szeros(3),
    side,
    side,
    2 * side,
    0.5,
)

friction_parameterization = SA{Float64}[
    1.0  0.0
    0.0  1.0
]
body_body_contact = NonlinearContact{Float64,8}(0.5, friction_parameterization, collision)

contacts = [
        # corner_contact1,
        # corner_contact2,
        # corner_contact3,
        # corner_contact4,
        # corner_contact5,
        # corner_contact6,
        # corner_contact7,
        # corner_contact8,
        box_contacts...,
        sphere_contact,
        ContactConstraint((body_body_contact, pbody.id, cbody.id), name=:body_body)
]

mech = Mechanism(origin, bodies, joints, contacts,
            gravity=0.0 * -9.81,
            timestep=0.01)

mech.bodies[1].state.x2 = [2.0, 0.1, 2.0]
mech.bodies[1].state.v15 = [-1.0, 0.1, 0.0]

mech.bodies[2].state.x2 = [0.0, 0.0, side]
# q = rand(4)
# q ./= norm(q)
mech.bodies[2].state.q2 = RotZ(0.25 * π)

function ctrl!(m, k)
    u = [
          0.0; 0.0; 0.0 * -mech.bodies[1].mass * 9.81 * mech.timestep; 0.0; 0.0; 0.0; 0.0; 0.0; -1.0 * mech.bodies[2].mass * 9.81 * mech.timestep; 0.0; 0.0; 0.0;
    ]
    set_input!(m, u)
    return nothing
end

storage = simulate!(mech, 5.0,
    ctrl!,
    # verbose=true,
    record=true,
    opts=SolverOptions(verbose=true, rtol=1.0e-5, btol=1.0e-5));

mech.bodies[1].name
mech.bodies[2].name
vis = Visualizer()
render(vis)
visualize(mech, storage,
    vis=vis)
mech.bodies[2].state.vsol
x31 = next_position(mech.bodies[1].state.x2, mech.bodies[1].state.vsol[1], mech.timestep)
x32 = next_position(mech.bodies[2].state.x2, mech.bodies[2].state.vsol[2], mech.timestep)
q31 = next_orientation(mech.bodies[1].state.q2, mech.bodies[1].state.ϕsol[1], mech.timestep)
q32 = next_orientation(mech.bodies[2].state.q2, mech.bodies[2].state.ϕsol[2], mech.timestep)

distance(collision, x31, q31, x32, q32)
∂distance∂x(:parent, collision, x31, q31, x32, q32)
∂distance∂x(:child, collision, x31, q31, x32, q32)

pp = contact_point(:parent, collision, x31, q31, x32, q32)
pc = contact_point(:child, collision, x31, q31, x32, q32)

cop = contact_point_origin(x31, q31, collision.origin_sphere)
coc = contact_point_box(cop, x32, q32, collision.width, collision.depth, collision.height)


text: abstract type BaseStructureLearning end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

"""
    write_to_file(
        model::Model,
        filename::String;
        format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
    )

Write the JuMP model `model` to `filename` in the format `format`.

If the filename ends in `.gz`, it will be compressed using Gzip.
If the filename ends in `.bz2`, it will be compressed using BZip2.
"""
function write_to_file(
    model::Model,
    filename::String;
    format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC,
)
    dest = MOI.FileFormats.Model(format = format, filename = filename)
    # We add a `full_bridge_optimizer` here because MOI.FileFormats models may not
    # support all constraint types in a JuMP model.
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    MOI.copy_to(bridged_dest, model)
    # `dest` will contain the underlying model, with constraints bridged if
    # necessary.
    MOI.write_to_file(dest, filename)
    return
end

"""
    Base.write(
        io::IO,
        model::Model;
        format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_MOF
    )

Write the JuMP model `model` to `io` in the format `format`.
"""
function Base.write(
    io::IO,
    model::Model;
    format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_MOF,
)
    if format == MOI.FileFormats.FORMAT_AUTOMATIC
        error("Unable to infer the file format from an IO stream.")
    end
    dest = MOI.FileFormats.Model(format = format)
    # We add a `full_bridge_optimizer` here because MOI.FileFormats models may not
    # support all constraint types in a JuMP model.
    bridged_dest = MOI.Bridges.full_bridge_optimizer(dest, Float64)
    MOI.copy_to(bridged_dest, model)
    # `dest` will contain the underlying model, with constraints bridged if
    # necessary.
    write(io, dest)
    return
end

"""
    read_from_file(
        filename::String;
        format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC
    )

Return a JuMP model read from `filename` in the format `format`.

If the filename ends in `.gz`, it will be uncompressed using Gzip.
If the filename ends in `.bz2`, it will be uncompressed using BZip2.
"""
function read_from_file(
    filename::String;
    format::MOI.FileFormats.FileFormat = MOI.FileFormats.FORMAT_AUTOMATIC,
)
    src = MOI.FileFormats.Model(format = format, filename = filename)
    MOI.read_from_file(src, filename)
    model = Model()
    MOI.copy_to(model, src)
    return model
end

"""
    Base.read(io::IO, ::Type{Model}; format::MOI.FileFormats.FileFormat)

Return a JuMP model read from `io` in the format `format`.
"""
function Base.read(io::IO, ::Type{Model}; format::MOI.FileFormats.FileFormat)
    if format == MOI.FileFormats.FORMAT_AUTOMATIC
        error("Unable to infer the file format from an IO stream.")
    end
    src = MOI.FileFormats.Model(format = format)
    read!(io, src)
    model = Model()
    MOI.copy_to(model, src)
    return model
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Parameters
polynomial(x::Vector, P::Int) = hcat([x.^p for p in 1:P]...)
polynomial(x::AbstractRange, P::Int)= polynomial(collect(x),P)
polynomial(x::Vector, P::Int, f::Function) = hcat([f.(p*x) for p in 1:P]...)
polynomial(x::AbstractRange, P::Int, f::Function) = polynomial(collect(x), P, f)

function polyfit(x,y;x_pred=x, P::Int=3)
    X = polynomial(x,P)
    X_pred=polynomial(x_pred,P)
    β̂ = multiple_regression(X,y)
    y_pred = insert_ones(X_pred) * β̂
end

@with_kw mutable struct Spline{T<:Number}
    xmin::T = -1.0
    xmax::T = 1.0
    K::Int = 5
end

function (s::Spline)(x::Vector, y::Vector, x_pred::Vector;
          splinematrix::Function=spline_matrix)
    @unpack xmin, xmax, K = s
    X = splinematrix(x,K,xmin,xmax)
    β = (X'X)\X'y
    X_pred = splinematrix(x_pred,K,xmin,xmax)
    ypred = X_pred*β
end

(s::Spline)(x::Vector, y::Vector, x_pred::AbstractRange;kwargs...) = s(x,y,collect(x_pred);kwargs...)

function spline_matrix(x::Vector,K::Int,xmin,xmax;use_x=false)
    n = length(x)
    Knots = use_x ? x : range(xmin,stop=xmax,length=K)
    K = length(Knots)
    X = zeros(n,K+4)
    X[:,1] .= 1
    X[:,2] .= x
    X[:,3] .= x.^2
    X[:,4] .= x.^3
    for j in 1:K
        X[:,j+4] = @. max(x - Knots[j], 0)^3
    end
    return X
end

function natural_spline_matrix(x::Vector,K::Int,xmin,xmax)
    n = length(x)
    Knots  = range(xmin,stop=xmax,length=K)
    X = zeros(n,K)
    X[:,1] .= 1
    X[:,2] .= x
    d(x::Vector,a::Number,b::Number) =  @. (max(x-a,0)^3 - max(x-b,0)^3) / (b-a)
    for j in 1:K-2
        X[:,j+2] = d(x,Knots[j],Knots[end]) - d(x,Knots[end-1],Knots[end])
    end
    return X
end

struct SmoothingSpline
    x::Vector
end

function (s::SmoothingSpline)(y::Vector,x_pred::Union{Vector,AbstractRange};λ=0)
    x = s.x
    X = smoothing_spline_matrix(x,x)
    G = green_silverman(x)
    γ = (X'X + λ*G)\X'y
    X_pred = smoothing_spline_matrix(x_pred,x)
    ypred = X_pred*γ
end

function smoothing_spline_matrix(x::Union{Vector,AbstractRange},Knots::Union{Vector, AbstractRange})
    n = length(x)
    K = length(Knots)
    X = zeros(n,K)
    X[:,1] .= 1
    X[:,2] .= x
    d(x::Union{Vector,AbstractRange},a::Number,b::Number) =  @. (max(x-a,0)^3 - max(x-b,0)^3) / (b-a)
    for j in 1:K-2
        X[:,j+2] = d(x,Knots[j],Knots[end]) - d(x,Knots[end-1],Knots[end])
    end
    return X
end

function green_silverman(x::Vector)
    n = length(x)
    g = zeros(n,n)
    for i in 3:n ,j in i:n
        g[j,i] =(12(x[end]-x[end-1])*(x[end-1]-x[j-2])*(x[end-1]-x[i-2])+
                (x[end-1]-x[j-2])^2*(12x[end-1]+6x[j-2]-18x[i-2]))/
                (x[end]-x[i-2])/(x[end]-x[j-2])
        g[i,j] = g[j,i]
    end
    return g
end

function cv_ss_fast(X::AbstractArray,y::Vector,G::AbstractArray,K::Int;λ=0)
    n = length(y); m = round(Int,n/K)
    H = X*((X'X + λ*G)\X')
    e = (I(n) - H)*y
    S = 0
    @views for j in 1:K
        test = j*m-m+1:j*m;
        err = (I(m)-H[test,test]) \ e[test]
        S += dot(err,err)
    end
    return S/n,tr(H)
end

#epanechnikovカーネル
epanechnikov(x,y,λ) = norm(x-y)/λ |> c -> max(0.75*(1-c^2),0)
#nadaraya-watson推定量
function nadaraya_watson_estimator(x_observed::Vector,y_observed::Vector,kk::Function,x_pred::TP;λ=1.0) where {TP<:Number}
    n = length(y_observed)
    S = 0; T = 0;
    for i in 1:n
        S += kk(x_observed[i],x_pred,λ)*y_observed[i]
        T += kk(x_observed[i],x_pred,λ)
    end
    result = T == 0 ? 0 : S/T
    return result
end

function nadaraya_watson_estimator(x_observed::Vector,y_observed::Vector,kk::Function,x_pred::Union{Vector,AbstractRange};λ=1.0)
    y_pred = [nadaraya_watson_estimator(x_observed,y_observed,kk,xp; λ) for xp in x_pred]
end

function local_regression(x::Matrix{T},
                          y::Vector{T},
                          kernel::Function;
                          x_pred=x, λ=1) where {T<:Number}
    N,_ = size(x)
    @assert N==length(y)
    X = insert_ones(x)
    m,_ = size(x_pred)
    y_pred = Vector{T}(undef,m)
    for j in 1:m
        W = [kernel(x[i,:],x_pred[j,:],λ) for i in 1:N] |> diagm
        β̂ = (X'*W*X)\X'*W*y
        y_pred[j] = dot(insert_ones(x_pred[j,:]),β̂)
    end
    return y_pred
end
function local_regression(x::Vector{T}, y::Vector{T},kernel::Function;x_pred=x,λ=1) where {T<:Number}
    return local_regression(x[:,:],y,kernel;x_pred=x_pred[:,:],λ=λ)
end


text:
# This file is pat of RLinearAlgebra.jl
# 1. Specifies type
# 2. Implements sample function
# 3. Exports Type

# using Random

"""
    LinSysVecRowHopRandCyclic <: LinSysVecRowSelect

A mutable structure that specifies a cycling through the rows of a linear system, where the
cycling order is determined randomly once the current cycling order has been used `hop`
number of times. The solver randomly chooses the cycling order whenever necessary.

# Fields
- `order::Union{Vector{Int64},Nothing}`
- `hop::Int64`

# Constructors
- `LinSysVecRowHopRandCyclic()` defaults to setting the `order` to `nothing` and the `hop`
    to `5` (i.e., each ordering is used five times before sampling a new ordering).
- `LinSysVecRowHopRandCyclic(hop::Int64)` defaults to setting the `order` to `nothing` and
    the `hop` to whatever is specified by the argument.
"""
mutable struct LinSysVecRowHopRandCyclic <: LinSysVecRowSelect
    order::Union{Vector{Int64},Nothing}
    hop::Int64
end
LinSysVecRowHopRandCyclic() = LinSysVecRowHopRandCyclic(nothing, 5)
LinSysVecRowHopRandCyclic(hop::Int64) = LinSysVecRowHopRandCyclic(nothing, hop)

# Common sample interface for linear systems
function sample(
    type::LinSysVecRowHopRandCyclic,
    A::AbstractArray,
    b::AbstractVector,
    x::AbstractVector,
    iter::Int64
)

    # Determine whether the current ordering has been used `hop` times
    if mod(iter, 1:(length(b) * type.hop)) == 1
        type.order = randperm(length(b))
    end

    eqn_ind = type.order[mod(iter, 1:length(b))]
    return A[eqn_ind, :], b[eqn_ind]
end

#export LinSysVecRowHopRandCyclic


text:
using DSGE
using HDF5

path = dirname(@__FILE__)

# Test in generic case
# Global minimum at (a, a^2)
a = 1
b = 100
function rosenbrock_csminwel(x::Vector)
    return (a-x[1])^2.0 + b*(x[2]-x[1]^2.0)^2.0
end

# this is the actual gradient of the Rosenbrock function
function rosenbrock_grad(x::Array)
    dr = similar(x)
    dr[1] = -2*(a-x[1]) - 2b*(2x[1])*(x[2]-x[1]^2.0)
    dr[2] = 2b*(x[2]-x[1]^2.0)
    badg = false
    return dr, badg
end

# A really bad guess
x_init = [10.0, -9.0]
x_expected = [a, a^2.0]

@testset "Check csminwel gradient calculation for correctness" begin
    res_real_grad,  = csminwel(rosenbrock_csminwel, rosenbrock_grad, x_init)
    @test x_expected ≈ res_real_grad.minimizer
    res_numeric_grad,  = csminwel(rosenbrock_csminwel, x_init)
    @test x_expected ≈ res_numeric_grad.minimizer atol=1.0e-8
end

nothing


text:
function fit(m::MultinomialNB, X::Matrix{Int64}, y::Vector{C}) where C
    ensure_data_size(X, y)
    for j=1:size(X, 2)
        c = y[j]
        m.c_counts[c] += 1
        m.x_counts[c] .+= X[:, j]
        m.x_totals += X[:, j]
        m.n_obs += 1
    end
    return m
end

"""Calculate log P(x|C)"""
function logprob_x_given_c(m::MultinomialNB, x::Vector{Int64}, c::C) where C
    x_priors_for_c = m.x_counts[c] ./ sum(m.x_counts[c])
    x_probs_given_c = x_priors_for_c .^ x
    logprob = sum(log(x_probs_given_c))
    return logprob
end

"""Calculate log P(x|C)"""
function logprob_x_given_c(m::MultinomialNB, X::Matrix{Int64}, c::C) where C
    x_priors_for_c = m.x_counts[c] ./ sum(m.x_counts[c])
    x_probs_given_c = x_priors_for_c .^ X
    logprob = sum(log.(x_probs_given_c), dims=1)
    return dropdims(logprob, dims=1)
end


text:
#=
    reference
       http://szudzik.com/ElegantPairing.pdf
               x    .    .    .
           z  : 0 |  1 |  2 | 3   ...
         -----:---|----|----|----
     y     0  : 0 |  2 |  6 | 12
     .     1  : 1 |  3 |  7 | 13
     .     2  : 4 |  5 |  8 | 14
     .     3  : 9 | 10 | 11 | 15
          ...
    for x,y,z >= 0
    pair(x,y) =
        if x == max(x,y)
            x*x + x + y
        else
            y*y + x
        end

    unpair(z) =
        if z - isqrt(z)^2 < isqrt(z)
            (z - isqrt(z)^2, isqrt(z))
        else
            (isqrt(z), z - isqrt(z)^2 - isqrt(z))
        end
=#

const SignedInt   = Union{Int8, Int16, Int32, Int64, Int128}

# elegantpair(pairablemax(T), pairablemax(T)) == pairedmax(T)

pairablemax(::Type{Int128}) = Int128(13_043_817_825_332_782_211)
pairablemax(::Type{Int64})  = Int64(3_037_000_498)
pairablemax(::Type{Int32})  = Int32(46_339)
pairablemax(::Type{Int16})  = Int16(180)
pairablemax(::Type{Int8})   = Int8(10)

pairedmax(::Type{Int128}) = Int128(170_141_183_460_469_231_722_567_801_800_623_612_943)
pairedmax(::Type{Int64})  = Int128(9_223_372_030_926_249_000)
pairedmax(::Type{Int32})  = Int128(2_147_395_599)
pairedmax(::Type{Int16})  = Int128(32_760)
pairedmax(::Type{Int8})   = Int128(120)

function type_for_pair_indexing(xmax::T, ymax::T) where T<:SignedInt
    if min(xmax, ymax) <= zero(T) || (T===Int128 && max(xmax, ymax) > pairablemax(Int128))
         throw(ErrorException("($xmax,$ymax) is outside of domain for paired indexing"))
    end

    x = Int128(xmax)
    y = Int128(ymax)
    xy = elegantpair(x,y)

    if xy <= pairablemax(Int8)
        itype = Int8
    elseif xy <= pairablemax(Int16)
        itype = Int16
    elseif xy <= pairablemax(Int32)
        itype = Int32
    elseif xy <= pairablemax(Int64)
         itype = Int64
    elseif xy <= pairablemax(Int128)
         itype = Int128
    else
         throw(ErrorException("($xmax,$ymax) is outside of domain for paired indexing"))
    end
    return itype
end

# x, y >= 0
function elegantpair(x::T, y::T) where T<:SignedInt
    mx2 = mx = max(x,y)
    mx2 *= mx
    mx2 += x
    return mx2 + (x===mx ? y : zero(T))
end

# z >= 0
function elegantunpair(z::T) where T<:SignedInt
    isz  = isqrt(z)
    isz2 = isz*isz
    x = z - isz2
    if x < isz
        y = isz
    else
        y = x > isz ? x - isz : x
        x = isz
    end
    return x, y
end

function safe_elegantpair(x::T, y::T) where T<:SignedInt
    mx = max(x,y)
    (x < 0 || y < 0 || mx > pairablemax(T)) && throw(DomainError())
    return elegantpair(x, y, mx)
end

function safe_elegantunpair(z::T) where T<:SignedInt
    z < zero(T) && throw(DomainError())
    return elegantunpair(z)
end

# x, y >= 0, mx = max(x,y)
function elegantpair(x::T, y::T, mx::T) where T<:SignedInt
    mx2 = mx
    mx2 *= mx
    mx2 += x
    return mx2 + (x===mx ? y : zero(T))
end

#=
    transform to 1-based indexing
    elegantpair(lo, hi) < elegantpair(hi, lo)
       where lo, hi > 0 and lo < hi
    assume that more often nrows > ncols
       so pair(colidx, rowidx)
=#
@inline two(::Type{T}) where T<:Integer = one(T)+one(T)

#=
pair(x::T, y::T) where T<:Signed = elegantpair(y,x) - two(T)
function pair(x::T, y::T) where T<:Signed
    mx2 = mx = max(x,y)
    mx2 *= mx
    mx2 += y
    mx2 += (y===mx ? x : zero(T))
    mx2 -= two(T)
    return mx2
end
One one machine, he version below has better consistant median time.
=#

function pair(x::T, y::T) where T<:Signed
    mx2 = mx = max(x,y)
    mx2 = mx2*mx + y
    mx2 += (y===mx ? x : zero(T))
    mx2 -= two(T)
    return mx2
end

function unpair(z::T) where T<:Signed
    y, x = elegantunpair(z + two(T))
    return x, y
end