text:
using LinearAlgebra
using LinearAlgebra: chkstride1, BlasInt
using LinearAlgebra.LAPACK: liblapack, chklapackerror
using LinearAlgebra.BLAS: @blasfunc
using Base: has_offset_axes

for (gesvj, elty, relty) in
    ((:dgesvj_, :Float64, :Float64),
     (:sgesvj_, :Float32, :Float32),
     (:zgesvj_, :ComplexF64, :Float64),
     (:cgesvj_, :ComplexF32, :Float32))
    @eval begin
#       SUBROUTINE ZGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
#                          LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
#
#       .. Scalar Arguments ..
#       INTEGER            INFO, LDA, LDV, LWORK, LRWORK, M, MV, N
#       CHARACTER*1        JOBA, JOBU, JOBV
#       ..
#       .. Array Arguments ..
#       COMPLEX*16         A( LDA, * ),  V( LDV, * ), CWORK( LWORK )
#       DOUBLE PRECISION   RWORK( LRWORK ),  SVA( N )
        function gesvj!(joba::AbstractChar, jobu::AbstractChar,
                        jobv::AbstractChar, A::AbstractMatrix{$elty};
                        Vinit = similar(A, $elty, (0, size(A,2))),
                        ctol=zero($relty))
            @assert !has_offset_axes(A)
            @assert !has_offset_axes(Vinit)
            chkstride1(A)
            m, n   = size(A)
            if m < n
                throw(ArgumentError("matrix A must be tall or square"))
            end
            minmn  = min(m, n)
            S      = similar(A, $relty, minmn)
            if (jobv != 'N') && (size(Vinit,2) != n)
                throw(DimensionMismatch("matrix Vinit must have n columns"))
            end
            mv = 0
            if jobv ∈ ('A','N')
                mv = size(Vinit,1)
                V = Vinit
            elseif jobv ∈ ('V','J')
                V = similar(A, $elty, (n,n))
            end
            cmplx  = eltype(A) <: Complex
            if cmplx
                lcwork = -1
                cwork = Vector{$elty}(undef, 1)
                lwork = -1
                work   = Vector{$relty}(undef, 1)
            else
                lwork  = max(6,m+n)
                work   = Vector{$relty}(undef, lwork)
            end
            if jobu == 'C'
                if (ctol < 1)
                    throw(ArgumentError("appropriate ctol must be provided for jobu = 'C'"))
                end
                work[1] = ctol
            end
            info = Ref{BlasInt}()
            if cmplx
                for i = 1:2 # first call returns workspace requirements
                    ccall((@blasfunc($gesvj), liblapack), Cvoid,
                          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobA/U/V
                           Ref{BlasInt}, Ref{BlasInt}, # m,n
                           Ptr{$elty}, Ref{BlasInt}, # A,lda
                           Ptr{$relty}, Ref{BlasInt}, # SVA, mv
                           Ptr{$elty}, Ref{BlasInt}, # V, ldv
                           Ptr{$elty}, Ref{BlasInt}, # cwork, lcwork
                           Ptr{$relty}, Ref{BlasInt}, # rwork, lwork
                           Ptr{BlasInt}), # info
                          joba, jobu, jobv, m, n, A, max(1,stride(A,2)), S, mv, V, max(1,stride(V,2)),
                          cwork, lcwork, work, lwork, info)
                    chklapackerror(info[])
                    if (i == 1)
                        lwork = BlasInt(work[1])
                        resize!(work, lwork)
                        lcwork = BlasInt(real(cwork[1]))
                        resize!(cwork, lcwork)
                    end
                end
            else
                ccall((@blasfunc($gesvj), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobA/U/V
                       Ref{BlasInt}, Ref{BlasInt}, # m,n
                       Ptr{$elty}, Ref{BlasInt}, # A,lda
                       Ptr{$relty}, Ref{BlasInt}, # SVA, mv
                       Ptr{$elty}, Ref{BlasInt}, # V, ldv
                       Ptr{$elty}, Ref{BlasInt}, # work, lwork
                       Ptr{BlasInt}), # info
                      joba, jobu, jobv, m, n, A, max(1,stride(A,2)), S, mv, V, max(1,stride(V,2)),
                      work, lwork, info)
                chklapackerror(info[])
            end

            scale = work[1]
            ranka = Int(work[2])
            extra = work[3:6]

            if jobu ∈ ('U','F','C')
                U = A[:,1:ranka]
            else
                U = similar(A, (m,0))
            end
            if jobv ∉ ('V','J','A')
                    V = similar(A, (n,0))
            end
            return U, S, V, scale, ranka, extra
        end # function
    end # eval
end # loop over types

"""
    gesvj!(joba, jobu, jobv, A; Vinit = similar(A, (0,n))) -> (U, S, V, scale, ranka, extras)

Compute the singular value decomposition of `A`, `A = U * (scale * Diagonal(S)) * V'`.

Uses the single-sided Jacobi scheme from LAPACK.

`ranka` is the number of significant singular values computed.

`extras` is `rwork[3:6]` returned from the LAPACK routine.

See the LAPACK documentation (`dgesvj` etc.) for details.
"""
function gesvj! end

for (gejsv, elty, relty) in
    ((:dgejsv_, :Float64, :Float64),
     (:sgejsv_, :Float32, :Float32),
     (:zgejsv_, :ComplexF64, :Float64),
     (:cgejsv_, :ComplexF32, :Float32))
    @eval begin
#     SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
#                         M, N, A, LDA, SVA, U, LDU, V, LDV,
#                         CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
#
#     .. Scalar Arguments ..
#     IMPLICIT    NONE
#     INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
#     ..
#     .. Array Arguments ..
#     COMPLEX*16     A( LDA, * ),  U( LDU, * ), V( LDV, * ), CWORK( LWORK )
#     DOUBLE PRECISION   SVA( N ), RWORK( LRWORK )
#     INTEGER     IWORK( * )
#     CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV

        function gejsv!(joba::AbstractChar, jobu::AbstractChar,
                        jobv::AbstractChar, A::AbstractMatrix{$elty};
                        jobr::AbstractChar = 'R',
                        jobt::AbstractChar = 'N', jobp::AbstractChar = 'N')
            @assert !has_offset_axes(A)
            chkstride1(A)
            m, n   = size(A)
            if m < n
                throw(ArgumentError("matrix A must be tall or square"))
            end
            minmn  = min(m, n)
            S      = similar(A, $relty, minmn)
            if jobv == 'N'
                V = similar(A, $elty, (n,0))
            else
                V = similar(A, $elty, (n,n))
            end
            if jobu == 'U'
                U = similar(A, $elty, (m,n))
            elseif jobu == 'F'
                U = similar(A, $elty, (m,m))
            elseif jobu == 'N'
                U = similar(A, $elty, (m,0))
            else
                # use as workspace
                U = similar(A, $elty, (m,n))
            end
            cmplx  = eltype(A) <: Complex
            if cmplx
                lcwork = BlasInt(-1)
                cwork = Vector{$elty}(undef, 7)
                liwork = 1
                lwork  = BlasInt(-1)
            else
                liwork = m+3*n;
                lwork = max(n+m*32,2*m+n,6*n+2*n*n)
            end
            work   = Vector{$relty}(undef, max(7,lwork))
            iwork = Vector{BlasInt}(undef, liwork)
            info   = Ref{BlasInt}()
            if cmplx
                for i = 1:2 # first call returns workspace requirements
                    ccall((@blasfunc($gejsv), liblapack), Cvoid,
                          (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobA/U/V
                           Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobR/T/P
                           Ref{BlasInt}, Ref{BlasInt}, # m,n
                           Ptr{$elty}, Ref{BlasInt}, # A,lda
                           Ptr{$relty}, # SVA,
                           Ptr{$elty}, Ref{BlasInt}, # U, ldu
                           Ptr{$elty}, Ref{BlasInt}, # V, ldv
                           Ptr{$elty}, Ref{BlasInt}, # cwork, lwork
                           Ptr{$relty}, Ref{BlasInt}, # rwork, lrwork
                           Ptr{BlasInt}, # iwork
                           Ptr{BlasInt}), # info
                          joba, jobu, jobv, jobr, jobt, jobp, m, n,
                          A, max(1,stride(A,2)), S,
                          U, max(1,stride(U,2)), V, max(1,stride(V,2)),
                          cwork, lcwork, work, lwork, iwork, info)
                    chklapackerror(info[])
                    if (i == 1)
                        lwork = BlasInt(work[1])
                        resize!(work, lwork)
                        lcwork = BlasInt(real(cwork[1]))
                        resize!(cwork, lcwork)
                        resize!(iwork, iwork[1])
                    end
                end
            else
                ccall((@blasfunc($gejsv), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobA/U/V
                       Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, # jobR/T/P
                       Ref{BlasInt}, Ref{BlasInt}, # m,n
                       Ptr{$elty}, Ref{BlasInt}, # A,lda
                       Ptr{$relty},  # SVA,
                       Ptr{$elty}, Ref{BlasInt}, # U, ldu
                       Ptr{$elty}, Ref{BlasInt}, # V, ldv
                       Ptr{$elty}, Ref{BlasInt}, # work, lwork
                       Ptr{BlasInt}, # iwork
                       Ptr{BlasInt}), # info
                      joba, jobu, jobv, jobr, jobt, jobp, m, n,
                      A, max(1,stride(A,2)), S,
                      U, max(1,stride(U,2)), V, max(1,stride(V,2)),
                      work, lwork, iwork, info)
                chklapackerror(info[])
            end
            if iwork[3] != 0
                @warn "loss of accuracy: some column norms were denormals"
            end
            scales = (work[1], work[2])
            if joba ∈ ('E','G')
                sconda = work[3]
            else
                # absurd value as a passive-aggressive warning
                sconda = -one($relty)
            end
            extra = work[4:7]
            return U,S,V,scales,sconda,extra
        end # function
    end # eval
end # loop over types

"""
    gejsv!(joba, jobu, jobv, A; jobr='R', jobt='N', jobp='N') -> (U, S, V, scales, sconda, extras)

Finds the singular value decomposition of `A`,
`A = U * (σ * Diagonal(S)) * V'`, using the
preconditioned Jacobi algorithm implemented in LAPACK.

`scales` is a 2-tuple specifying a scale factor for `S`
to avoid over/underflow: `σ=scales[2]/scales[1]` in the above formula.
(`rwork[1:2]` in the LAPACK result).

`sconda` is an estimate of equilibrated condition number (`rwork[3]` in the LAPACK result) or -1 if not requested via `joba`.

`extras` are `rwork[4:7]` from the LAPACK result.

Briefly, `jobr='N'` allows possibly unreliably extreme singular values;
`jobt='T'` allows transposition of `A` (if square), to speed convergence;
`jobp='P'` attempts to quash denormals for efficiency.

See the LAPACK documentation (`dgejsv` etc.) for details.
"""
function gejsv! end


text:
using GridCreation

@testset "zernike zero separation" begin
    for i in 1:100
        sep = zernike_zero_separation(i)
        @test 1/sep ≈ Float64(i)^2
    end
end

@testset "zernike fringe and nm round trip" begin
    for i in 1:100
        n, m = zernike_fringe_to_nm(i)
        j = zernike_nm_to_fringe(n,m)
        @test j == i
    end
end

@testset "zernike ANSI round trips" begin
    for i in 1:100
        n, m = zernike_ansi_j_to_nm(i)
        j = zernike_nm_to_ansi_j(n,m)
        @test j == i
    end
end

@testset "zernike gives correct values" begin
    for i in 1:100
        # a visual test proves the correctness of the zernike function in
        # general.  The values could only really be tested in a regression
        # fashion, generally, which is undesireable and assumes another
        # implementation is de-facto correct.

        # the test here, therefore, verifies only the analytic value of one
        # at the edge for the maxima of the azumithal term.
        n, m = zernike_noll_to_nm(i)
        out = zernike(n, m, [1.], [0.], norm=false)
        @test (out ≈ [1.]) || (out ≈ [0.])
        out = zernike(n, m, [1.], [0], norm=true)
        norm = [zernike_norm(n,m)]
        @test (out ≈ norm) || (out ≈ [0.])
    end
end

@testset "zernike series contains proper modes" begin
    x, y = mkCartVecs(1/8, 16)
    r, t = cartVecsToPolarGrid(x,y)
    modes = 1:11
    nms = zernike_noll_to_nm.(modes)
    out = zernike_series(nms, r, t)
    for mode in modes
        n, m = zernike_noll_to_nm(mode)
        truth = zernike(n, m, r, t)
        @test out[:, :, mode] ≈ truth
    end
end

@testset "zernike sum contains proper modes" begin
    x, y = mkCartVecs(1/8, 16)
    r, t = cartVecsToPolarGrid(x,y)
    nms = zernike_noll_to_nm.(1:11)
    weights = zeros(11)
    weights[5]=1.
    out = zernike_sum(nms, weights, r, t)
    truth = zernike(2, -2, r, t)
    @test out ≈ truth
end

ext: jl


text:
# generate summarizers for each term and collect into a Vector to save on schema
# metadata.

"""
    make_summarizer(term::T, types::Dict{Symbol,Type}) where T<:Terms.Term

Create a `Summarizer` for a term, given the types of the data.  If no summaries
are needed, this method should return an empty array `[]`.
"""
make_summarizer(term::Terms.Term, types::Dict) = []
make_summarizer(term::Terms.Eval, types::Dict) =
    is_categorical(types[term.name]) ?
    CategoricalSummarizer{types[term.name]}(term.name) :
    []
make_summarizer(term::Terms.Interaction, types::Dict) =
    mapreduce(t -> make_summarizer(t, types), vcat, [], term.terms)


text:
Base.eltype(::AbstractDataGraph{T}) where {T} = T

## Edge access

"""
    has_vertex(g, label)

Check whether a vertex with label `label` exists by calling [`haskey(g, label)`](@ref).
"""
function Graphs.has_vertex(g::AbstractDataGraph{T,VL}, label::VL) where {T,VL}
    return haskey(g, label)
end

"""
    has_edge(g, label_s, label_d)

Check whether an edge exists between the vertices with labels `label_s` and `label_d`.
"""
function Graphs.has_edge(g::AbstractDataGraph{T,VL}, label_s::VL, label_d::VL) where {T,VL}
    if has_vertex(g, label_s) && has_vertex(g, label_d)
        s = get_vertex(g, label_s)
        d = get_vertex(g, label_d)
        return has_edge(g, s, d)
    else
        return false
    end
end

## Vertex and edge modification

"""
    add_edge!(g::DataDiGraph, label_s, label_d, data)

Add an edge to `g` by specifying the labels `label_s` and `label_d` of both vertices, along with the associated `data`.
"""
function Graphs.add_edge!(
    g::AbstractDataGraph{T,VL,VD,ED}, label_s::VL, label_d::VL, data::ED
) where {T,VL,VD,ED}
    s, d = get_vertex(g, label_s), get_vertex(g, label_d)
    return add_edge!(g, s, d, data)
end

"""
    add_vertex!(g, label)

Shortcut for `add_vertex(g, label, nothing)`, to use in graphs where vertices have no metadata.
"""
function Graphs.add_vertex!(
    g::AbstractDataGraph{T,VL,VD,ED}, label::VL
) where {T,VL,VD<:Nothing,ED}
    return add_vertex!(g, label, nothing)
end

"""
    add_edge!(g, s, d)

Shortcut for `add_edge!(g, s, d, nothing)`, to use in graphs where edges have no metadata.
"""
function Graphs.add_edge!(
    g::AbstractDataGraph{T,VL,VD,ED}, s::Integer, d::Integer
) where {T,VL,VD,ED<:Nothing}
    return add_edge!(g, s, d, nothing)
end

"""
    add_edge!(g, label_s, label_d)

Shortcut for `add_edge!(g, label_s, label_d, nothing)`, to use in graphs where edges have no metadata.
"""
function Graphs.add_edge!(
    g::AbstractDataGraph{T,VL,VD,ED}, label_s::VL, label_d::VL
) where {T,VL,VD,ED<:Nothing}
    return add_edge!(g, label_s, label_d, nothing)
end


text:
for N in [Float64, Rational{Int}, Float32]
    # random ball
    rand(BallInf)

    # 1D BallInf
    b = BallInf(N[0], N(1))
    # Test Dimension
    @test dim(b) == 1
    # Test Support Vector
    d = N[1]
    @test σ(d, b) == N[1]
    d = N[-1]
    @test σ(d, b) == N[-1]

    # 2D BallInf
    b = BallInf(N[0, 0], N(1))
    # Test Dimension
    @test dim(b) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, b) == N[1, 1]
    d = N[-1, 1]
    @test σ(d, b) == N[-1, 1]
    d = N[-1, -1]
    @test σ(d, b) == N[-1, -1]
    d = N[1, -1]
    @test σ(d, b) == N[1, -1]

    # 2D BallInf not 0-centered
    b = BallInf(N[1, 2], N(1))
    # Test Dimension
    @test dim(b) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, b) == N[2, 3]
    d = N[-1, 1]
    @test σ(d, b) == N[0, 3]
    d = N[-1, -1]
    @test σ(d, b) == N[0, 1]
    d = N[0, -1]
    @test σ(d, b) == N[2, 1]

    # 2D BallInf radius =/= 1
    b = BallInf(N[0, 0], N(2))
    # Test Dimension
    @test dim(b) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, b) == N[2, 2]
    d = N[-1, 1]
    @test σ(d, b) == N[-2, 2]
    d = N[-1, -1]
    @test σ(d, b) == N[-2, -2]
    d = N[1, -1]
    @test σ(d, b) == N[2, -2]

    # support vector for single entry vector
    svec = σ(SingleEntryVector(2, 3, N(2)), BallInf(zeros(N, 3), N(2)))
    @test svec[1] ∈ Interval(N[-2, 2]) && svec[2] == N(2) && svec[3] ∈ Interval(N[-2, 2])

    # support function
    B = BallInf(N[1, 2], N(1))
    @test ρ(N[1, 1], B) == N(5)
    @test ρ(N[1, 0], B) == N(2)
    @test ρ(N[0, 1], B) == N(3)
    @test ρ(N[-1, -1], B) == N(-1)
    @test ρ(N[-1, 0], B) == N(0)
    @test ρ(N[-1, 1], B) == N(3)

    # boundedness
    @test isbounded(b)

    # isempty
    @test !isempty(b)

    # isuniversal
    answer, w = isuniversal(b, true)
    @test !isuniversal(b) && !answer && w ∉ b

    # membership
    b = BallInf(N[1, 1], N(1))
    @test N[0.5, -0.5] ∉ b
    @test N[0.5, 1.5] ∈ b

    # an_element function
    b = BallInf(N[1, 2], N(3))
    @test an_element(b) ∈ b

    # check that vertices_list for zero radius doesn't repeat vertices
    b = BallInf(N[1, 2], N(0))
    vl = vertices_list(b)
    @test vl == [center(b)]

    # high and low
    b = BallInf(N[1, 2], N(1))
    @test high(b) == N[2, 3]
    @test low(b) == N[0, 1]

    # isflat
    @test !isflat(BallInf(N[1, 1], N(1))) && isflat(BallInf(N[1, 1], N(0)))

    # split
    b = BallInf(N[3, 3], N(1))
    @test split(b, [1, 1]) == [Hyperrectangle(N[3, 3], N[1, 1])]
    @test ispermutation(split(b, [2, 1]),
        [Hyperrectangle(N[2.5, 3], N[0.5, 1]),
         Hyperrectangle(N[3.5, 3], N[0.5, 1])])

    # translation
    b = BallInf(N[1, 2], N(1))
    @test translate(b, N[1, 2]) == BallInf(N[2, 4], N(1))

    # affine map
    M = N[2 1; 0 2]
    B = BallInf(N[0, 0], N(1))
    v = N[1, -1]
    am = affine_map(M, B, v)
    @test ispermutation(vertices_list(am),
                        [N[4, 1], N[0, 1], N[-2, -3], N[2, -3]])

    # volume
    B = BallInf(N[0, 0], N(1))
    @test volume(B) ≈ N(4)
    if N <: AbstractFloat
        B = BallInf(zeros(N, 100), N(1/2 + 1e-5))
        @test volume(B) ≈ N(1.0020019812942185)
    end
end

# tests that only work with Float64
for N in [Float64]
    # concrete Minkowski sum
    b = BallInf(N[1, 2], N(1))
    p = minkowski_sum(b, N[2 0; 0 1] * b)
    @test p isa HPolytope{N} && ispermutation(constraints_list(p),
        [HalfSpace(N[0, -1], N(-2)), HalfSpace(N[0, 1], N(6)),
         HalfSpace(N[-1, 0], N(0)), HalfSpace(N[1, 0], N(6))])
end


text:
# tên file bản đồ
maps = [
    "012-1x2.csv",
    "013-1x2.csv",
    "014-1x2.csv",
    "015-1x2.csv",
    "016-1x2.csv",
    "017-1x2.csv",
    "018-1x2.csv",
    "019-1x2.csv",
    "020-1x2.csv",
    "12-3x3.csv",
    "13-3x3.csv",
    "14-3x3.csv",
    "15-3x3.csv",
    "16-3x3.csv",
    "17-3x3.csv",
    "18-3x3.csv",
    "19-3x3.csv",
    "20-3x3.csv",
]
maps = joinpath.("src", "maps", maps)

# positions
# - :corner
# - :random_c1
# - :random_c2
# - :border
positions = [:corner, :border]

# nb_bph_init
# - bất kì số tự nhiên
nb_bph_inits = [200, 20]

# pr_killed
# - bất kì số thực
pr_killeds = [0.075, 0.09, 0.12, 0.15]

# Đầu vào là tích descartes của những mảng trên
params1 = map(Iterators.product(maps, nb_bph_inits, pr_killeds, positions)) do (m, n, pr, p)
    (envmap=m, nb_bph_init=n, pr_killed=pr, init_position=p)
end
params1 = params1[:]

const no_flower_map = joinpath("src", "maps", "no-flower.csv")
const params2 = [
    (envmap=no_flower_map, pr_killed=0.15, nb_bph_init=20, init_position=:corner),
    (envmap=no_flower_map, pr_killed=0.15, nb_bph_init=200, init_position=:corner),
    (envmap=no_flower_map, pr_killed=0.15, nb_bph_init=20, init_position=:border),
    (envmap=no_flower_map, pr_killed=0.15, nb_bph_init=200, init_position=:border),
]

const params3 = let maps = ["019-1x2.csv", "19-3x3.csv"]
    local p = map(Iterators.product(maps, nb_bph_inits, positions)) do (envmap, n, p)
        (envmap=joinpath("src", "maps", envmap), pr_killed=0.0, nb_bph_init=n, init_position=p)
    end
    collect(p)[:]
end

params = [params1; params2; params3]

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
# If you want to run this matrix-benchmarks, you should install
# a dependency `Benchmarkplotting` which is not registered on JuliaRegistries.
#    pkg> add https://github.com/thautwarm/Benchmarkplotting.jl
using MLStyle
using ArgParse
include("matrix-benchmark/sampler.jl")
export ArbitrarySampler
versus_items = ("misc", "tuple", "array", "structfields", "vs-match", "datatype")

function parse_cmd()
    s = ArgParseSettings()
    @add_arg_table s begin
        "versus"
        help = join(", ", versus_items)
        required = true
        action => :store_arg
        nargs = '+'
    end
    parse_args(ARGS, s)
end

check_versus(x) = x in versus_items

function benchmark(x)
    filename = "matrix-benchmark/bench-$x.jl"
    open(filename) do f
        include_string(Main, read(f, String), filename)
    end
end

action = @λ begin
    [] -> nothing
    [GuardBy(check_versus) && hd, tl...] -> begin
        benchmark(hd)
        action(tl)
    end
    [hd, tl...] -> begin
        @warn "Unknown versus item: $hd"
        action(hd)
    end
end

function main()
    parsed_args = parse_cmd()
    Dict([a => b for (a, b) in parsed_args])["versus"] |> action
end

main()


text:
include("regex_parse.jl")
include("regex_transform.jl")
include("regex_build.jl")

function regex_rules(regex::AbstractString, datatype::DataType=String, rulenameprefix="")
	ast = parse_regex(regex, datatype)
	transform_regex_ast(ast)
	transform_ast(ast) # this standard transform (which analysis reachability) isn't really needed for a regex, but included for consistency with other translators
	build_regex_rules(ast, rulenameprefix)
end

function regex_generator(io::IO, genname::Symbol, regex::AbstractString, datatype::DataType=String)
	rules = regex_rules(regex, datatype)
	description = "$(datatype) satisfying regular expression " * escape_string(regex)
	output_generator(io, genname, description, rules)
end

regex_generator(genname::Symbol, regex::AbstractString, datatype::DataType=String) = include_generator(genname, regex_generator, regex, datatype)


text:
module SmartGameFormat

using DataStructures
using Crayons

export

    SGFNode,
    SGFGameTree,

    load_sgf,
    save_sgf,
    parse_sgf,
    print_sgf

include("types.jl")
include("lexer.jl")
include("parser.jl")
include("prettyprint.jl")
include("interface.jl")

end # module


text:
using Pkg; Pkg.activate("../")
using WebAssemblyText
using Test

@testset "WebAssemblyText.jl" begin
    a=1
    b=1
    @test a == b
end


text:
for surrogate_type in surrogate_methods
    surrtype = Symbol("$(surrogate_type)")

    ############################
    # User provided time series
    ############################
    fname1 = Symbol("$(surrtype)_plot")
    @eval begin
        """
                `$($surrtype)_plot(ts)`

            Visualize a `$($surrtype)` surrogate.
        """
        function $(fname1)(ts)
            surrplot(ts, $(surrtype)(ts))
        end
    end

    fname2 = Symbol("$(surrtype)_anim")
    @eval begin
        """
                `$($surrtype)_anim(ts; [n_iters::Int])`

            Create `n_iters` `$($surrtype)` surrogates for `ts` and animate them. Returns a `Plots.Animation` instance.
        """
        function $(fname2)(ts; n_iters = 15)
            anim = Plots.@animate for i = 1:n_iters
                surrplot(ts, $(surrtype)(ts))
            end
        end
    end

    fname2b = Symbol("$(surrtype)_gif")
    @eval begin
        """
                `$($surrtype)_gif(ts; [n_iters::Int])`

            Create `n_iters` `$($surrtype)` surrogates for `ts`, animate them and create
            a gif from the animation.
        """
        function $(fname2b)(ts; fps = 2, n_iters = 15)
            anim = Plots.@animate for i = 1:n_iters
                surrplot(ts, $(surrtype)(ts))
            end
            gif(anim, fps = fps)
        end
    end

    eval(Expr(:export, fname1))
    eval(Expr(:export, fname2))
    eval(Expr(:export, fname2b))

    ################################
    # Synthetic time series examples
    ################################

    for process in processes
        proc = Symbol("$(process)")

        # Surrogate panel plots for combinations of processes and surrogate types
        fname3 = Symbol("$(surrtype)_$(process)_plot")
        @eval begin
            """
                    `$($fname3)(; [fps, n_iters, n_steps, new_realization_every_iter::Bool])`

                Create a `$($surrtype)` surrogate for a realization of a `$($proc)` process and visualize it.
            """
            function $(fname3)(;n_iters = 15,
                                n_steps = 500,
                                new_realization_every_iter::Bool = false)
                # generate time series
                ts = $(proc)(n_steps = n_steps)
                surrplot(ts, $(surrtype)(ts))
            end
        end

        # Animate surrogate generation for different surrogate types and processes
        fname4 = Symbol("$(surrtype)_$(process)_anim")
        @eval begin
            """
                    `$($fname4)(; [fps, n_iters, n_steps, new_realization_every_iter::Bool])`

                Create `n_iters` `$($surrtype)` surrogates for realizations of a `$($proc)` process and animate them. Returns a `Plots.Animation` instance.

                If `new_realization_every_iter = true`, then a fresh `$($proc)` realization is used at each iteration. Otherwise, the same time series, generating only the surrogates anew at each iteration.
            """
            function $(fname4)(;n_iters = 15,
                                n_steps = 500,
                                new_realization_every_iter::Bool = false)
                # generate time series
                ts = $(proc)(n_steps = n_steps)

                # animate
                anim = Plots.@animate for i = 1:n_iters
                    if new_realization_every_iter
                        ts = $(proc)(n_steps = n_steps)
                        surrplot(ts, $(surrtype)(ts))
                    else
                        surrplot(ts, $(surrtype)(ts))
                    end
                end
            end
        end

        # Animate surrogate generation for different surrogate types and processes
        fname4b = Symbol("$(surrtype)_$(process)_gif")
        @eval begin
            """
                    `$($fname4b)(; [fps, n_iters, n_steps, new_realization_every_iter::Bool])`

                Create `n_iters` `$($surrtype)` surrogates for realizations of a `$($proc)`
                process, animate them and create a gif from the animation.

                If `new_realization_every_iter = true`, then a fresh `$($proc)` realization is used at each iteration. Otherwise, the same time series, generating only the surrogates anew at each iteration.
            """
            function $(fname4b)(;fps = 2, n_iters = 15,
                                n_steps = 500,
                                new_realization_every_iter::Bool = false)
                # generate time series
                ts = $(proc)(n_steps = n_steps)

                # animate
                anim = Plots.@animate for i = 1:n_iters
                    if new_realization_every_iter
                        ts = $(proc)(n_steps = n_steps)
                        surrplot(ts, $(surrtype)(ts))
                    else
                        surrplot(ts, $(surrtype)(ts))
                    end
                end

                gif(anim, fps = fps)
            end
        end

        eval(Expr(:export, fname3))
        eval(Expr(:export, fname4))
        eval(Expr(:export, fname4b))
    end
end


text:
# Autogenerated wrapper script for ROCmOpenCLRuntime_jll for x86_64-linux-musl-cxx11
export libOpenCL, libamdocl

using hsa_rocr_jll
using hsakmt_roct_jll
using ROCmDeviceLibs_jll
using ROCmCompilerSupport_jll
using Libglvnd_jll
using Xorg_libX11_jll
using Xorg_xorgproto_jll
JLLWrappers.@generate_wrapper_header("ROCmOpenCLRuntime")
JLLWrappers.@declare_library_product(libOpenCL, "libOpenCL.so.1")
JLLWrappers.@declare_library_product(libamdocl, "libamdocl64.so")
function __init__()
    JLLWrappers.@generate_init_header(hsa_rocr_jll, hsakmt_roct_jll, ROCmDeviceLibs_jll, ROCmCompilerSupport_jll, Libglvnd_jll, Xorg_libX11_jll, Xorg_xorgproto_jll)
    JLLWrappers.@init_library_product(
        libOpenCL,
        "lib64/libOpenCL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libamdocl,
        "lib/libamdocl64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# REPL tests
using TestHelpers

# Writing ^C to the repl will cause sigint, so let's not die on that
ccall(:jl_exit_on_sigint, Void, (Cint,), 0)
# These are integration tests. If you want to unit test test e.g. completion, or
# exact LineEdit behavior, put them in the appropriate test files.
# Furthermore since we are emulating an entire terminal, there may be control characters
# in the mix. If verification needs to be done, keep it to the bare minimum. Basically
# this should make sure nothing crashes without depending on how exactly the control
# characters are being used.
begin
# Use pipes so we can easily do blocking reads
stdin_read,stdin_write = (Base.Pipe(C_NULL), Base.Pipe(C_NULL))
stdout_read,stdout_write = (Base.Pipe(C_NULL), Base.Pipe(C_NULL))
stderr_read,stderr_write = (Base.Pipe(C_NULL), Base.Pipe(C_NULL))
Base.link_pipe(stdin_read,true,stdin_write,true)
Base.link_pipe(stdout_read,true,stdout_write,true)
Base.link_pipe(stderr_read,true,stderr_write,true)

repl = Base.REPL.LineEditREPL(TestHelpers.FakeTerminal(stdin_read, stdout_write, stderr_write))
# In the future if we want we can add a test that the right object
# gets displayed by intercepting the display
repl.specialdisplay = Base.REPL.REPLDisplay(repl)
repl.no_history_file = true

repltask = @async begin
    Base.REPL.run_repl(repl)
end

sendrepl(cmd) = write(stdin_write,"inc || wait(b); r = $cmd; notify(c); r\r")

inc = false
b = Condition()
c = Condition()
sendrepl("\"Hello REPL\"")
inc=true
begin
    notify(b)
    wait(c)
end
# Latex completions
write(stdin_write, "\x32\\alpha\t")
readuntil(stdout_read, "α")
# Bracketed paste in search mode
write(stdin_write, "\e[200~paste here ;)\e[201~")
# Abort search (^C)
write(stdin_write, '\x03')
# Test basic completion in main mode
write(stdin_write, "Base.REP\t")
readuntil(stdout_read, "Base.REPL")
write(stdin_write, '\x03')
write(stdin_write, "\\alpha\t")
readuntil(stdout_read,"α")
write(stdin_write, '\x03')
# Close REPL ^D
write(stdin_write, '\x04')
wait(repltask)
end

ccall(:jl_exit_on_sigint, Void, (Cint,), 1)

let exename=joinpath(JULIA_HOME,(ccall(:jl_is_debugbuild,Cint,())==0?"julia":"julia-debug"))

# Test REPL in dumb mode
@unix_only begin

const O_RDWR = Base.FS.JL_O_RDWR
const O_NOCTTY = Base.FS.JL_O_NOCTTY

fdm = ccall(:posix_openpt,Cint,(Cint,),O_RDWR|O_NOCTTY)
fdm == -1 && error("Failed to open PTY master")
rc = ccall(:grantpt,Cint,(Cint,),fdm)
rc != 0 && error("grantpt failed")
rc = ccall(:unlockpt,Cint,(Cint,),fdm)
rc != 0 && error("unlockpt")

fds = ccall(:open,Cint,(Ptr{Uint8},Cint),ccall(:ptsname,Ptr{Uint8},(Cint,),fdm), O_RDWR|O_NOCTTY)

# slave
slave   = RawFD(fds)
master = Base.TTY(RawFD(fdm); readable = true)

nENV = copy(ENV)
nENV["TERM"] = "dumb"
p = spawn(setenv(`$exename -f --quiet`,nENV),slave,slave,slave)
start_reading(master)
Base.wait_readnb(master,1)
write(master,"1\nquit()\n")

wait(p)

output = readavailable(master)
@test output == "julia> 1\r\nquit()\r\n1\r\n\r\njulia> "

close(master)
ccall(:close,Cint,(Cint,),fds)

end

# Test stream mode
outs, ins, p = readandwrite(`$exename -f --quiet`)
write(ins,"1\nquit()\n")
@test readall(outs) == "1\n"
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using ProgressMeter
using Printf

function main()
    PARAMDAT = "PARAMDAT.json"

    # 膨張率 [%]
    expantion = 10

    # 構造体メモリ
    memori = 20000

    # ----------------------------------------------------------
    # 変数はないはず
    # ----------------------------------------------------------

    xmax,ymax,nodes,nodes_vtk,elements,readQ = read_allgrid()

    out_file,in_file,n_div,s_x,e_x,n_x,s_y,e_y,n_y = input_para(PARAMDAT)

    main_pre(n_x,n_y,s_x,e_x,s_y,e_y)

    # 空間分割
    divnum = 2^n_div

    # 辺の長さ
    lx = abs(e_x - s_x)
    ly = abs(e_y - s_y)
    lx_exp = lx *expantion/200
    ly_exp = ly *expantion/200

    # モートン単位長
    unitx = abs(e_x - s_x) / divnum *(1 + expantion/100)
    unity = abs(e_y - s_y) / divnum *(1 + expantion/100)

    # 空間分割中心算出
    dx = (e_x - s_x) / n_x
    dy = (e_y - s_y) / n_y

    new_cell_center = zeros(n_x,n_y,2)
    new_morton_num = zeros(n_x,n_y,2)
    for i in 1:n_x
        for j in 1:n_y
            # 左上の点のモートン空間番号
            ulx = floor((s_x + dx * (i-1) -s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
            uly = floor((s_y + dy *  j    -s_y + ly_exp) / unity)

            # 右下の点のモートン空間番号
            lrx = floor((s_x + dx *  i    -s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
            lry = floor((s_y + dy * (j-1) -s_y + ly_exp) / unity)

            # 0<=s, 0<=t
            #new_morton_num[i,j,1],new_morton_num[i,j,2] = cal_morton(ulx,uly,lrx,lry)
            s,t = cal_morton(ulx,uly,lrx,lry,n_div)
            new_morton_num[i,j,1] = s
            new_morton_num[i,j,2] = t

            new_cell_center[i,j,1] = s_x +dx/2+dx*(i-1)
            new_cell_center[i,j,2] = s_y +dy/2+dy*(j-1)
        end
    end
    println(" fin kuukan bunkatu\n")

    # 登録管理構造体
    s = Int64(0)
    for i in 1:(n_div+1)
        s = s + 4^(i-1)
    end
    morton_tree = zeros(Int64,s,memori)
    morton_tree_num = ones(Int64,s)

    # セルを囲む四角形
    # これでモートン空間番号を算出
    cellnum = size(elements)[1]
    cell_center = zeros(cellnum,2)
    #morton_num = zeros(cellnum,2)

    for i in 1:cellnum
        x1 = nodes_vtk[elements[i,1],1]
        y1 = nodes_vtk[elements[i,1],2]
        x2 = nodes_vtk[elements[i,2],1]
        y2 = nodes_vtk[elements[i,2],2]
        x3 = nodes_vtk[elements[i,3],1]
        y3 = nodes_vtk[elements[i,3],2]
        x4 = nodes_vtk[elements[i,4],1]
        y4 = nodes_vtk[elements[i,4],2]

        temp = [x1,x2,x3,x4]
        maxx = maximum(temp)
        minx = minimum(temp)
        temp = [y1,y2,y3,y4]
        maxy = maximum(temp)
        miny = minimum(temp)

        # 左上の点のモートン空間番号
        ulx = floor((minx-s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
        uly = floor((maxy-s_y + ly_exp) / unity)

        # 右下の点のモートン空間番号
        lrx = floor((maxx-s_x + lx_exp) / unitx) # 原点に合わせて番号を付けるため（-s_x + lx_exp）
        lry = floor((miny-s_y + ly_exp) / unity)

        # =-------------異常値をはじく事-----------
        if (ulx < 0 || uly <0) || (lrx < 0 || lry <0)
            continue
        elseif (ulx < 0 || uly <0) || (lrx < 0 || lry <0)
            continue
        else
            #morton_num[i,1],morton_num[i,2] = cal_morton(ulx,uly,lrx,lry)
            # 0<=s, 0<=t
            s,t = cal_morton(ulx,uly,lrx,lry,n_div)

            num = Int64(t + (4^s-1)/3   +1)     # 線形空間番号

            s = morton_tree_num[num]
            morton_tree[num,s] = i
            morton_tree_num[num] += 1
        end

        cell_center[i,1] = 0.25*(x1+x2+x3+x4)
        cell_center[i,2] = 0.25*(y1+y2+y3+y4)
    end
    println(" fin kuukan bunkatu2 \n")

    # 線形探索用
    s,search0,search1,search2,search3,search4 = liner_morton(n_div)

    # 分割したセル中心と近い点を三点計算
    close_point_num = zeros(n_x,n_y,3)
    distance = zeros(n_x,n_y,3)
    distance[:,:,:] .= 100

    prog = Progress(n_x,1)
    for i in 1:n_x
        next!(prog)
        for j in 1:n_y
            if new_morton_num[i,j,1] ==0      # 所属空間
                temp = copy(search0)
                loop = Int(s[1] * 4^0)
            elseif new_morton_num[i,j,1] ==1      # 所属空間
                temp = copy(search1)
                loop = Int(s[2] * 4^1)
            elseif new_morton_num[i,j,1] ==2      # 所属空間
                temp = copy(search2)
                loop = Int(s[3] * 4^2)
            elseif new_morton_num[i,j,1] ==3      # 所属空間
                temp = copy(search3)
                loop = Int(s[4] * 4^3)
            elseif new_morton_num[i,j,1] ==4      # 所属空間
                temp = copy(search4)
                loop = Int(s[5] * 4^4)
            else
                " stop program !!"
                throw(UndefVarError(:x))
            end

            # 空間事に探査
            for k in 1:loop
                for l in 1:memori

                    if morton_tree[temp[k],l] == 0
                        break
                    end
                    bangou = morton_tree[temp[k],l]

                    x1 = new_cell_center[i,j,1]
                    y1 = new_cell_center[i,j,2]
                    x2 = cell_center[bangou,1]
                    y2 = cell_center[bangou,2]
                    temp_distance = ((x1-x2)^2 + (y1-y2)^2)^0.5

                    for ii in 1:3
                        if temp_distance < distance[i,j,ii]
                            for m in 1:(4-ii)-1
                                distance[i,j,4-m] = distance[i,j,3-m]
                                close_point_num[i,j,4-m] = close_point_num[i,j,3-m]
                            end
                            distance[i,j,ii] = temp_distance
                            close_point_num[i,j,ii] = morton_tree[temp[k],l]
                            break
                        end
                    end
                end
            end
        end
    end

    println(" fin kuukan bunkatu3 \n")

    println("\n ---------------------------------- \n")
    println(" start interpolation ")
    println("\n ---------------------------------- \n")

    # 補間 loop
    Qcell = zeros(n_x,n_y,5)
    for i in 1:n_x
        for j in 1:n_y
            # cell中心ベクトル p
            p  = [new_cell_center[i,j,1], new_cell_center[i,j,2]]

            # 四点ベクトル p
            x0 = cell_center[Int64(close_point_num[i,j,1]),1]
            y0 = cell_center[Int64(close_point_num[i,j,1]),2]
            x1 = cell_center[Int64(close_point_num[i,j,2]),1]
            y1 = cell_center[Int64(close_point_num[i,j,2]),2]
            x2 = cell_center[Int64(close_point_num[i,j,3]),1]
            y2 = cell_center[Int64(close_point_num[i,j,3]),2]

            p0 = [x0, y0]
            p1 = [x1, y1]
            p2 = [x2, y2]

            # 局所座標ベクトル e
            eu = [p1[1] - p0[1], p1[2] - p0[2]]
            ev = [p2[1] - p0[1], p2[2] - p0[2]]
            inv_A = inverse_matrix(eu,ev)

            # 補間距離ベクトル u
            up = zeros(2)
            for l in 1:2
                up[l] = inv_A[l,1]*(p[1]-p0[1]) + inv_A[l,2]*(p[2]-p0[2])
            end
            for k in 1:5
                Qcell[i,j,k] = readQ[Int64(close_point_num[i,j,1]),k]
                            + up[1]*readQ[Int64(close_point_num[i,j,2]),k]
                            + up[2]*readQ[Int64(close_point_num[i,j,3]),k]
                #Qcell[i,j] = Qcell[i,j]*avoga
            end
        end
    end

    println("\n ---------------------------------- \n")
    println(" start atari ")
    println("\n ---------------------------------- \n")

    # 分割した空間がセルと当たるか否かの判定
    atari = zeros(Int64,n_x,n_y)
    for i in 1:n_x
        for j in 1:n_y
            for k in 1:3
                # ---------- 過大評価 ---------------
                cnum = Int64(close_point_num[i,j,k])

                x1 = nodes_vtk[elements[cnum,1],1]
                y1 = nodes_vtk[elements[cnum,1],2]
                x2 = nodes_vtk[elements[cnum,2],1]
                y2 = nodes_vtk[elements[cnum,2],2]
                x3 = nodes_vtk[elements[cnum,3],1]
                y3 = nodes_vtk[elements[cnum,3],2]
                x4 = nodes_vtk[elements[cnum,4],1]
                y4 = nodes_vtk[elements[cnum,4],2]

                temp_x = [x1,x2,x3,x4]
                maxx = maximum(temp_x)
                minx = minimum(temp_x)
                temp_y = [y1,y2,y3,y4]
                maxy = maximum(temp_y)
                miny = minimum(temp_y)

                x = new_cell_center[i,j,1]
                y = new_cell_center[i,j,2]

                if minx < x && x < maxx
                    if miny < y && y < maxy
                        atari[i,j] += 1
                    end
                end
            end
        end
    end

    # output

    fff = "grid_morton/atari"
    open(fff,"w") do f
        write(f,"result:atari\n")
        for i in 1:n_x
            for j in 1:n_y
                a1 = string(atari[i,j])
                write(f, a1*"\n")
            end
        end
    end
    println("\nwrite "*fff)

    fff = "result_morton/"* out_file
    open(fff,"w") do f
        write(f,"result:rho[kg/m^3], u[m/s], v[m/s], p[Pa], T[K]\n")
        for i in 1:n_x
            for j in 1:n_y

                a1 = @sprintf("%8.8f", Qcell[i,j,1])
                a2 = @sprintf("%8.8f", Qcell[i,j,2])
                a3 = @sprintf("%8.8f", Qcell[i,j,3])
                a4 = @sprintf("%8.8f", Qcell[i,j,4])
                a5 = @sprintf("%8.8f", Qcell[i,j,5])

                #=
                a1 = string(Int64(close_point_num[i,j,1]))
                a2 = string(Int64(close_point_num[i,j,2]))
                a3 = string(Int64(close_point_num[i,j,3]))
                =#

                write(f, a1*" "*a2*" "*a3*" "*a4* " "*a5*"\n")
            end
        end
    end
    println("\nwrite "*fff)

end

# ----------------
main()
# ----------------

max_stars_repo_name: hide-dog/Qtree

max_issues_repo_name: hide-dog/Qtree

max_forks_repo_name: hide-dog/Qtree


text:
module BraidChains

import PeaceCypher: sign, hash
using PeaceCypher

using SynchronicBallot
import SynchronicBallot: Mix, GateKeeper

using Tables
using Sockets
using Pkg.TOML

# For keychain I will likely need a method `newsigner`

abstract type Agent end
abstract type Action end
abstract type Braid end

include("core.jl")
include("sbraids.jl")
include("sbraider.jl")
include("keychains.jl")

# I could capitalize every letter to sepperate Agents from everyone else
# export GUARDIAN, REGISTRATOR, MEMBER, BRAIDER

export Guardian, Registrator, Member, Braider, Agent
export Transaction, pseudonyms, BraidChain, Contract, state

export BraiderConfig, Mix, GateKeeper, BraidOfficer, braid

export TransactionLog, TransactionHole, Trigger, TransactionVector

export KeyChain, braid!, update!

end # module


text:
# Types

export
    Interaction,
    GeneralInteraction,
    SpecificInteraction,
    Simulator,
    Thermostat,
    NeighbourFinder,
    Logger,
    Atom,
    Simulation

const DefaultFloat = Float64

"An interaction between atoms that contributes to forces on the atoms."
abstract type Interaction end

"""
A general interaction that will apply to all atom pairs.
Custom general interactions should sub-type this type.
"""
abstract type GeneralInteraction <: Interaction end

"""
A specific interaction between sets of specific atoms, e.g. a bond angle.
Custom specific interactions should sub-type this type.
"""
abstract type SpecificInteraction <: Interaction end

"""
A type of simulation to run, e.g. leap-frog integration or energy minimisation.
Custom simulators should sub-type this type.
"""
abstract type Simulator end

"""
A way to keep the temperature of a simulation constant.
Custom thermostats should sub-type this type.
"""
abstract type Thermostat end

"""
A way to find near atoms to save on simulation time.
Custom neighbour finders should sub-type this type.
"""
abstract type NeighbourFinder end

"""
A way to record a property, e.g. the temperature, throughout a simulation.
Custom loggers should sub-type this type.
"""
abstract type Logger end

"""
    Atom(; <keyword arguments>)

An atom and its associated information.
Properties unused in the simulation or in analysis can be left with their
default values.

# Arguments
- `attype::AbstractString=""`: the type of the atom.
- `name::AbstractString=""`: the name of the atom.
- `resnum::Integer=0`: the residue number if the atom is part of a polymer.
- `resname::AbstractString=""`: the residue name if the atom is part of a
    polymer.
- `charge::T=0.0`: the charge of the atom, used for electrostatic interactions.
- `mass::T=0.0`: the mass of the atom.
- `σ::T=0.0`: the Lennard-Jones finite distance at which the inter-particle
    potential is zero.
- `ϵ::T=0.0`: the Lennard-Jones depth of the potential well.
"""
struct Atom{T}
    attype::String
    name::String
    resnum::Int
    resname::String
    charge::T
    mass::T
    σ::T
    ϵ::T
end

# We define constructors rather than using Base.@kwdef as it makes conversion
#   more convenient with the parametric type
function Atom(;
                attype="",
                name="",
                resnum=0,
                resname="",
                charge=0.0,
                mass=0.0,
                σ=0.0,
                ϵ=0.0)
    return Atom{typeof(charge)}(attype, name, resnum, resname, charge, mass, σ, ϵ)
end

function Base.show(io::IO, a::Atom)
    print(io, "Atom{", typeof(a.charge), "} with name \"", a.name,
                "\", type \"", a.attype, "\", charge=", a.charge, ", mass=",
                a.mass, ", σ=", a.σ, ", ϵ=", a.ϵ)
end

"""
    Simulation(; <keyword arguments>)

The data needed to define and run a molecular simulation.
Properties unused in the simulation or in analysis can be left with their
default values.

# Arguments
- `simulator::Simulator`: the type of simulation to run.
- `atoms::Vector{A}`: the atoms, or atom equivalents, in the simulation. Can be
    of any type.
- `specific_inter_lists::SI=()`: the specific interactions in the simulation,
    i.e. interactions between specific atoms such as bonds or angles. Typically
    a `Tuple`.
- `general_inters::GI=()`: the general interactions in the simulation, i.e.
    interactions between all or most atoms such as electrostatics. Typically a
    `Tuple`.
- `coords::C`: the coordinates of the atoms in the simulation. Typically a
    `Vector` of `SVector`s of any dimension and type `T`, where `T` is an
    `AbstractFloat` type.
- `velocities::C=zero(coords)`: the velocities of the atoms in the simulation,
    which should be the same type as the coordinates. The meaning of the
    velocities depends on the simulator used, e.g. for the `VelocityFreeVerlet`
    simulator they represent the previous step coordinates for the first step.
- `temperature::T=0.0`: the temperature of the simulation.
- `box_size::T`: the size of the cube in which the simulation takes place.
- `neighbour_finder::NeighbourFinder=NoNeighbourFinder()`: the neighbour finder
    used to find close atoms and save on computation.
- `thermostat::Thermostat=NoThermostat()`: the thermostat which applies during
    the simulation.
- `loggers::Dict{String, <:Logger}=Dict()`: the loggers that record properties
    of interest during the simulation.
- `timestep::T=0.0`: the timestep of the simulation.
- `n_steps::Integer=0`: the number of steps in the simulation.
- `n_steps_made::Vector{Int}=[]`: the number of steps already made during the
    simulation. This is a `Vector` to allow the `struct` to be immutable.
"""
struct Simulation{T, A, C, GI, SI}
    simulator::Simulator
    atoms::Vector{A}
    specific_inter_lists::SI
    general_inters::GI
    coords::C
    velocities::C
    temperature::T
    box_size::T
    neighbour_finder::NeighbourFinder
    thermostat::Thermostat
    loggers::Dict{String, <:Logger}
    timestep::T
    n_steps::Int
    n_steps_made::Vector{Int}
end

function Simulation(;
                    simulator,
                    atoms,
                    specific_inter_lists=(),
                    general_inters=(),
                    coords,
                    velocities=zero(coords),
                    temperature=0.0,
                    box_size,
                    neighbour_finder=NoNeighbourFinder(),
                    thermostat=NoThermostat(),
                    loggers=Dict{String, Logger}(),
                    timestep=0.0,
                    n_steps=0,
                    n_steps_made=[0])
    T = typeof(timestep)
    A = eltype(atoms)
    C = typeof(coords)
    GI = typeof(general_inters)
    SI = typeof(specific_inter_lists)
    return Simulation{T, A, C, GI, SI}(
                simulator, atoms, specific_inter_lists, general_inters, coords,
                velocities, temperature, box_size, neighbour_finder, thermostat,
                loggers, timestep, n_steps, n_steps_made)
end

function Base.show(io::IO, s::Simulation)
    print(io, "Simulation with ", length(s.coords), " atoms, ",
                typeof(s.simulator), " simulator, ", s.timestep, " timestep, ",
                s.n_steps, " steps, ", first(s.n_steps_made), " steps made")
end


text:
using AutoHashEquals
using Base.Test
using Compat
VERSION >= v"0.4" && using Base.Markdown: plain

function sausage(x)
    buf = IOBuffer()
    serialize(buf, x)
    seekstart(buf)
    deserialize(buf)
end

@auto_hash_equals type A
    a::Int
    b
end
@test typeof(A(1,2)) == A
@test hash(A(1,2)) == hash(2,hash(1,hash(:A)))
@test A(1,2) == A(1,2)
@test A(1,2) == sausage(A(1,2))
@test A(1,2) != A(1,3)
@test A(1,2) != A(3,2)

abstract B
@auto_hash_equals immutable C<:B x::Int end
@auto_hash_equals immutable D<:B x::Int end
@test isa(C(1), B)
@test isa(D(1), B)
@test C(1) != D(1)
@test hash(C(1)) != hash(D(1))
@test C(1) == C(1)
@test C(1) == sausage(C(1))
@test hash(C(1)) == hash(C(1))

if VERSION < v"0.4-" typealias Void Nothing end
@compat abstract E{N<:Union{Void,Int}}
@auto_hash_equals type F{N}<:E{N} e::N end
@auto_hash_equals type G{N}<:E{N}
    e::N
end
G() = G{Void}(nothing)
@test hash(F(1)) == hash(1, hash(:F))
@test hash(F(1)) != hash(F(2))
@test F(1) == F(1)
@test F(1) == sausage(F(1))
@test F(1) != F(2)
@test hash(G()) == hash(nothing, hash(:G))
@test G() == G()

macro dummy(x)
    x
end
@test @dummy(1) == 1

@auto_hash_equals type H
    @dummy h
    i
    @dummy j::Int
end
@test H(1,2,3) != H(2,1,3)
@test H(1,2,3) == H(1,2,3)
@test H(1,2,3) == sausage(H(1,2,3))
@test hash(H(1,2,3)) == hash(H(1,2,3))
@test hash(H(1,2,3)) != hash(H(2,1,3))

@auto_hash_equals immutable I{A,B}
    a::A
    b::B
end
@test I{AbstractString,Int}("a", 1) == I{AbstractString,Int}("a", 1)
@test I{AbstractString,Int}("a", 1) == sausage(I{AbstractString,Int}("a", 1))

macro cond(test, block)
    if eval(test)
        block
    end
end

@cond VERSION >= v"0.4" && VERSION < v"0.5-" begin
    println("0.4 specific doc test")
    @doc """this is my data type""" ->
    @auto_hash_equals type MyType
        field::Int
    end
    @test plain(@doc MyType) == "this is my data type\n"
end

@cond VERSION >= v"0.5-pre" begin
    println("0.5 specific doc test")
    """this is my data type"""
    @auto_hash_equals type MyType
        field::Int
    end
    @test_broken plain(@doc MyType) == "this is my data type\n"
end

println("ok")


text:
__precompile__()

module QPControl

# Contact-related
export
    ContactPoint,
    disable!

# Task-related
export
    AbstractMotionTask,
    SpatialAccelerationTask,
    AngularAccelerationTask,
    LinearAccelerationTask,
    PointAccelerationTask,
    JointAccelerationTask,
    MomentumRateTask,
    LinearMomentumRateTask,
    setdesired!

# Low-level
export
    MomentumBasedController,
    SE3PDController,
    addtask!,
    addcontact!,
    regularize!,
    centroidal_frame

# High-level
export
    StandingController

using LinearAlgebra
using RigidBodyDynamics
using RigidBodyDynamics.Graphs
using RigidBodyDynamics.Contact
using RigidBodyDynamics.PDControl
using Parametron
using StaticArrays
using Rotations

import MathOptInterface

const MOI = MathOptInterface
const RBD = RigidBodyDynamics

include("contacts.jl")
include("tasks.jl")
include("exceptions.jl")
include(joinpath("lowlevel", "momentum.jl"))
include(joinpath("lowlevel", "se3pdcontroller.jl"))
include(joinpath("highlevel", "standing.jl"))
include(joinpath("trajectories", "trajectories.jl"))

end # module


text:
module Mercury

import ZMQ
import Sockets
import ProtoBuf
import Logging
import StaticArrays

# Import correct library name for specific system
libhg_library_filename = ""
if Sys.islinux()
    libhg_library_filename = joinpath(@__DIR__, "..", "deps", "build", "libhg.so")
elseif Sys.isapple()
    libhg_library_filename = joinpath(@__DIR__, "..", "deps", "build", "libhg.dylib")
end
const libhg = libhg_library_filename

const MercuryMessage = Union{ProtoBuf.ProtoType,AbstractVector{UInt8}}

include("utils.jl")
include("serial_relay.jl")
include("rate_limiter.jl")

include("publishers/abstract_publisher.jl")
include("publishers/zmq_publisher.jl")

include("subscribers/abstract_subscriber.jl")
include("subscribers/zmq_subscriber.jl")

include("node.jl")

end # module


text:
using BSON
using Flux

# Export the Neural Network given their weights
function export_NN(weights_str, NN_str, uw_NN, vw_NN, wT_NN)
    PATH = pwd()

    _, re_uw = Flux.destructure(uw_NN)
    _, re_vw = Flux.destructure(vw_NN)
    _, re_wT = Flux.destructure(wT_NN)

    uw_weights = BSON.load(joinpath(PATH, "Output", "$(weights_str[1]).bson"))[:weights]
    vw_weights = BSON.load(joinpath(PATH, "Output", "$(weights_str[2]).bson"))[:weights]
    wT_weights = BSON.load(joinpath(PATH, "Output", "$(weights_str[3]).bson"))[:weights]

    uw_NN = re_uw(uw_weights)
    vw_NN = re_vw(vw_weights)
    wT_NN = re_wT(wT_weights)

    bson(joinpath(PATH, "NDEs", "$(NN_str[1]).bson"), Dict(:neural_network => uw_NN))
    bson(joinpath(PATH, "NDEs", "$(NN_str[2]).bson"), Dict(:neural_network => vw_NN))
    bson(joinpath(PATH, "NDEs", "$(NN_str[3]).bson"), Dict(:neural_network => wT_NN))
end

# uw_NN = BSON.load(joinpath(PATH, "NDEs", "uw_NDE_SWNH_100.bson"))[:neural_network]
# vw_NN = BSON.load(joinpath(PATH, "NDEs", "vw_NDE_SWNH_100.bson"))[:neural_network]
# wT_NN = BSON.load(joinpath(PATH, "NDEs", "wT_NDE_SWNH_100.bson"))[:neural_network]

N_inputs = 96
hidden_units = 400
N_outputs = 31
uw_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs))
vw_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs))
wT_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs))

# File name of the weights
weights_str = ["uw_weights_large_NN", "vw_weights_large_NN", "wT_weights_large_NN"]

# Names of the output files
NN_str = ["uw_NN_large", "vw_NN_large", "wT_NN_large"]

export_NN(weights_str, NN_str, uw_NN, vw_NN, wT_NN)


text:
using Random: GLOBAL_RNG

abstract type SemiMartingale{dType, tType, N} end

"""
Type for holding the increments of a Wiener process

### Fields:
* `nd`: dimension of the Wiener process
* `nt`: number of increments in the DataSeries
* `ns`: number of sample paths of the Wiener process
* `Δt`: time increment of the TimeSeries
* `ΔW`: variable storing the increments of the Wiener process over Δt, or the discrete random variable \\hat{I}
* `ΔZ`: variable holding the time integral of the Wiener process \\int_{tk}^{tk+1} (W(t)-W(tk))dt, or the discrete random variable \\tilde{I}

### Parameters:
* `dType`: type of the elements of the increments of the Wiener process
* `tType`: type of the time steps
* `N`:     the number of dimensions of the arrays holding data in ΔW and ΔZ
* `CONV`:  mode of convergence: `:strong` or `:weak`
"""
struct WienerProcess{dType, tType, N, CONV} <: SemiMartingale{dType, tType, N}
    nd::Int
    nt::Int
    ns::Int
    Δt::tType
    ΔW::Array{dType,N}
    ΔZ::Array{dType,N}

    function WienerProcess{dType, tType, N, CONV}(nd, nt, ns, Δt,
                    ΔW = (ns == 1 ? zeros(dType, nd, nt) : zeros(dType, nd, nt, ns)),
                    ΔZ = (ns == 1 ? zeros(dType, nd, nt) : zeros(dType, nd, nt, ns))) where {dType, tType, N, CONV}

        @assert CONV==:strong || CONV==:weak || CONV==:null
        @assert nd > 0
        @assert nt ≥ 1
        @assert ns > 0
        @assert N ∈ (2,3)

        if CONV != :null
            @assert nd == size(ΔW, 1) == size(ΔZ, 1)
            @assert nt == size(ΔW, 2) == size(ΔZ, 2)
            @assert ns == size(ΔW, 3) == size(ΔZ, 3)
        end

        new(nd, nt, ns, Δt, ΔW, ΔZ)
    end

end

function WienerProcess(dType, nd, nt, ns, Δt::tType, conv=DEFAULT_SCONV; rng=GLOBAL_RNG) where {tType <: Number}
    wp = WienerProcess{dType, tType, ns==1 ? 2 : 3, conv}(nd, nt, ns, Δt)
    generate_wienerprocess!(wp, rng)
    return wp
end

function WienerProcess(Δt::tType, dW::Array{dType,1}, dZ::Array{dType,1}, conv=DEFAULT_SCONV) where {dType, tType}
    @assert size(dW,1) == size(dZ,1)

    nd = 1
    nt = size(dW,1)
    ns = 1

    return WienerProcess{dType, tType, 2, conv}(nd, nt, ns, Δt, reshape(dW, (1,nt)), reshape(dZ, (1,nt)))
end

function WienerProcess(Δt::tType, dW::Array{dType,N}, dZ::Array{dType,N}, conv=DEFAULT_SCONV) where {dType, tType, N}
    @assert size(dW) == size(dZ)

    nd = size(dW,1)
    nt = size(dW,2)
    ns = size(dW,3)

    return WienerProcess{dType, tType, N, conv}(nd, nt, ns, Δt, dW, dZ)
end

Base.:(==)(wp1::WienerProcess{DT1, TT1, N1, C1}, wp2::WienerProcess{DT2, TT2, N2, C2}) where {DT1, TT1, N1, C1, DT2, TT2, N2, C2} = (
                                DT1 == DT2
                             && TT1 == TT2
                             && N1  == N2
                             && C1  == C2
                             && wp1.nd == wp2.nd
                             && wp1.nt == wp2.nt
                             && wp1.ns == wp2.ns
                             && wp1.Δt == wp2.Δt
                             && wp1.ΔW == wp2.ΔW
                             && wp1.ΔZ == wp2.ΔZ)

function set_chi_and_eta!(W::WienerProcess{dType, tType, N, :strong}, chi::AbstractArray{dType}, eta::AbstractArray{dType}) where {dType, tType, N}
    @. W.ΔW = chi * √W.Δt
    @. W.ΔZ = W.Δt^(3/2)/2 * (chi+eta/√3)
    return nothing
end

function set_chi_and_eta!(W::WienerProcess{dType, tType, N, :weak}, chi::AbstractArray{dType}, eta::AbstractArray{dType}) where {dType, tType, N}
    W.ΔW .= 0
    W.ΔZ .= sqrt(W.Δt)

    # Generating the random variable \\hat I
    # P(\\hat I=+-sqrt(3*Δt)) = 1/6
    # P(\\hat I=0) = 2/3
    indx        = findall(x->x<1. / 6., chi)
    W.ΔW[indx] .= -sqrt(3*W.Δt)
    indx        = findall(x->x>5. / 6., chi)
    W.ΔW[indx] .= sqrt(3*W.Δt)

    # Generating the random variable \\tilde I
    # P(\\tilde I=+-sqrt(Δt)) = 1/2
    indx        = findall(x->x<0.5, eta)
    W.ΔZ[indx] .= -sqrt(W.Δt)

    return nothing
end

# Generates a new series of increments for the strong Wiener process W
function generate_wienerprocess!(W::WienerProcess{DT, TT, 2, :strong}, rng=GLOBAL_RNG) where {DT,TT}
    chi = randn(rng, DT, W.nd, W.nt)
    eta = randn(rng, DT, W.nd, W.nt)
    set_chi_and_eta!(W, chi, eta)
end

# Generates a new series of increments for the strong Wiener process W
function generate_wienerprocess!(W::WienerProcess{DT, TT, 3, :strong}, rng=GLOBAL_RNG) where {DT,TT}
    chi = randn(rng, DT, W.nd, W.nt, W.ns)
    eta = randn(rng, DT, W.nd, W.nt, W.ns)
    set_chi_and_eta!(W, chi, eta)
end

# Generates a new series of increments for the weak Wiener process W
function generate_wienerprocess!(W::WienerProcess{DT, TT, 2, :weak}, rng=GLOBAL_RNG) where {DT,TT}
    chi = rand(rng, DT, W.nd, W.nt)
    eta = rand(rng, DT, W.nd, W.nt)
    set_chi_and_eta!(W, chi, eta)
end

# Generates a new series of increments for the weak Wiener process W
function generate_wienerprocess!(W::WienerProcess{DT, TT, 3, :weak}, rng=GLOBAL_RNG) where {DT,TT}
    chi = rand(rng, DT, W.nd, W.nt, W.ns)
    eta = rand(rng, DT, W.nd, W.nt, W.ns)
    set_chi_and_eta!(W, chi, eta)
end

Base.ndims(ds::WienerProcess{DT,TT,N}) where {DT,TT,N} = N


text:
using Documenter, LMDB

makedocs(
    modules = [LMDB],
    clean = false,
    format = Documenter.HTML(),
    sitename = "LMDB.jl",
    authors = "Art Wild, Fabian Gans",
    pages = [
        "Home" => "index.md",
        "Manual" => "manual.md",
        "API" => [
            "Index"=>"api/index.md",
        ]
    ]
)

deploydocs(
    repo = "github.com/wildart/LMDB.jl.git",
)


text:
"Derive Effective Concentration with Regression.jl"

include("../src/Regression.jl")
include("../src/Model.jl")
include("../src/Utils.jl")
include("../src/Opt.jl")
import PyPlot
using Printf
using Test
using DifferentialEquations
using Parameters

const plt = PyPlot

## Model and Domain
sample_model = Model.hill(k=0.5, n=4.)
domain = [0.0,1.0]
data = Utils.screen_domain(sample_model, domain)

## Get Thresholds
thr_50 = Regression.measure_EC(data.x, data.y, 0.5)[1]
thr_90 = Regression.measure_EC(data.x, data.y, 0.9)[1]
thr_10 = Regression.measure_EC(data.x, data.y, 0.1)[1]
@test Regression.measure_EC(sample_model, domain,0.5,0.01)[1] ≈ thr_50 # Simulation method and data method

### Get extremes
argmin = Regression.measure_EC(data.x, data.y, 0.)[1]
argmax = Regression.measure_EC(data.x, data.y, 1-1e-10)[1]

### Visualization

fig, ax = plt.subplots(figsize=(10,7))
color = ["black", "green", "red", "blue", "gray"]
ax.plot(data.x, data.y, label="data", color=color[1])
ax.axvline(thr_50, label= string("EC50 = $(@sprintf("%.2f", thr_50))"), color=color[2])
ax.axvline(thr_90, label="EC90 = $(@sprintf("%.2f", thr_90))", color=color[3])
ax.axvline(thr_10, label="EC10 = $(@sprintf("%.2f", thr_10))", color=color[4])

ax.axvline(argmin, label="argMin = $(@sprintf("%.2f", min))", color=color[5])
ax.axvline(argmax, label="argMax = $(@sprintf("%.2f", max))", color=color[5])

ax.legend()
ax.set_title("Estimation of Effective Concentration (EC)")
ax.set_xlabel("Input")
ax.set_ylabel("Output (Monotonous)")
plt.display_figs()

## Positive Feedback

input = LinRange(0.1,1.0,1000)
output = zeros(length(input))

for (i, s) in enumerate(input)
    local mod = Model.positive_feedback(s=s)

    local prob = SteadyStateProblem(mod, mod.x)

    local sol = solve(prob, DynamicSS(Rosenbrock23()))

    output[i] = sol.u[1]
end

thr_50 = Regression.measure_EC(input, output, 0.5)[1]
thr_90 = Regression.measure_EC(input, output, 0.9)[1]
thr_10 = Regression.measure_EC(input, output, 0.1)[1]

### Get extremes
argmin = Regression.measure_EC(input, output, 0.)[1]
argmax = Regression.measure_EC(input, output, 1-1e-10)[1]

fig, ax = plt.subplots(figsize=(10,7))
color = ["black", "green", "red", "blue", "gray"]
ax.plot(input, output, label="data", color=color[1])
ax.axvline(thr_50, label= string("EC50 = $(@sprintf("%.2f", thr_50))"), color=color[2])
ax.axvline(thr_90, label="EC90 = $(@sprintf("%.2f", thr_90))", color=color[3])
ax.axvline(thr_10, label="EC10 = $(@sprintf("%.2f", thr_10))", color=color[4])

ax.axvline(argmin, label="argMin = $(@sprintf("%.2f", min))", color=color[5])
ax.axvline(argmax, label="argMax = $(@sprintf("%.2f", max))", color=color[5])

ax.set_xscale("log")
ax.legend()
ax.set_title("Estimation of Effective Concentration (EC)")
ax.set_xlabel("Input")
ax.set_ylabel("Output (Monotonous)")
plt.display_figs()

##

#=
Utils.screen_domain(Model.positive_feedback, 0.1, :s, [0.0,1.0])

df = Model.positive_feedback()

a = :s
df.a
df.:s

eval(quote df.$a=4 end)
=#


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    IDW(var₁=>param₁, var₂=>param₂, ...)

Inverse distance weighting estimation solver.

## Parameters

* `neighbors` - Number of neighbors (default to all the data)
* `distance`  - A distance defined in Distances.jl (default to `Euclidean()`)
* `power`     - Power of the distances (default to `1`)

### References

Shepard 1968. *A two-dimensional interpolation function for irregularly-spaced data.*
"""
@estimsolver IDW begin
  @param neighbors = nothing
  @param distance = Euclidean()
  @param power = 1
end

function solve(problem::EstimationProblem, solver::IDW)
  # retrieve problem info
  pdata = data(problem)
  pdomain = domain(problem)

  mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))

  # result for each variable
  μs = []; σs = []

  for covars in covariables(problem, solver)
    for var in covars.names
      # get user parameters
      varparams = covars.params[(var,)]

      # determine value type
      V = mactypeof[var]

      # retrieve non-missing data
      locs = findall(!ismissing, pdata[var])
      𝒟 = view(pdata, locs)
      n = nelements(𝒟)

      # determine number of nearest neighbors to use
      k = isnothing(varparams.neighbors) ? n : varparams.neighbors

      # determine distance type
      D = varparams.distance

      # determine power of distances
      p = varparams.power

      @assert n > 0 "estimation requires data"

      @assert k ≤ n "invalid number of neighbors"

      @assert p > 0 "power must be positive"

      # fit search tree
      X = [coordinates(centroid(𝒟, i)) for i in 1:nelements(𝒟)]
      if D isa NearestNeighbors.MinkowskiMetric
        tree = KDTree(X, D)
      else
        tree = BallTree(X, D)
      end

      # lookup non-missing values
      z = 𝒟[var]

      # pre-allocate memory for results
      varμ = Vector{V}(undef, nelements(pdomain))
      varσ = Vector{V}(undef, nelements(pdomain))

      # estimation loop
      for loc in traverse(pdomain, LinearPath())
        x = coordinates(centroid(pdomain, loc))

        is, ds = knn(tree, x, k)
        ws = one(V) ./ ds.^p
        Σw = sum(ws)

        if isinf(Σw) # some distance is zero?
          j = findfirst(iszero, ds)
          μ = z[is[j]]
          σ = zero(V)
        else
          ws /= Σw
          vs  = view(z, is)
          μ = sum(ws[i]*vs[i] for i in eachindex(vs))
          σ = minimum(ds)
        end

        varμ[loc] = μ
        varσ[loc] = σ
      end

      push!(μs, var => varμ)
      push!(σs, Symbol(var,"_distance") => varσ)
    end
  end

  georef((; μs..., σs...), pdomain)
end


text:
###########################################################
# Solution types

###
# Trajectory
#
#   Enhanced Vector type to store temporal information for
#   convenience, and auto-interpolate solution temporally
##

struct Trajectory{T} <: AbstractVector{T}
    v :: AbstractVector{T}
    t :: AbstractRange{Time}
    interpolation
    # Input auto-interpolation
    function Trajectory{T}(t::AbstractRange{Time}, v::AbstractVector{T}) where T
        i = interpolate( (collect(t),), v, Gridded(Constant()) )
        new(v, t, i)
    end
end
@inline Base.@propagate_inbounds Base.getindex(T::Trajectory, i::Int...) = T.v[i...]
(T::Trajectory)(t::Time) = T.interpolation[t]
size(T::Trajectory) = size(T.v)
length(T::Trajectory) = length(T.v)

###
# Ensemble
#
#   Enhanced Array type to store temporal information for
#   convenience, and easily give statistics over the collection
##

struct Ensemble{T} <: AbstractArray{T, 2}
    a :: AbstractArray{T, 2}
    t :: AbstractRange{Time}
    n :: Integer
    function Ensemble{T}(t::AbstractRange{Time}, n::Integer, ::Type{T}) where {T}
        l = length(t)
        a = zeros(T, l, n)
        new(a, t, n)
    end
    function Ensemble{T}(t::AbstractRange{Time}, a::AbstractArray{T, 2}) where {T}
        length(t) == first(size(a)) || error("Incompatible array size")
        n = last(size(a))
        new(a, t, n)
    end
    function Ensemble{T}(tr::Trajectory{T}) where {T}
	let e = Ensemble(tr.t, 1, T)
	    e[1] = tr.v
            e
        end
    end
end
Ensemble(t::R, n::Integer, T::Type) where R<:AbstractRange{Time} = Ensemble{T}(t, n, T)
Ensemble(t::R, a::AbstractVector{T}) where {T,R<:AbstractRange{Time}} = Ensemble{T}(t,reshape(a,length(a),1))
Ensemble(t::R, a::AbstractArray{T,2}) where {T,R<:AbstractRange{Time}} = Ensemble{T}(t,a)
Ensemble(tr::Trajectory{T}) where {T} = Ensemble{T}(tr)
@inline Base.@propagate_inbounds Base.getindex(e::Ensemble, i::Int...) = e.a[i...]
@inline Base.@propagate_inbounds Base.setindex!(e::Ensemble, v::AbstractVector, n::Int) = Base.setindex!(e.a,v,:,n)
(e::Ensemble)(n::Integer) = Trajectory{eltype(e)}(e.t, e.a[:,n])
size(e::Ensemble) = size(e.a)
length(e::Ensemble) = length(e.a)
Base.ndims(e::Ensemble) = ndims(e.a)
#Base.indices(e::Ensemble, d::Integer) = indices(e.a, d)

function map(f::Function, e::Ensemble{T}) where {T}
    max = length(e.t)
    tr = zeros(T, max)
    for i in 1:max
        @inbounds tr[i] = f(e.a[i,:])
    end
    Trajectory{T}(e.t, tr)
end
mean(e::Ensemble) = map(mean, e)
median(e::Ensemble) = map(median, e)
std(e::Ensemble) = map(std, e)

###########################################################
# Simple propagators
##

# Hamiltonian propagation
"""
    ham(dt::Time, H::QOp; ket=false)

Return increment function for Hamiltonian evolution generated
by `H` over a time step `dt`.

Uses an exact matrix exponential, assuming no time-dependence.

### Returns:
  - ket=true  : (t::Time, ψ::QKet) -> u * ψ
  - ket=false : (t::Time, ρ::QOp)  -> u * ρ * u'

"""
function ham(dt::Time, H::QOp; ket=false)
    u::QOp = sparse(exp( -im * dt * Array(H)))
    ut = u'
    if ket
        (t::Time, ψ::QKet) -> u * ψ
    else
        (t::Time, ρ::QOp) -> u * ρ * ut
    end
end
function ham(dt::Time, H::Function; ket=false)
    (t::Time, state) -> ham(dt, H(t), ket=ket)(t, state)
end

# Superoperator Hamiltonian evolution
"""
    sham(dt::Time, H::QOp)

Return increment function using a superoperator for Hamiltonian
evolution generated by `H` over a time step `dt`.

Uses an exact matrix exponential, assuming no time-dependence.

### Returns:
  - (t::Time, ρvec) -> u*ρvec : Evolution superoperator

"""
function sham(dt::Time, H::QOp)
    u::QOp = sparse(exp( -im * dt * Array(H)))
    l = superopl(u)*superopr(u')
    (t::Time, ρvec) -> l * ρvec
end
function sham(dt::Time, H::Function)
    (t::Time, ρvec) -> sham(dt, H(t))(t, ρvec)
end

# Runge-Kutta Hamiltonian evolution
"""
    rinc(t::Time, ρ::QOp, inc::Function, dt::Time)

4th-order Runge-Kutta algorithm given increment `inc`.

### Returns:
  - dρ :: QOp
"""
function rinc(t::Time, ρ::QOp, inc::Function, dt::Time)::QOp
    dρ1::QOp = inc(t, ρ)
    dρ2::QOp = inc(t + dt/2, ρ + dρ1*dt/ 2)
    dρ3::QOp = inc(t + dt/2, ρ + dρ2*dt/ 2)
    dρ4::QOp = inc(t + dt, ρ + dρ3*dt)
    dt*(dρ1 + 2*dρ2 + 2*dρ3 + dρ4)/6
end

"""
    ham_rk4(dt::Time, H::QOp; ket=false)

Return increment function for Hamiltonian evolution generated
by Hamiltonian `H` over a time step `dt`.

Uses a 4th-order Runge-Kutta integration method to construct the state
increment from the first-order differential (master) equation.

### Returns:
  - ket=true  : (t::Time, ψ::QKet) -> ψnew
  - ket=false : (t::Time, ρ::QOp)  -> ρnew
"""
function ham_rk4(dt::Time, H::Function; ket=false)
    inc(t::Time, ρ::QOp)::QOp = - im * comm(H(t),ρ)
    if ket
        inc(t::Time, ψ::QKet)::QKet = - im * H(t) * ψ
    end
    (t::Time, ρ) -> ρ + rinc(t, ρ, inc, dt)
end
function ham_rk4(dt::Time, H::QOp; ket=false)
    h(t) = H
    ham_rk4(dt, h, ket=ket)
end

# Jump-nojump Lindblad propagator
"""
    lind(dt::Time[, H]; clist=QOp[], flist=Function[])

Return increment function over a time step `dt` for Lindblad dissipative
evolution generated by (an optional) Hamiltonian `H` (operator or function),
a list `clist` of constant dissipative operators, and a list `flist` of
time-dependent functions returning dissipative operators.

Uses the "jump no-jump" method to efficiently approximate the exact
Lindblad propagator as a composition of Hamiltonian evolution, jumps,
and no-jump informational backaction. Assumes small dt.
[Physical Review A **92**, 052306 (2015)]

### Returns:
  - (t::Time, ρ(t)::QOp) -> ρ(t+dt)
"""
function lind(dt::Time; clist=QOp[], flist=Function[])
    ns = Function[]
    ds = Function[]
    # Construct operations for constant operators
    if isempty(clist)
        push!(ns, (t, ρ) -> ρ)
    else
        n::QOp = sparse(sqrt(Array(QComp(1.0)LinearAlgebra.I,size(first(clist))) -
                       dt * Array(mapreduce(a -> a' * a, +, clist))))
        push!(ns, (t, ρ) -> n * ρ * n)
        push!(ds, (t, ρ) -> mapreduce(a -> a * ρ * a', +, clist) * dt)
    end
    # Construct operations for time-dependent operators
    if isempty(flist)
        push!(ns, (t, ρ) -> ρ)
    else
        nf(t::Time)::QOp = sparse(sqrt(Array(QComp(1.0)LinearAlgebra.I,size(first(flist)(t))) -
                           dt * Array(mapreduce(a -> a(t)' * a(t), +, flist))))
        push!(ns, (t, ρ) -> nf(t) * ρ * nf(t))
        push!(ds, (t, ρ) -> mapreduce(a -> a(t) * ρ * a(t)', +, flist) * dt)
    end
    push!(ds, (t, ρ) -> last(ns)(t, first(ns)(t, ρ)))
    (t::Time, ρ) -> mapreduce(f -> f(t, ρ), +, ds)
end
function lind(dt::Time, H; clist=QOp[], flist=Function[])
    # Rely on Hamiltonian to specify type of H
    h = ham(dt, H)
    # Apply Hamiltonian first, then the Lindblad increment
    (t::Time, ρ) -> lind(dt, clist=clist, flist=flist)(t, h(t, ρ))
end

# Runge-Kutta Lindblad propagator
"""
    lind_rk4(dt::Time[, H]; clist=QOp[], flist=Function[])

Return increment function over a time step `dt` for Lindblad dissipative
evolution generated by (an optional) Hamiltonian `H` (operator or function),
a list `clist` of constant dissipative operators, and a list `flist` of
time-dependent functions returning dissipative operators.

Uses a 4th-order Runge-Kutta integration method to construct the state
increment from the first-order Lindblad differential (master) equation.

### Returns:
  - (t::Time, ρ(t)::QOp) -> ρ(t) + dρ
"""
lind_rk4(dt::Time, H) = ham_rk4(dt, H)
function lind_rk4(dt::Time; clist=QOp[], flist=Function[])
    inc(t::Time, ρ::QOp)::QOp =
        mapreduce(a -> diss(a)(ρ), +, clist; init=zero(ρ)) +
        mapreduce(a -> diss(a(t))(ρ), +, flist; init=zero(ρ))
    (t::Time, ρ) -> ρ + rinc(t, ρ, inc, dt)
end
function lind_rk4(dt::Time, H::Function; clist=QOp[], flist=Function[])
    inc(t::Time, ρ::QOp)::QOp = - im * comm(H(t),ρ) +
        mapreduce(a -> diss(a)(ρ), +, clist; init=zero(ρ)) +
        mapreduce(a -> diss(a(t))(ρ), +, flist; init=zero(ρ))
    (t::Time, ρ) -> ρ + rinc(t, ρ, inc, dt)
end
function lind_rk4(dt::Time, H::QOp; clist=QOp[], flist=Function[])
    h(t) = H
    lind_rk4(dt, h, clist=clist, flist=flist)
end

# Superoperator Lindblad propagator
"""
    slind(dt::Time, H::QOp; clist=QOp[], flist=Function[])

Return increment function over a time step `dt` for Lindblad dissipative
evolution generated by (an optional) Hamiltonian `H` (operator or function),
a list `clist` of constant dissipative operators, and a list `flist` of
time-dependent functions returning dissipative operators.

Uses direct matrix exponentiation of the total superoperator for
the increment.

### Returns:
  - (t::Time, ρvec) -> u * ρvec : Superoperator for total evolution over dt
"""
slind(dt::Time, H) = sham(dt, H)
function slind(dt::Time, H::QOp; clist=QOp[], flist=Function[])
    h = -im*scomm(H)
    if isempty(flist)
        let l = h + mapreduce(sdiss, +, clist; init=zero(h)),
            u = sparse(exp(dt*Array(l)));
            (t::Time, ρvec) -> u * ρvec
        end
    else
        let l(t::Time) = h + mapreduce(sdiss, +, clist; init=zero(h)) +
                             mapreduce(f -> sdiss(f(t)), +, flist; init=zero(h)),
            u(t::Time) = sparse(exp(dt*Array(l(t))))
            (t::Time, ρvec) -> u(t) * ρvec
        end
    end
end
function slind(dt::Time, H::Function; clist=QOp[], flist=Function[])
    (t::Time, ρvec) -> slind(dt, H(t), clist=clist, flist=flist)(t, ρvec)
end

###
# Diffusive stochastic evolution
###

"""
    meas(dt::Time, H; mclist=Tuple{QOp,Time,Float64}[],
                      mflist=Tuple{Function,Time,Float64}[],
                      clist=QOp[], flist=Function[])

Return increment function over a time step `dt` for diffusive monitored
evolution generated by Hamiltonian `H`, a list `mclist` of tuples specifying
constant measurement operators, and a list `mflist` of tuples specifying
time-dependent measurement operators. The tuples have the form `(m, τ, η)`,
where `m` is a (generally non-Hermitian) operator (or time-dependent function)
that specifies the measurement backaction (see below), `τ` is the timescale
of the measurement collapse, and `η` is the measurement efficiency (bounded
between 0 and 1).

These quantities are related as follows:
  - Γ = 1/(2τη) : ensemble measurement dephasing rate
  - Γm = 1/(2τ) : dephasing rate produced by averaging the collected signal
  - γ = Γ - Γm = (1-η)/(2τη) : residual dephasing rate from signal loss
  - ``m = m_o - i m_p``
  - ``m_o = (m + m')/2``  : measured observable (Hermitian part of `m`)
  - ``m_p = i(m - m')/2`` : phase-backaction generator (anti-Hermitian part of `m`)

The backaction takes the form of a three-step process:
  1. Sample `r` from a Gaussian distribution with mean ``⟨m_o⟩`` and variance τ/dt
  2. Conditionally update the state with a purity-preserving Kraus operator:
  - ``M_r = (dt/2πτ)^{1/4} e^{ -i m_p r dt/(2τ) - dt(r - m_o)^2/(4τ) }``
  3. Apply residual Lindblad dephasing evolution, including
  - `m` with rate ``γ = (1-η)/(2τη)``
  - the dephasing operators in `clist` and `flist`
  - natural Hamiltonian evolution `H`

Uses the "jump no-jump" method to efficiently approximate the residual
Lindblad dephasing as a composition of Hamiltonian evolution, jumps,
and no-jump informational backaction. Assumes no time-dependence in `alist`,
and small dt.  [Physical Review A **92**, 052306 (2015)]

### Returns:
  - (t::Time, ρ(t)::QOp) -> (ρ(t+dt)::QOp, rlist::Float64...)

"""
function meas(dt::Time, H; mclist=Tuple{QOp,Time,Float64}[], mflist=Tuple{Any,Time,Float64}[], clist=QOp[], flist=Function[])
    # Assemble readout generating functions and Kraus operators
    ros = Function[]
    gks = Function[]
    as = QOp[]
    afs = Function[]
    for (m, τ, η) in mclist
        push!(ros, readout(dt, m, τ))
        push!(gks, gausskraus(dt, m, τ))
        # For inefficient measurements append to Lindblad dephasing
        if abs(1.0 - η) > 1e-4
            push!(as, sqrt((1.0-η)/(4*τ*η)).*m)
        end
    end
    for (m, τ, η) in mflist
        push!(ros, readout(dt, m, τ))
        push!(gks, gausskraus(dt, m, τ))
        # For inefficient measurements append to Lindblad dephasing
        if abs(1.0 - η) > 1e-4
            push!(afs, t -> sqrt((1.0-η)/(4*τ*η)).*m(t))
        end
    end
    # Assemble total Lindblad dephasing using jump/no-jump method
    l = lind(dt, H, clist=vcat(clist, as), flist=vcat(flist, afs))
    # Increment that samples each readout, applies all Kraus operators
    # then applies Lindblad dephasing (including Hamiltonian evolution)
    (t::Time, ρ) -> let rs = map(ro -> ro(t, ρ), ros),
                        gs = map(z -> z[2](t, z[1]), zip(rs,gks)),
                        ρ1 = foldr(sand, gs; init=ρ);
                    (l(t, ρ1/tr(ρ1)), rs...) end
end

"""
    readout(dt::Time, m::QOp, τ::Time)
    readout(dt::Time, m::Function, τ::Time)

Return a function that accepts a (time t, state ρ) tuple and returns a random
number that is Gaussian-distributed about the mean ⟨(m + m')/2⟩ with standard
deviation σ = sqrt(τ/dt). ```m``` may be a function of time t that returns a QOp.
"""
function readout(dt::Time, m::QOp, τ::Time)
    mo = (m .+ m')./2
    σ = sqrt(τ/dt)
    (t::Time, ρ) -> σ*randn() + real(expect(ρ, mo))
end
function readout(dt::Time, m::Function, τ::Time)
    σ = sqrt(τ/dt)
    (t::Time, ρ) -> let mo = (m(t) .+ m(t)')./2;
                        σ*randn() + real(expect(ρ, mo)) end
end

"""
    gausskraus(dt::Time, m::QOp, τ::Time)

Return a function that accepts a Gaussian random sample `r` and returns a
Kraus operator for the associated state back-action. Note that `m` is an
operator of the form `m = m_o - i m_p`, where `m_o` is the Hermitian
observable being monitored, and `m_p` is the generator for phase-backaction.

The generated Kraus operator has the following Gaussian form:

  - ``M_r = (dt/2πτ)^(1/4) e^{ -i m_p r dt/(2τ) - dt(r - m_o)^2/(4τ) }``

However, the constant factors have been eliminated in anticipation of
subsequent state renormalization, to simply yield:

  - ``M_r = e^{ m r dt/(2τ) - m_o^2 dt/(4τ) }``

"""
function gausskraus(dt::Time, m::QOp, τ::Time)
    mo = (m .+ m') ./ 2
    mo2 = Array(mo^2 ./ 2)
    mf = Array(m)
    v = dt/(2*τ)
    (t::Time, r) -> sparse(exp((r*v).*mf .- v.*mo2))
end
function gausskraus(dt::Time, m::Function, τ::Time)
    v = dt/(2*τ)
    (t::Time, r) -> let mo = (m(t) .+ m(t)') ./ 2
                        mo2 = Array(mo^2 ./ 2)
                        mf = Array(m(t))
                        sparse(exp((r*v) .* mf .- v .* mo2)) end
end

###
# Crude trajectory integrator
###

# Return trajectory array [f1(now), f2(now), ...]
"""
    trajectory(inc::Function, init,
                    tspan::Tuple{Time,Time}, fs::Function...;
                    dt::Time=1/10^4, points::Int=1000, verbose=true, readout=true)

Compute time-stepped trajectory, starting from state `init`, incrementing with `inc`
by time step `dt` for `t` in `tspan`, keeping `points` intermediate values
of `f(ρ(t))` for each `f` in `fs`. If the trajectory is stochastic, optionally store
the simulated readout as well.

### Returns:
  - Array of Trajectories for each function (and readout) value:
    [Trajectory(f(ρ(t)))..., Trajectory(r(t))...]

"""
@inline function trajectory(inc::Function, init,
                    tspan::Tuple{Time,Time}, fs::Function...;
                    dt::Time=1/10^4, points::Int=1000, verbose=true, readout=true)
    # Simulate and collect data
    s = simulate(inc, init, tspan, fs...,
                 dt=dt, points=points, verbose=verbose, readout=readout)
    # Return interpolated trajectory objects
    ts = range(first(tspan), length=length(first(s)), stop=last(tspan))
    map(v -> Trajectory{typeof(first(v))}(ts, v), s)
end

"""
    ensemble(n::Integer, inc::Function, init,
             tspan::Tuple{Time,Time}, fs::Function...;
             dt::Time=1/10^4, points::Int=1000, verbose=true, readout=true)

Compute time-stepped ensemble of `n` trajectories, starting from state `init`,
incrementing with `inc` by time step `dt` for `t` in `tspan`, keeping `points`
intermediate values of `f(ρ(t))` for each `f` in `fs`. If the trajectory is
stochastic, optionally store the simulated readouts as well.

### Returns:
  - Array of Ensembles for each function (and readout) value:
    [Ensemble(f(ρ(t)))..., Ensemble(r(t))...]

"""
@inline function ensemble(n::Integer, inc::Function, init,
                  tspan::Tuple{Time,Time}, fs::Function...;
                  dt::Time=1/10^4, points::Int=1000, verbose=true, readout=true)
    # Information if desired
    if verbose
        pinit = inc(first(tspan), init)
        N = Int(fld(abs(last(tspan)-first(tspan)), dt))
        @info("Trajectories: ",n,", steps each: ",N,
             ", points each: ",min(N,points),", values each = ", length(fs))
        readout && typeof(pinit) <: Tuple &&
            @info("Readouts: values each = ", length(pinit)-1)
    end

    inittime=time()
    # Simulate and collect data
    # Uses pmap to optionally support parallelization
    s = hcat(pmap( _ -> simulate(inc, init, tspan, fs...,
                      dt=dt, points=points, verbose=false, readout=readout),
                1:n)...)
    elapsed = time() - inittime

    # Return ensemble of each stored value
    ts = range(first(tspan), length=length(first(s)), stop=last(tspan))
    out = collect(Ensemble(ts, hcat(s[i,:]...)) for i in 1:first(size(s)))

    # Info if desired
    if verbose
        steps = length(first(out))
        @info("Time elapsed: ",elapsed," s, Steps: ",steps,
             ", Steps per second: ",steps/elapsed)
    end
    out
end

function simulate(inc::Function, init,
                    tspan::Tuple{Time,Time}, fs::Function...;
                    dt::Time=1/10^4, points::Int=1000, verbose=true, readout=true)
    # Constants for integration
    t0 = first(tspan)              # Initial time
    tmax = last(tspan)             # Final time
    ts = t0:dt:tmax                      # simulated time range
    N = length(ts)                 # total num of points
    if points > N
        Ns = N                     # reset points if needed
    else
        Ns = points                # stored points
    end
    Nf = length(fs)                # stored f values per point
    Nl = Int(cld(N-1, points))     # steps per stored point
    Nldt = Nl*dt                   # time-step per stored point

    # Preallocate trajectory arrays for speed
    traj = [let fi=f(init); Array{typeof(fi)}(undef,Ns) end for f in fs]
    # Functions to update values
    function update!(i::Int, ρ)
        for k in 1:Nf
            @inbounds traj[k][i] = fs[k](ρ)
        end
    end
    next(t::Time, ρ) = inc(t, ρ)

    # Modify preallocation and updates for stochastic readout
    # Trial point to test for readout
    pinit = inc(first(tspan), init)
    if readout && typeof(pinit) <: Tuple
        # Stored readout variables per point
        Nr = length(pinit) - 1
        # Preallocate readout arrays for speed
	rtraj = [Array{typeof(r)}(undef,Ns) for r in pinit[2:end]]
        append!(traj, rtraj)
	# Specialize updates to handle readout tuples
        function update!(i::Int, t::Tuple)
            update!(i, first(t))
            for k in 1:Nr
                @inbounds traj[Nf+k][i] = t[k+1]
            end
        end
        next(t::Time, tup::Tuple) = inc(t, first(tup))
    else
        Nr = 0
    end

    # Seed loop
    verbose && @info("Trajectory: steps = ",N-1,", points = ",Ns,", values = ",Nf)
    verbose && readout && Nr > 0 && @info("Readout: values = ",Nr)
    inittime = time()
    now = init
    tnow = t0
    update!(1, now)
    # loop
    for i in 2:Ns
        # inner loop without storage
        for _ in 1:Nl
            tnow += dt
            now = next(tnow, now)
        end
        # store point
        update!(i, now)
    end
    elapsed = time() - inittime
    # Performance summary
    verbose && @info("Time elapsed: ",elapsed," s, Steps per second: ",(N-1)/elapsed)
    traj
end


text:
"""
NueralNetwork module

MIT License
Copyright (c) 2019 amenoyoya
"""

# 数値微分関数
numeric_diff(f, x) = begin
    h = 1e-4 # 10^(-4)
    (f(x + h) - f(x - h)) / 2h # 中心差分から微分値を計算
end

# 数値微分による勾配計算
## 関数 f(Float64配列), Float64配列 x -> Float64配列 勾配
numeric_gradient(f, x::Array{Float64,1})::Array{Float64,1} = begin
    h = 1e-4 # 10^(-4)
    grad = Array{Float64, 1}(undef, length(x)) # xと同じ長さの配列 [undef, undef, ...] を生成
    # 各変数ごとの数値微分を配列にまとめる
    for i in 1:length(x)
        # 指定indexの変数に対する中心差分を求める
        org = x[i]
        x[i] = org + h
        f1 = f(x) # f([..., x[i] + h, ...])
        x[i] = org - h
        f2 = f(x) # f([..., x[i] - h, ...])
        grad[i] = (f1 - f2) / 2h # i番目の変数に対する数値微分
        x[i] = org # x[i]の値をもとに戻す
    end
    return grad
end

# 数値微分による勾配計算（バッチ対応版）
## 関数(Array{Float64,2})::Float64, 入力値::Array{Float64,2} -> 勾配::Array{Float64,2}
numeric_gradient(f, x::Array{Float64,2})::Array{Float64,2} = begin
    h = 1e-4 # 10^(-4)
    grad = Array{Float64, 2}(undef, size(x, 1), size(x, 2)) # xと同じ次元の行列を生成
    # 各変数ごとの数値微分を行列にまとめる
    for row in 1:size(x, 1), col in 1:size(x, 2)
        # 指定indexの変数に対する中心差分を求める
        org = x[row, col]
        x[row, col] = org + h
        f1 = f(x) # f([..., x[row, col] + h, ...]) -> Float64
        x[row, col] = org - h
        f2 = f(x) # f([..., x[row, col] - h, ...]) -> Float64
        grad[row, col] = (f1 - f2) / 2h # (row, col)番目の変数に対する数値微分
        x[row, col] = org # x[i]の値をもとに戻す
    end
    return grad
end

# 勾配降下法
## 損失関数 f, ベクトル x, 学習率 lr, 試行回数 ste_num -> 最適化ベクトル
gradient_descent(f, x::Array{Float64,1}, lr::Float64=0.01, step_num::Int=100)::Array{Float64,1} = begin
    step_num > 0 ? gradient_descent(f, x - lr * numeric_gradient(f, x), lr, step_num - 1) : x
end

max_forks_repo_name: amenoyoya/scala-ml


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV -t 6
#     language: julia
#     name: julia-1.8-t4
# ---

# %%
a = Threads.Atomic{Int}(2)
b = a[]

# %%
using Printf
@show Threads.nthreads()
@show a = Threads.Atomic{Int}(99)
Threads.@threads for i in 1:10
    a[] = i
    @printf "(id: %2d)  i = %2d,  a[] = %2d\n" Threads.threadid() i a[]
    b = a[]
    @printf "(id: %2d)  i = %2d,  b = %2d\n" Threads.threadid() i b
end

# %%

max_stars_repo_path: 0015/Atomic example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/Atomic example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/Atomic example.jl

max_forks_repo_name: genkuroki/public


text:
export asqp, solveKKT, getLinearIndependentConstraints

"""
function asqp(G,c,Ae,be,Ai,bi,x=0*c; maxIter, tolP, out)

Find optimizer for quadratic programming problem using active set method

min_x q(x) = 0.5\*x'\*G*x + c'\*x
subject to Ax >= b

HERE: Assume that active constraints are linearly independent.

	Input:
        1. G 	 - matrix of weights for x, nxn
        2. c 	 - linear vector of weights, nx1
        3. Ai/Ae - matrix of inequality/equality constraints, mxn or pxn
        4. bi/be - inequality/equality constraints to satisfy, mx1 or px1
        5. x 	 - initial starting point if provided or default to origin [0,0,...,0], nx1

	Keyword arguments:
        6. maxIter - number of iterations for finding feasible set and direction before ending asqp
        7. tolP    - tolerance for ending search for direction and asqp
        8. out	   - print intermediate outputs if out > 1

    Output:
        1. x 	- minimizer for quadratic programming problem, nx1
        2. iter - number of iterations to find minimizer
        3. Wci  - inequality working conditions at minimizer, m'xn

"""
function asqp(G,c,Ae,be,Ai,bi,xc=0*c; maxIter::Int64=5, tolP::Float64=1e-10,out::Int64 = 2)

	# if rank(G) < size(G,2)
	# 	Ginv = x->(G+1e-6*eye(size(G,2)))\x
	# else
	# 	Ginv = x->G\x
	#  end

    # println("########## ASQP start ########### \n Ae $Ae \n be $be \n Ai $Ai \n bi $bi \n G $G \n c $c \n x0 $xc")
    Wci = find(Ai*xc-bi .< 0)
	if (length(Wci)>0 && size(Ae,1)>0)  || length(Wci)>1
    	Wci = getLinearIndependentConstraints([Ae;Ai],[vec(collect(1:size(Ae,1)));  size(Ae,1)+Wci])
		Wci = Wci[size(Ae,1)+1:end]-size(Ae,1)
	end
    Aic = Ai[Wci,:]
    # println("Aic $Aic")

    iter = 0
    while true
        if out>1
            qc = .5*xc'*G*xc + c'*xc
            nlwc = length(Wci)
            println("i=$iter,\t q(x)= $qc, \t numel(Wci)=$nlwc")
        end

        # println(" Current xc = $xc")

        # update gradient of QP
        gc  = G*xc+c

        # update pc and lambda of inequality constraints
        pc,lami = solveKKT(xc,G,gc,Ae,be,Aic,bi[Wci])
        # pc,lami = solveSchur(xc,G,gc,Ae,be,Aic,bi[Wci],Ginv)

        # println("\n ##### iter $iter ##### \n pc $pc \t lami ", lami)

        if norm(pc)<tolP # if stepsize to x* is small
            # check Lagrange multipliers
            if all(lami.>=0) # we are done!
                break #return xc,iter,Wc
            else # remove constraint
                l_, idm = findmin(lami); # find inequality working set index corresponding to smallest value of multiplier
                Wci = deleteat!(Wci,idm) # remove one inequality constraint
                # println("removed one inequality constraint, new Wci $Wci")
                Aic  = Ai[Wci,:] # update active constraints
            end
        else
            # add one constraint to active set
            t2 = map(x -> min(0,x), (bi-Ai*xc)) # keep only negative entries that satisfy inequaity constraint, t2 has only entries <= 0, necessary to ensure alpha > 0
            t1 = (Ai*pc);
            t  = t2./t1;
            # println("bi-Ai*xc ", bi-Ai*xc)
            # println("\n t = t2./t1  $t \n t1 = Ai*pc $t1")
            # println("\n t2 = min(0,bi-Ai*xc) $t2")
            # println("t1 >= 0 \t", 1*(t1.>=0))
            t[find(t1 .>= 0)] = 1; # set entries in t equal 1 if t1 entries >= 0, necessary to ensure alpha > 0
            t[Wci]    = 1;
            # println("t[Wci]=1 $t")
            alpha = min(1.0,minimum(t)); # minimum for finding array min only, find step size
            if alpha<0
                error("alpha negative");
            end
            xc = xc + alpha*pc;
            # println("xc + alpha*pc $xc")
            if alpha<1 # add one blocking constraint
                idm = findin(t, alpha); # find indicies of entries in t == alpha
				# Wci = sort([Wci;idm[1]]); # sort to ensure index for inequality constraints are above equality constraints
				Wci = [Wci;idm[1]];
                Aic  = Ai[Wci,:] # found working set
            end
        end
        # check stopping
        if iter > maxIter
            break
        end
        iter = iter+1;
    end
    return xc, iter, Wci
end

"""
 function solveKKT(xc,G,gc,Ae,be,Ai,bi)

 Solves KKT system

       | G   Ae'  Ai' | |  p   |   | -gc |
 KKT = | Ae   0   0   | | -lame| = | -he |
       | Ai   0   0   | | -lami|   | -hi |

 gc = c+G*xc, he = Ae*xc-be, hi = Ai*xc-bi

 Input:
	 1. xc	  - current point, nx1
	 1. G 	  - matrix of weights for x, nxn
	 2. gc 	  - gradient of quadratic problem, nx1
	 3. Ai/Ae - matrix of inequality/equality constraints, mxn or pxn
	 4. bi/be - inequality/equality constraints to satisfy, mx1 or px1

Output:
	1. p 	- stepsize in direction of satisfying both inequality and equality constraints to update xc
	2. lami - Lagrange multipliers associated with inequaity constraints

"""
function solveKKT(xc,G,gc,Ae,be,Ai,bi)

    hi = Ai*xc-bi # h = 0, x0 = 0, b = 0
    he = Ae*xc-be

    KKT = [G Ae' Ai';
           Ae zeros(size(Ae,1), size(Ae,1)+size(Ai,1));
           Ai zeros(size(Ai,1), size(Ae,1)+size(Ai,1)) ]

	if rank(KKT) < size(KKT,1)
	   KKT = [G+1e-6*eye(size(G,1)) Ae' Ai';
	           Ae zeros(size(Ae,1), size(Ae,1)+size(Ai,1));
	           Ai zeros(size(Ai,1), size(Ae,1)+size(Ai,1)) ]
	end
    # println("solveKKT $KKT \n rank(KKT) ", rank(KKT))
    # println("\n ###### solveKKT ###### \n rank(KKT) $(rank(KKT)) \t size(KKT) ", size(KKT), " \t cond(KKT) ", cond(KKT), "\n xc $xc \n gc $gc \n G $G \n cond(G) ", cond(G)) #1e6 or larger

    rhs = [-gc; -he; -hi];

    try
        xsol = KKT\rhs;
        # println("inside try xsol $xsol")
        catch E
        println("ERROR!!!! \n rank(KKT) $(rank(KKT)) \t size(KKT) ", size(KKT), " \t cond(KKT) ", cond(KKT), "\n xc $xc \n gc $gc \n Ai $Ai \n bi $bi \n be $be \n Ae $Ae \n G $G \n rank(G) ", rank(G)) #1e6 or larger
        # JLD.save("asqp error case.jld","G", G, "xc", xc, "gc", gc, "Ai", Ai, "bi", bi, "Ae", Ae, "be", be)
        # MAT.matwrite("test.mat", Dict(
                #"G" => G,
                #"xc" => xc,
                #"gc" => gc,
                #"Ai" => Ai,
                #"bi" => bi,
                #"Ae" => Ae,
                #"be" => be,
                #"KKT" => KKT
                #))
        throw(E)
    end

    xsol = KKT\rhs
    p = xsol[1:length(gc)]
    lame = -xsol[length(gc)+(1:size(Ae,1))];
    lami = -xsol[length(gc)+size(Ae,1)+1:end];

    # println("\n p and lam $xsol \n Ae $Ae \n Ai ", Ai)

    return p,lami
end
#
#
# """
#  solveSchur(xc,G,gc,Ae,be,Ai,bi,Ginv=x->G\\x)
#
#  Solves KKT system
#
#        | G   Ae'  Ai' | |  p   |   | -gc |
#  KKT = | Ae   0   0   | | -lame| = | -he |
#        | Ai   0   0   | | -lami|   | -hi |
#
#  using Schur complement
#
#  he = Ae*xc-be, gc = c+G*xc, hi = Aic*xc-bic
#
# """
# function solveSchur(xc,G,gc,Ae,be,Ai,bi,Ginv=x->G\x)
#     # println("##### solve Schur #######")
#
#     A = [Ae;Ai]
#     b = [be;bi]
#     elen = length(be)
#     #println("A $A \t x $xc \t b $b \n G $G \t Wc $Wc")
#
#     h   = A*xc-b
#     if isempty(A)
#         lami = zeros(0)
#         p    = Ginv(-gc)
#     else
#         lam  = (A*(Ginv(A')))\(A*Ginv(gc)-h)  # note that lam has a minus sign in KKT system
# # println("lam = $lam")
# # println("A'*lam-gc=$(A'*lam-gc)")
#         p    = Ginv(A'*lam-gc)
#         lami = lam[elen+1:end]
#     end
#     return p, lami
# end

"""
function getLinearIndependentConstraints(A,Wc)

Find linearly independent rows of A (working set of constraints are independent of each other)

Input:
	1. A  - set of constraints in matrix form, matrix of integers
	2. Wc - current working set of constraints, array of integers

Output:
	1. Wc - updated working set of constraints, array of integers

"""

function getLinearIndependentConstraints(A,Wc)
    if length(Wc)==1 # only working constraint
        return Wc
    end
    # save('/tmp/getLinearIndependentConstraints.mat')

    # println("##### getLinearIndependentConstraints ###### \n A $A \n Wc $Wc")
    # make matrix smaller using its sparsity
    A  = A[Wc,:];
    id = find(0.< sum(abs.(A),1)); # sum along colums of A > 0
    # println("id $id \n A $A")
    A = A[:,id[:]]; # restrict to nonzero columns
    # println("A ", A)

    # check if rank is full
    if maximum(size(A)) <=3000
        if rank(full(A)) == length(Wc)
            #println("full rank A")
            return Wc
        end
    end

    debit = 1;
    credit = 2:size(A,1);

    # eliminate linear dependent constraints row by row
    Aac = A[debit,:];
    # println("Aac before credit ", Aac)

    while !(isempty(credit))
        # tic()
        At  = A[credit,:]';
        # println("At ", At, "\n Aac ", size(Aac)) #row vector always as columns
        xx = (Aac)\At; # solve system after adding new row (constraints)
        # println("xx ", xx)
        res = full(Aac*xx-At); # convert sparse to full matrix
        # println("res ", res)
        nres = sum(res.^2,1); # residual sum of squares along column
        mx,ii = findmax(nres); # find column with largest RSS
        # println("max is $mx")

        if mx>1e-3
            debit  = [debit;credit[ii]];
            # println("debit ", debit)
            credit = deleteat!(Vector(credit),ii) # remove row of constraints from inactive set
            # println("credit ", credit, "\n At[:,ii]' ", size(At[:,ii],2), "\ Aac ", Aac)
            Aac = [Aac'; At[:,ii]']; # update active set of constraints
            # println("credit ", credit, "\n Aac ", Aac)
        else
            break;
        end
        # toc()
    end

    Wc = Wc[debit]
    return Wc
end