text:
include("utils.jl")

using HybridArrays, Random
using Manifolds: default_metric_dispatch
using StaticArrays: Dynamic

Random.seed!(42)

@testset "Power manifold" begin
    Ms = Sphere(2)
    Ms1 = PowerManifold(Ms, 5)
    @test power_dimensions(Ms1) == (5,)
    @test manifold_dimension(Ms1) == 10
    @test injectivity_radius(Ms1) == π
    Ms2 = PowerManifold(Ms, 5, 7)
    @test power_dimensions(Ms2) == (5, 7)
    @test manifold_dimension(Ms2) == 70
    Ms2n = PowerManifold(Ms1, NestedPowerRepresentation(), 7)
    @test power_dimensions(Ms2n) == (7,)
    @test manifold_dimension(Ms2n) == 70

    Mr = Manifolds.Rotations(3)
    Mr1 = PowerManifold(Mr, 5)
    Mrn1 = PowerManifold(Mr, Manifolds.NestedPowerRepresentation(), 5)
    @test PowerManifold(PowerManifold(Mr, 2), 3) == PowerManifold(Mr, 2, 3)
    @test PowerManifold(Torus(2), 3) isa PowerManifold{ℝ,Torus{2}}
    @test PowerManifold(Mrn1, 3) ==
          PowerManifold(Mr, Manifolds.NestedPowerRepresentation(), 5, 3)
    @test PowerManifold(Mrn1, Manifolds.ArrayPowerRepresentation(), 3) ==
          PowerManifold(Mr, Manifolds.ArrayPowerRepresentation(), 5, 3)
    @test manifold_dimension(Mr1) == 15
    @test manifold_dimension(Mrn1) == 15
    Mr2 = PowerManifold(Mr, 5, 7)
    Mrn2 = PowerManifold(Mr, Manifolds.NestedPowerRepresentation(), 5, 7)
    @test manifold_dimension(Mr2) == 105
    @test manifold_dimension(Mrn2) == 105

    @test repr(Ms1) == "PowerManifold(Sphere(2, ℝ), 5)"
    @test repr(Mrn1) == "PowerManifold(Rotations(3), NestedPowerRepresentation(), 5)"

    @test Manifolds.allocation_promotion_function(Ms, exp, ([1],)) ==
          Manifolds.allocation_promotion_function(Ms1, exp, (1,))

    @test Ms^5 === Oblique(3, 5)
    @test Ms^(5,) === Ms1
    @test Mr^(5, 7) === Mr2

    @test is_default_metric(Ms1, PowerMetric())
    @test default_metric_dispatch(Ms1, PowerMetric()) === Val{true}()
    types_s1 = [Array{Float64,2}, HybridArray{Tuple{3,Dynamic()},Float64,2}]
    types_s2 = [Array{Float64,3}, HybridArray{Tuple{3,Dynamic(),Dynamic()},Float64,3}]

    types_r1 = [Array{Float64,3}, HybridArray{Tuple{3,3,Dynamic()},Float64,3}]

    types_rn1 = [Vector{Matrix{Float64}}]
    TEST_STATIC_SIZED && push!(types_rn1, Vector{MMatrix{3,3,Float64,9}})

    types_r2 = [Array{Float64,4}, HybridArray{Tuple{3,3,Dynamic(),Dynamic()},Float64,4}]
    types_rn2 = [Matrix{Matrix{Float64}}]

    retraction_methods = [Manifolds.PowerRetraction(ManifoldsBase.ExponentialRetraction())]
    inverse_retraction_methods =
        [Manifolds.InversePowerRetraction(ManifoldsBase.LogarithmicInverseRetraction())]

    sphere_dist = Manifolds.uniform_distribution(Ms, @SVector [1.0, 0.0, 0.0])
    power_s1_pt_dist =
        Manifolds.PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5))
    power_s2_pt_dist =
        Manifolds.PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7))
    sphere_tv_dist =
        Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0)
    power_s1_tv_dist = Manifolds.PowerFVectorDistribution(
        TangentBundleFibers(Ms1),
        rand(power_s1_pt_dist),
        sphere_tv_dist,
    )
    power_s2_tv_dist = Manifolds.PowerFVectorDistribution(
        TangentBundleFibers(Ms2),
        rand(power_s2_pt_dist),
        sphere_tv_dist,
    )

    id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0)
    power_r1_pt_dist =
        Manifolds.PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5))
    power_rn1_pt_dist = Manifolds.PowerPointDistribution(
        Mrn1,
        rotations_dist,
        [randn(Float64, 3, 3) for i in 1:5],
    )
    power_r2_pt_dist =
        Manifolds.PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7))
    power_rn2_pt_dist = Manifolds.PowerPointDistribution(
        Mrn2,
        rotations_dist,
        [randn(Float64, 3, 3) for i in 1:5, j in 1:7],
    )
    rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0)
    power_r1_tv_dist = Manifolds.PowerFVectorDistribution(
        TangentBundleFibers(Mr1),
        rand(power_r1_pt_dist),
        rotations_tv_dist,
    )
    power_rn1_tv_dist = Manifolds.PowerFVectorDistribution(
        TangentBundleFibers(Mrn1),
        rand(power_rn1_pt_dist),
        rotations_tv_dist,
    )
    power_r2_tv_dist = Manifolds.PowerFVectorDistribution(
        TangentBundleFibers(Mr2),
        rand(power_r2_pt_dist),
        rotations_tv_dist,
    )
    power_rn2_tv_dist = Manifolds.PowerFVectorDistribution(
        TangentBundleFibers(Mrn2),
        rand(power_rn2_pt_dist),
        rotations_tv_dist,
    )

    @testset "get_component, set_component!, getindex and setindex!" begin
        p1 = randn(3, 5)
        @test get_component(Ms1, p1, 1) == p1[:, 1]
        @test p1[Ms1, 1] == p1[:, 1]
        @test p1[Ms1, 1] isa Vector
        p2 = [10.0, 11.0, 12.0]
        set_component!(Ms1, p1, p2, 2)
        @test get_component(Ms1, p1, 2) == p2
        p1[Ms1, 2] = 2 * p2
        @test p1[Ms1, 2] == 2 * p2
        p1[Ms1, 2] += p2
        @test p1[Ms1, 2] ≈ 3 * p2
        p1[Ms1, 2] .+= p2
        @test p1[Ms1, 2] ≈ 4 * p2
        @test view(p1, Ms1, 1) == p1[Ms1, 1]
        @test view(p1, Ms1, 1) isa SubArray

        Msn1 = PowerManifold(Ms, Manifolds.NestedPowerRepresentation(), 5)
        pn1 = [randn(3) for _ in 1:5]
        @test get_component(Msn1, pn1, 1) == pn1[1]
        @test pn1[Msn1, 1] == pn1[1]
        @test pn1[Msn1, 1] isa Vector
        set_component!(Msn1, pn1, p2, 2)
        @test get_component(Msn1, pn1, 2) == p2
        pn1[Msn1, 2] = 2 * p2
        @test pn1[Msn1, 2] == 2 * p2
        pn1[Msn1, 2] += p2
        @test pn1[Msn1, 2] ≈ 3 * p2
        pn1[Msn1, 2] .+= p2
        @test pn1[Msn1, 2] ≈ 4 * p2
        @test view(pn1, Msn1, 1) == pn1[Msn1, 1]
        @test view(pn1, Msn1, 1) isa SubArray
    end

    @testset "ComponenException" begin
        M = PowerManifold(Sphere(2), NestedPowerRepresentation(), 2)
        p = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
        X = [[0.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
        @test_throws ComponentManifoldError is_point(M, X, true)
        @test_throws ComponentManifoldError is_vector(M, p, X, true)
    end

    @testset "power vector transport" begin
        m = PowerVectorTransport(ParallelTransport())
        p = repeat([1.0, 0.0, 0.0], 1, 5)
        q = repeat([0.0, 1.0, 0.0], 1, 5)
        X = log(Ms1, p, q)
        Y = vector_transport_to(Ms1, p, X, q, m)
        Z = -log(Ms1, q, p)
        @test isapprox(Ms1, q, Y, Z)
    end
    trim(s::String) = s[1:min(length(s), 20)]

    basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))
    for T in types_s1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_s1_pt_dist)) for _ in 1:3]
            @test injectivity_radius(Ms1, pts1[1]) == π
            basis_diag = get_basis(
                Ms1,
                pts1[1],
                DiagonalizingOrthonormalBasis(log(Ms1, pts1[1], pts1[2])),
            )
            basis_arb = get_basis(Ms1, pts1[1], DefaultOrthonormalBasis())
            test_manifold(
                Ms1,
                pts1;
                test_reverse_diff=true,
                test_musical_isomorphisms=true,
                test_injectivity_radius=false,
                test_default_vector_transport=true,
                test_project_point=true,
                test_project_tangent=true,
                vector_transport_methods=[
                    PowerVectorTransport(ParallelTransport()),
                    PowerVectorTransport(SchildsLadderTransport()),
                    PowerVectorTransport(PoleLadderTransport()),
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                test_vee_hat=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[power_s1_pt_dist],
                tvector_distributions=[power_s1_tv_dist],
                basis_types_to_from=(basis_diag, basis_arb, basis_types...),
                rand_tvector_atol_multiplier=6.0,
                retraction_atol_multiplier=12.0,
                is_tangent_atol_multiplier=12.0,
                exp_log_atol_multiplier=2 * prod(power_dimensions(Ms1)),
                test_inplace=true,
            )
        end
    end
    for T in types_s2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_s2_pt_dist)) for _ in 1:3]
            test_manifold(
                Ms2,
                pts2;
                test_reverse_diff=true,
                test_musical_isomorphisms=true,
                test_injectivity_radius=false,
                test_vee_hat=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[power_s2_pt_dist],
                tvector_distributions=[power_s2_tv_dist],
                rand_tvector_atol_multiplier=6.0,
                retraction_atol_multiplier=12,
                is_tangent_atol_multiplier=12.0,
                exp_log_atol_multiplier=2 * prod(power_dimensions(Ms2)),
                test_inplace=true,
            )
        end
    end

    for T in types_r1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_r1_pt_dist)) for _ in 1:3]
            test_manifold(
                Mr1,
                pts1;
                test_reverse_diff=false,
                test_injectivity_radius=false,
                test_musical_isomorphisms=true,
                test_vee_hat=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[power_r1_pt_dist],
                tvector_distributions=[power_r1_tv_dist],
                basis_types_to_from=basis_types,
                rand_tvector_atol_multiplier=8.0,
                retraction_atol_multiplier=12,
                is_tangent_atol_multiplier=12.0,
                exp_log_atol_multiplier=2e2 * prod(power_dimensions(Mr2)),
                test_inplace=true,
            )
        end
    end

    for T in types_rn1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_rn1_pt_dist)) for _ in 1:3]
            test_manifold(
                Mrn1,
                pts1;
                test_reverse_diff=false,
                test_injectivity_radius=false,
                test_musical_isomorphisms=true,
                test_vee_hat=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[power_rn1_pt_dist],
                tvector_distributions=[power_rn1_tv_dist],
                basis_types_to_from=basis_types,
                rand_tvector_atol_multiplier=5.0,
                retraction_atol_multiplier=12,
                is_tangent_atol_multiplier=12.0,
                exp_log_atol_multiplier=4e2 * prod(power_dimensions(Mrn1)),
                test_inplace=true,
            )
        end
    end
    for T in types_r2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_r2_pt_dist)) for _ in 1:3]
            test_manifold(
                Mr2,
                pts2;
                test_reverse_diff=false,
                test_injectivity_radius=false,
                test_musical_isomorphisms=true,
                test_vee_hat=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[power_r2_pt_dist],
                tvector_distributions=[power_r2_tv_dist],
                rand_tvector_atol_multiplier=8.0,
                retraction_atol_multiplier=12,
                is_tangent_atol_multiplier=12.0,
                exp_log_atol_multiplier=4e3 * prod(power_dimensions(Mr2)),
                test_inplace=true,
            )
        end
    end
    for T in types_rn2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_rn2_pt_dist)) for _ in 1:3]
            test_manifold(
                Mrn2,
                pts2;
                test_reverse_diff=false,
                test_injectivity_radius=false,
                test_musical_isomorphisms=true,
                test_vee_hat=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[power_rn2_pt_dist],
                tvector_distributions=[power_rn2_tv_dist],
                rand_tvector_atol_multiplier=8.0,
                retraction_atol_multiplier=12,
                is_tangent_atol_multiplier=12.0,
                exp_log_atol_multiplier=4e3 * prod(power_dimensions(Mrn2)),
                test_inplace=true,
            )
        end
    end

    @testset "Power manifold of Circle" begin
        pts_t = [[0.0, 1.0, 2.0], [1.0, 1.0, 2.4], [0.0, 2.0, 1.0]]
        MT = PowerManifold(Circle(), 3)
        @test representation_size(MT) == (3,)
        test_manifold(
            MT,
            pts_t;
            test_reverse_diff=false,
            test_forward_diff=false,
            test_injectivity_radius=false,
            test_musical_isomorphisms=true,
            test_vee_hat=true,
            retraction_methods=retraction_methods,
            inverse_retraction_methods=inverse_retraction_methods,
            rand_tvector_atol_multiplier=5.0,
            retraction_atol_multiplier=12,
            is_tangent_atol_multiplier=12.0,
            exp_log_atol_multiplier=1.0,
            test_inplace=true,
        )
    end

    @testset "Basis printing" begin
        p = hcat([[1.0, 0.0, 0.0] for i in 1:5]...)
        Bc = get_basis(Ms1, p, DefaultOrthonormalBasis())
        @test sprint(show, "text/plain", Bc) == """
        DefaultOrthonormalBasis(ℝ) for a power manifold
        Basis for component (1,):
        $(sprint(show, "text/plain", Bc.data.bases[1]))
        Basis for component (2,):
        $(sprint(show, "text/plain", Bc.data.bases[2]))
        Basis for component (3,):
        $(sprint(show, "text/plain", Bc.data.bases[3]))
        Basis for component (4,):
        $(sprint(show, "text/plain", Bc.data.bases[4]))
        Basis for component (5,):
        $(sprint(show, "text/plain", Bc.data.bases[5]))
        """
    end

    @testset "Power manifold of Circle" begin
        pts_t = [[0.0, 1.0, 2.0], [1.0, 1.0, 2.4], [0.0, 2.0, 1.0]]
        MT = PowerManifold(Circle(), 3)
        @test representation_size(MT) == (3,)
        @test pts_t[1][MT, 2] == 1.0
        @test HybridVector{3}(pts_t[1])[MT, 2] == 1.0
        test_manifold(
            MT,
            pts_t;
            test_reverse_diff=false,
            test_forward_diff=false,
            test_injectivity_radius=false,
            test_musical_isomorphisms=true,
            retraction_methods=retraction_methods,
            inverse_retraction_methods=inverse_retraction_methods,
            rand_tvector_atol_multiplier=5.0,
            retraction_atol_multiplier=12,
            is_tangent_atol_multiplier=12.0,
            test_inplace=true,
        )
    end

    @testset "Atlas & Induced Basis" begin
        M = PowerManifold(Euclidean(2), NestedPowerRepresentation(), 2)
        p = [zeros(2), ones(2)]
        X = [ones(2), 2 .* ones(2)]
        A = RetractionAtlas()
        a = get_parameters(M, A, p, p)
        p2 = get_point(M, A, p, a)
        @test all(p2 .== p)
    end
end


text:
# Autogenerated wrapper script for OpenSpecFun_jll for x86_64-w64-mingw32-libgfortran4
export libopenspecfun

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("OpenSpecFun")
JLLWrappers.@declare_library_product(libopenspecfun, "libopenspecfun.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libopenspecfun,
        "bin\\libopenspecfun.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
"""
    mpdist(A, B, m, d = ZEuclidean, k = ((length(A) + length(B)) - 2m) ÷ 20)

The MP distance between `A` and `B` using window length `M` and returning the `k`th smallest value.
"""
function mpdist(A,B,m,d=ZEuclidean(),k=(length(A)+length(B)-2m)÷20)
    p1 = matrix_profile(A,B,m,d)
    p2 = matrix_profile(B,A,m,d)
    partialsort([p1.P; p2.P], k)
end

"""
    mpdist(MP::AbstractArray, th::Real, N::Int)

The MP distance given a precomputed matrix profile, calculating `k` as `th*N` where `N` is the length of the data used to create `MP`.
"""
function mpdist(MP::AbstractArray, th::Real, N::Int)
    k = ceil(Int, th*N)
    if isempty(MP)
        return typemax(eltype(MP))
    elseif length(MP) >= k
        return partialsort(MP, k)
    else
        return maximum(filter(isfinite, MP)) # Do not inplace this
    end
end

"""
    mpdist_profile(T::AbstractVector, S::Int, m::Int, [d])

All MP distance profiles between subsequences of length `S` in `T` using internal window length `m`.
"""
function mpdist_profile(T::AbstractVector{TT},S::Int, m::Int, d::DT = ZEuclidean()) where {TT,DT}
    S >= m || throw(ArgumentError("S should be > m"))
    SlidingDistancesBase.DSP.FFTW.set_num_threads(1)
    n = length(T)
    pad = S * ceil(Int, n / S) - n
    append!(T,zeros(TT, pad)) # vcat was not type stable
    N = S-m+1
    D = Matrix{float(TT)}(undef, length(T)-m+1, N)
    sliding_means = similar(D)
    m_profile = similar(D, 2N)

    prog = Progress((n-S)÷S, dt=1, desc="MP dist profile")
    map(1:S:n-S) do i
        dp = mpdist_profile(getwindow(T,S,i), T, m, d, D, sliding_means, m_profile)
        next!(prog)
        dp
    end
end

"""
    mpdist_profile(A::AbstractVector, B::AbstractVector, m::Int, [d = ZEuclidean()])

MP distance profile between two time series using window length `m`.
"""
function mpdist_profile(
    B::AbstractVector,
    A::AbstractVector,
    m::Int,
    d = ZEuclidean(),
    D = Matrix{float(eltype(A))}(undef, lastlength(A)-m+1, lastlength(B)-m+1),
    sliding_means = similar(D),
    m_profile = similar(D, 2size(D,2))
)
    mpdistmat!(D,A,B,m,d)
    th = floattype(B)(0.05)
    lA = lastlength(A)
    lB = lastlength(B)
    lA >= lB ||
    throw(ArgumentError("The first argument should be shorter than the second."))

    rows, cols = size(D)
    right_margs = vec(minimum(D, dims = 2))
    for i = 1:cols
        sliding_mean!(sliding_means[!, i], D[!, i], cols)
    end

    l = lA - lB + 1
    map(1:l) do i
        m_profile[1:cols] .= @view sliding_means[i+cols÷2, :] # left marg
        copyto!(m_profile, cols+1, right_margs, i, cols)
        mpdist(m_profile, th, 2 * lastlength(B))
    end
end

function mpdistmat!(D, A::AbstractVector, B::AbstractVector, m::Int, d)
    N = lastlength(B)-m +1
    Threads.@threads for i = 1:N # Not thread safe for more than one FFTW thread
        distance_profile!(D[!,i], ZEuclidean(), getwindow(B, m, i), A)
    end
    D
end

"""
    snippets(T, k, S; m = max(S ÷ 10, 4))

Summarize time series `T` by extracting `k` snippets of length `S`
The parameter `m` controls the window length used internally.
"""
function snippets(T, k, S, d=ZEuclidean(); m = max(S÷10, 4), th=S)
    D      = mpdist_profile(T,S,m,d)
    INF    = typemax(floattype(T))
    Q      = fill(INF, length(D[1]))
    n      = lastlength(T)
    minI   = 0
    onsets = zeros(Int, k)
    snippet_profiles = similar(D, k)
    for j = 1:k
        minA = INF
        for i = 1:length(D)
            A = sum(min.(D[i], Q))
            if A < minA #&& all(abs.(((i-1)*S+1) .- onsets[1:j-1]) .> th)
                minA = A
                minI = i
            end
        end
        snippet_profiles[j] = D[minI]
        Q .= min.(D[minI], Q)
        onsets[j] = (minI-1)*S+1
    end
    tot_min = reduce((x,y)->min.(x,y), snippet_profiles, init=snippet_profiles[1])
    Cfracs = map(1:k) do j
        spj = snippet_profiles[j]
        f = count(spj .== tot_min)
        Cfrac = f/length(spj)
    end

    profile = Profile(T,tot_min,Int[],S,nothing)
    snips = Subsequence(profile, onsets, "Snippet")
    Snippets(T,tot_min, snips, Cfracs, snippet_profiles, S)
end

struct Snippets
    T
    minprofile
    snippets
    Cfracs
    profiles
    S
end

segment(s::Snippets) = [findfirst(==(s.minprofile[i]), getindex.(s.profiles, i)) for i in eachindex(s.minprofile)]


text:
module Pith

using StatsFuns

export pith
include("kernels.jl")

function pith{T}(x, X::Vector{T}; m::Kernel = Gaussian())
    n = length(X)
    ĝ = zeros(T, n)
    h = (1.364 * std(X))/(pi*n^(1/5))
    @inbounds for i in 1:length(x)
        xi = x[i]
        for j = 1:length(X)
            ĝ[i] += m((X[j]-xi)/h)
        end
    end
    ĝ/(n*h)
end

end # module


text:
using Clang.cindex
using Base.Test

top = cindex.parse_header("cxx/cxxbasic.h"; cplusplus = true)

funcs = cindex.search(top, "func")
@test length(funcs) == 1
f = funcs[1]

#@test map(spelling, cindex.function_args(f)) == ASCIIString["Int", "Int"]
@test map(spelling, cindex.function_args(f)) == ASCIIString["x", "y"]
@test isa(return_type(f), IntType)

defs = cindex.function_arg_defaults(f)
@test defs == (nothing, -10)

funcs = cindex.search(top, "func2")
@test length(funcs) == 1
f = funcs[1]

@test all(map(x->x.text, function_return_modifiers(f)) .== ["const", ])


text:
#**************************************************************************************
# test_calculations.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file contains all the required function to test calculate the Geometrical efficiency.
#
#**************************************************************************************

using GeoEfficiency
using Base.Test

info("Statrting `geoEff` test...")
info("Spectial case for cylinderical detector; very restrict test")
for raduis=[1,2,3,4,5,6,7,8,9,10.1,10.5,10.6]
	detector = CylDetector(raduis)
	@test_approx_eq 	geoEff(detector, Point(0)) 							0.5
	@test_approx_eq 	geoEff(detector, Point(0, raduis-eps(raduis)))		0.5
	@test_approx_eq 	geoEff(detector, Point(0, -raduis+eps(raduis)))		0.5
	@test_approx_eq 	geoEff(detector, Point(0, raduis/2-eps(raduis)))	0.5
	@test_approx_eq 	geoEff(detector, Point(0, -raduis/2+eps(raduis)))	0.5

	@test_approx_eq 	geoEff(detector, Point(eps()))					0.5
	@test_approx_eq 	geoEff(detector, Point(eps(), raduis-eps(raduis)))		0.5
	@test_approx_eq 	geoEff(detector, Point(eps(), -raduis+eps(raduis)))		0.5
	@test_approx_eq 	geoEff(detector, Point(eps(), raduis/2-eps(raduis)))	0.5
	@test_approx_eq 	geoEff(detector, Point(eps(), -raduis/2+eps(raduis)))	0.5

end #for
info("test pass")

info("Spectial case for well detector")
for raduis=[1,2,3,4,5,6,7,8,9,10.1,10.5,10.6], height=[1,2,3,4,5,6,7,8,9,10.1,10.5,10.6], k=[1.1,2,3,4,5,6,7,8,9,10.1,10.5,10.6]
	welldepth = height/k		# k > 1
	holeradius = raduis/k		# k > 1
	detector = WellDetector(raduis,height,holeradius, welldepth)

	@test_approx_eq 	geoEff(detector, Point(welldepth)) 									0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth, holeradius-eps(holeradius)))		0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth, -holeradius+eps(holeradius)))		0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth, holeradius/2-eps(holeradius)))	0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth, -holeradius/2+eps(holeradius)))	0.5

	@test_approx_eq 	geoEff(detector, Point(welldepth+eps(welldepth)))									0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth+eps(welldepth), holeradius-eps(holeradius)))		0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth+eps(welldepth), -holeradius+eps(holeradius)))		0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth+eps(welldepth), holeradius/2-eps(holeradius)))		0.5
	@test_approx_eq 	geoEff(detector, Point(welldepth+eps(welldepth), -holeradius/2+eps(holeradius)))	0.5

end #for
info("test pass")

info("Statrting scaling test cylinderical detector with point source...")
for raduis=[1,2,3,4,5,6,7,8,9,10.1,10.5,10.6], j=2:100, k=2:100
	detector = CylDetector(raduis); detectork = CylDetector(raduis*k)

	@test_approx_eq 	geoEff(detector, Point(raduis/j)) 				geoEff(detectork, Point(k*raduis/j))

	@test_approx_eq 	geoEff(detector, Point(raduis*j, raduis/j))		geoEff(detectork, Point(k*raduis*j, k*raduis/j))
	println("scale ",k," ; test number ",j-1, " | sucsseful")
end #for
info("test pass")


text:
# A server will report most errors to the client over a Socket, but
# errors in ZMQ state can't be reported because the socket may be
# corrupted. Therefore, we need an exception type for errors that
# should be reported locally.
struct StateError <: Exception
    msg::AbstractString
end
show(io, thiserr::StateError) = print(io, "ZMQ: ", thiserr.msg)

# Basic functions
zmq_errno() = ccall((:zmq_errno, libzmq), Cint, ())
function jl_zmq_error_str()
    errno = zmq_errno()
    c_strerror = ccall((:zmq_strerror, libzmq), Ptr{UInt8}, (Cint,), errno)
    if c_strerror != C_NULL
        strerror = unsafe_string(c_strerror)
        return strerror
    else
        return "Unknown error"
    end
end


text:
## *
grad!(dy::TAny, ::Val{1}, op::Call{typeof(*), Tuple{TReal, TReal}}) = dy * op.args[2]
grad!(dy::TAny, ::Val{2}, op::Call{typeof(*), Tuple{TReal, TReal}}) = dy * op.args[1]

## /
grad!(dy::TAny, ::Val{1}, op::Call{typeof(/), Tuple{TReal, TReal}}) = dy / op.args[2]
grad!(dy::TAny, ::Val{2}, op::Call{typeof(/), Tuple{TReal, TReal}}) =
    begin x, y = op.args; -x * dy / (y * y) end

## +
function grad!(dy::TAny, ::Val{1}, op::Call{typeof(+), Tuple{TReal, TReal}})
    x, y = op.args
    return record!(dy.tape, Assign, dy)
end
function grad!(dy::TAny, ::Val{2}, op::Call{typeof(+), Tuple{TReal, TReal}})
    x, y = op.args
    return record!(dy.tape, Assign, dy)
end

## -
function grad!(dy::TAny, ::Val{1}, op::Call{typeof(-), Tuple{TReal, TReal}})
    x, y = op.args
    return record!(dy.tape, Assign, dy)
end
function grad!(dy::TAny, ::Val{2}, op::Call{typeof(-), Tuple{TReal, TReal}})
    x, y = op.args
    return record!(dy.tape, Call, -, (dy,))
end

## ^
function grad!(dy::TAny, ::Val{1}, op::Call{typeof(^), Tuple{TReal, TReal}})
    x, y = op.args
    return y * x ^ (y-1) * dy
end
function grad!(dy::TAny, ::Val{2}, op::Call{typeof(^), Tuple{TReal, TReal}})
    x, y = op.args
    return log(x) * x ^ y * dy  # this seems to be broken: log of negative number is undefined
end

## single arg functions
grad!(dy::TAny, ::Val{1}, op::Call{typeof(-), Tuple{TReal}}) = -dy
grad!(dy::TAny, ::Val{1}, op::Call{typeof(sin), Tuple{TReal}}) = cos(op.args[1]) * dy
grad!(dy::TAny, ::Val{1}, op::Call{typeof(cos), Tuple{TReal}}) = -sin(op.args[1]) * dy
grad!(dy::TAny, ::Val{1}, op::Call{typeof(log), Tuple{TReal}}) = dy / op.args[1]
grad!(dy::TAny, ::Val{1}, op::Call{typeof(exp), Tuple{TReal}}) = exp(op.args[1]) * dy
grad!(dy::TAny, ::Val{1}, op::Call{typeof(abs), Tuple{TReal}}) = sign(op.args[1]) * dy
grad!(dy::TAny, ::Val{1}, op::Call{typeof(abs2), Tuple{TReal}}) = 2 * op.args[1] * dy
grad!(dy::TAny, ::Val{1}, op::Call{typeof(sign), Tuple{TReal}}) = 0
grad!(dy::TAny, ::Val{1}, op::Call{typeof(tanh), Tuple{TReal}}) =
    (x = op.args[1]; t = tanh(x); (1.0f0 - t * t)  * dy)
grad!(dy::TAny, ::Val{1}, op::Call{typeof(sqrt), Tuple{TReal}}) = 0.5 * op.args[1] ^ (-0.5) * dy


text:
using NESSie.TestModel

@testset "bornion" begin
    for T in testtypes
        @test typeof(bornion("Na", T)) == BornIon{T}
    end
end

@test_skip φΩ(LocalES)
@test_skip φΣ(LocalES)
@test_skip φΩ(NonlocalES)
@test_skip φΣ(NonlocalES)


text:
# ------------------------------------------------------------------------------------------
# # Formatting The Plots
#
# Since the plotting functionality is implemented as a recipe to Plots.jl, [all of the
# options open to Plots.jl can be used in our
# plots](https://juliaplots.github.io/supported/). In addition, there are special features
# specifically for [differential equation
# plots](http://docs.juliadiffeq.org/latest/basics/plot.html). This tutorial will teach some
# of the most commonly used options. Let's first get the solution to some ODE. Here I will
# use one of the Lorenz ordinary differential equation. As with all commands in
# DifferentialEquations.jl, I got a plot of the solution by calling `solve` on the problem,
# and `plot` on the solution:
# ------------------------------------------------------------------------------------------

using DifferentialEquations, Plots
gr()
lorenz = @ode_def Lorenz begin
  dx = σ*(y-x)
  dy = ρ*x-y-x*z
  dz = x*y-β*z
end σ β ρ

p = [10.0,8/3,28]
u0 = [1., 5., 10.]
tspan = (0., 100.)
prob = ODEProblem(lorenz, u0, tspan, p)
sol = solve(prob)

plot(sol)

# ------------------------------------------------------------------------------------------
# Now let's change it to a phase plot. As discussed in the [plot functions
# page](http://docs.juliadiffeq.org/latest/basics/plot.html), we can use the `vars` command
# to choose the variables to plot. Let's plot variable `x` vs variable `y` vs variable `z`:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(:x,:y,:z))

# ------------------------------------------------------------------------------------------
# We can also choose to plot the timeseries for a single variable:
# ------------------------------------------------------------------------------------------

plot(sol,vars=[:x])

# ------------------------------------------------------------------------------------------
# Notice that we were able to use the variable names because we had defined the problem with
# the macro. But in general, we can use the indices. The previous plots would be:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3))
plot(sol,vars=[1])

# ------------------------------------------------------------------------------------------
# Common options are to add titles, axis, and labels. For example:
# ------------------------------------------------------------------------------------------

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
xaxis="Time (t)",yaxis="u(t) (in mm)",label=["X","Y","Z"])

# ------------------------------------------------------------------------------------------
# Notice that series recipes apply to the solution type as well. For example, we can use a
# scatter plot on the timeseries:
# ------------------------------------------------------------------------------------------

scatter(sol,vars=[:x])

# ------------------------------------------------------------------------------------------
# This shows that the recipe is using the interpolation to smooth the plot. It becomes
# abundantly clear when we turn it off using `denseplot=false`:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3),denseplot=false)

# ------------------------------------------------------------------------------------------
# When this is done, only the values the timestep hits are plotted. Using the interpolation
# usually results in a much nicer looking plot so it's recommended, and since the
# interpolations have similar orders to the numerical methods, their results are trustworthy
# on the full interval. We can control the number of points used in the interpolation's plot
# using the `plotdensity` command:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3),plotdensity=100)

# ------------------------------------------------------------------------------------------
# That's plotting the entire solution using 100 points spaced evenly in time.
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3),plotdensity=10000)

# ------------------------------------------------------------------------------------------
# That's more like it! By default it uses `100*length(sol)`, where the length is the number
# of internal steps it had to take. This heuristic usually does well, but unusually
# difficult equations it can be relaxed (since it will take small steps), and for equations
# with events / discontinuities raising the plot density can help resolve the discontinuity.
#
# Lastly notice that we can compose plots. Let's show where the 100 points are using a
# scatter plot:
# ------------------------------------------------------------------------------------------

plot(sol,vars=(1,2,3))
scatter!(sol,vars=(1,2,3),plotdensity=100)

# ------------------------------------------------------------------------------------------
# We can instead work with an explicit plot object. This form can be better for building a
# complex plot in a loop.
# ------------------------------------------------------------------------------------------

p = plot(sol,vars=(1,2,3))
scatter!(p,sol,vars=(1,2,3),plotdensity=100)
title!("I added a title")

# ------------------------------------------------------------------------------------------
# You can do all sorts of things. Have fun!
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
Overview: collection of subroutines for the Parametrisation
"""

"""
> identity_tangency!(Cp::Parametrisation,nm::Int64,neq::Int64,Φₗᵢₛₜ::Array{Int64},ϕ::Matrix{ComplexF64},ω₀::Array{Float64},ζ₀::Array{Float64})
It fills the first order parametrization using eigenmodes and eigenvalues of the mechanical problem.
- Cp : Parametrisation data structure
- nm : number of master modes
- neq : total number of equations
- Φₗᵢₛₜ : master modes list
- ϕ : modes
- ω₀ : eigenfrequencies
- ζ₀ : damping coefficients
- U : displacement field
"""
function identity_tangency!(Cp::Parametrisation,nm::Int64,neq::Int64,
                            Φₗᵢₛₜ::Array{Int64},ϕ::Matrix{Float64},
                            ω₀::Array{Float64},ζ₀::Array{Float64},
                            U::Field)
  for i = 1:nm
    λ₁ = - ζ₀[i]*ω₀[i] + ω₀[i]*sqrt(1.0-ζ₀[i]^2.0)*im
    λ₂ = - ζ₀[i]*ω₀[i] - ω₀[i]*sqrt(1.0-ζ₀[i]^2.0)*im
    Cp.f[i,i]       = λ₁
    Cp.f[i+nm,i+nm] = λ₂
    @inbounds for j = 1:U.neq
      Cp.W[j,       i]    = ϕ[j,Φₗᵢₛₜ[i]]*λ₁
      Cp.W[j+neq,   i]    = ϕ[j,Φₗᵢₛₜ[i]]
      Cp.W[j,    i+nm]    = ϕ[j,Φₗᵢₛₜ[i]]*λ₂
      Cp.W[j+neq,i+nm]    = ϕ[j,Φₗᵢₛₜ[i]]
    end
  end
  #
  return nothing
  #
end

"""
> realification!(Cp,Wr,fr,p,ndofs,neq,rmat)
It realifies mappings and reduced dynamics for an order p of the parametrisation
- Cp : Parametrisation data structure
- Wr : real-valued mapping
- fr : real-valued reduced dynamics
- p : order of the asymptotic development
- ndofs : number of degrees of freedom
- neq : number of equations
- rmat : realification matrix R
"""
function realification!(Cp,Wr,fr,p,ndofs,neq,rmat)
  #
  # It makes mappings and reduced dynamics real-valued
  #
  Cic = zeros(Int64,p)
  Ric = zeros(Int64,p)
  pm_comb = zeros(Int64,p)
  #
  if (p==0)
    @inbounds Wr[:] = Cp.W[:]
    @inbounds fr[:] = rmat*Cp.f[:]
    return nothing
  end
  #
  for i = 1:Cp.nc
    for j = 1:p
      Cic[j] = Cp.comb[j,i]
    end
    recursive_pm!(Cp,Wr,fr,Cic,Ric,pm_comb,p,ndofs,neq,i,1)
  end
  for i = 1:Cp.nc
    @inbounds fr[:,i] = rmat*fr[:,i]
  end
  return nothing
end

"""
> recursive_pm!(Cp,Wr,fr,Cic,Ric,pm_comb,p,ndofs,neq,comb,cc)
It realifies mappings and reduced dynamics of a given index combination comb of order p.
- Cp : parametrisation data structure
- Wr : real-valued mappings
- fr : real-valued reduced dynamics
- Cic : complex-valued parametrisation index combination of a given monomial
- Ric : real-valued parametrisation index combination of a given monomial
- pm_comb : +- combinations
- p : order of the asymptotic development
- ndofs : number of degrees of freedom
- neq : number of equations
- comb : reference index combination of the realified monomial
- cc : recursion depth
"""
function recursive_pm!(Cp,Wr,fr,Cic,Ric,pm_comb,p,ndofs,neq,comb,cc)
  #
  # Conversion of complex combinations to real-valued combinations.
  #
  nm = Int(ndofs/2)
  bin = [0, Int(ndofs/2)]
  for i = 1:2
    pm_comb[cc] = bin[i]
    if (cc<p)
      cc += 1
      recursive_pm!(Cp,Wr,fr,Cic,Ric,pm_comb,p,ndofs,neq,comb,cc)
      cc -= 1
    else
      pm = 1.0+0im
      #
      @inbounds for j = 1:p
        if (Cic[j]<=bin[2])
          Ric[j] = Cic[j]
        else
          Ric[j] = Cic[j] - bin[2]
        end
      end
      #
      @inbounds for j = 1:p
        if (pm_comb[j]!=0)
          pm *= 1.0/1.0im
          if (Cic[j]>bin[2])
            pm *= -1.0
          end
        end
      end
      #
      Ric = Ric + pm_comb
      pos = find_index_global(Ric,ndofs,p)
      cm = Cp.cmap[pos]
      cma = abs(cm)
      for j = 1:neq*2
        @inbounds Wr[j,cma] += (pm/2.0^p) * Cp.W[j,comb]
      end
      for j = 1:ndofs
        @inbounds fr[j,cma] += (pm/2.0^p) * Cp.f[j,comb]
      end
    end
  end
  return nothing
end

"""
> init_realification_matrix(ndofs::Int64)
It initializes the matrix used to realify the system
-
"""
function init_realification_matrix(ndofs::Int64)
  rmat = zeros(ComplexF64,(ndofs,ndofs))
  nm = Int(ndofs/2)
  @inbounds for i = 1:nm
    rmat[i,i] = 1.0
    rmat[i,i+nm] = 1.0
    rmat[i+nm,i] = 1.0im
    rmat[i+nm,i+nm] = -1.0im
  end
  return rmat
end


text:
using Flux
using NearestNeighbors
using StatsBase
using Zygote: @nograd

cd(@__DIR__)

include("xconv.jl")

struct PointCNN
    n_p::Int
    k::Int
    dilation::Int
    xconv::Xconv
end

function get_knn(P, n_p, dilation, k)
    p_idx = sample(1:size(P)[1], n_p; replace=false)
    ball_tree = [BallTree(transpose(P[:, 1, :, i])) for i in 1:size(P)[end]]
    return P[p_idx, :, :, :], [knn(ball_tree[i], transpose(P[p_idx, 1, :, i]), dilation * k)[1] for i in 1:size(P)[end]]
end

@nograd get_knn

PointCNN(n_p, cin, cout, cmid, k; dilation=1) = PointCNN(n_p, k, dilation, Xconv(cin, cout, 3, k, n_p, cmid, Int(floor(cout/(cin + cmid)))))

function (m::PointCNN)(P, F)
    p2, idxs = get_knn(P, m.n_p, m.dilation, m.k)
    p1 = cat([cat([reshape(P[idxs[i][j][1:m.dilation:(m.k * m.dilation)],1,:,i], (1, m.k, :)) for j in 1:m.n_p]..., dims=1) for i in 1:size(P)[end]]..., dims=4)
    if F != nothing
        F1 = cat([cat([reshape(F[idxs[i][j][1:m.dilation:(m.k * m.dilation)],1,:,i], (1, m.k, :)) for j in 1:m.n_p]..., dims=1) for i in 1:size(P)[end]]..., dims=4)
    else
        F1 = F
    end
    return p2, m.xconv(p2, p1, F1)
end

Flux.@functor PointCNN


text:
module ClusterAnalysis

using Statistics, LinearAlgebra
using Tables, NearestNeighbors

export KmeansResult, kmeans
export DBSCAN, dbscan

include("kmeans.jl")
include("DBSCAN.jl")

end


text:
using Measurements

# Figure 2A-D should be ordered:
#   ccw, lin, cw
# to be consistent with Figure 4.

# N-matrix, manually counted from bins:
# rows: cw, lin, ccw
# cols: in, out-50-bp, out-78-bp
N = [  63  54  48
      125 109  96
       61  55  48 ]

a = [ 0.22  -0.29  -0.545
      0.0   -0.31  -0.572
     -0.22  -0.35  -0.608 ]

σ = [ 0.05  0.03  0.07
      0.05  0.04  0.058
      0.05  0.03  0.059 ]

σ_mean = σ ./ sqrt.(N)

ā = a .± σ_mean

d(a_in, a_out) = (1 - a_in*a_in) / (1 - a_in*a_out) * a_out - a_in
d̃(a_in, a_out) = a_out - a_in

D = d.(ā[:, 1], ā[:, 2:3])
D̃ = d̃.(ā[:, 1], ā[:, 2:3])

Δ = D[2, :] .- (D[1, :] .+ D[3, :]) ./ 2
Δ̃ = D̃[2, :] .- (D̃[1, :] .+ D̃[3, :]) ./ 2


text:
# ---
# title: 809. Expressive Words
# id: problem809
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/expressive-words/description/>
# hidden: true
# ---
#
# Sometimes people repeat letters to represent extra feeling, such as "hello" ->
# "heeellooo", "hi" -> "hiiii".  In these strings like "heeellooo", we have
# _groups_ of adjacent letters that are all the same:  "h", "eee", "ll", "ooo".
#
# For some given string `S`, a query word is _stretchy_ if it can be made to be
# equal to `S` by any number of applications of the following _extension_
# operation: choose a group consisting of characters `c`, and add some number of
# characters `c` to the group so that the size of the group is 3 or more.
#
# For example, starting with "hello", we could do an extension on the group "o"
# to get "hellooo", but we cannot get "helloo" since the group "oo" has size
# less than 3.  Also, we could do another extension like "ll" -> "lllll" to get
# "helllllooo".  If `S = "helllllooo"`, then the query word "hello" would be
# stretchy because of these two extension operations: `query = "hello" ->
# "hellooo" -> "helllllooo" = S`.
#
# Given a list of query words, return the number of words that are stretchy.
#
#
#
#
#
#     **Example:**
#     Input:
#     S = "heeellooo"
#     words = ["hello", "hi", "helo"]
#     Output: 1
#     Explanation:
#     We can extend "e" and "o" in the word "hello" to get "heeellooo".
#     We can't extend "helo" to get "heeellooo" because the group "ll" is not size 3 or more.
#
#
#
#
# **Constraints:**
#
#   * `0 <= len(S) <= 100`.
#   * `0 <= len(words) <= 100`.
#   * `0 <= len(words[i]) <= 100`.
#   * `S` and all words in `words` consist only of lowercase letters
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


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

"""
    SimulationSolution

A solution to a spatial simulation problem.
"""
struct SimulationSolution{𝒟,ℛ}
  domain::𝒟
  realizations::ℛ
  nreals::Int

  function SimulationSolution{𝒟,ℛ}(domain, realizations) where {𝒟,ℛ}
    n = [length(r) for (var, r) in realizations]
    @assert length(unique(n)) == 1 "number of realizations must be unique"
    new(domain, realizations, n[1])
  end
end

SimulationSolution(domain::𝒟, realizations::ℛ) where {𝒟,ℛ} =
  SimulationSolution{𝒟,ℛ}(domain, realizations)

# -------------
# VARIABLE API
# -------------

Base.getindex(solution::SimulationSolution, var::Symbol) =
  solution.realizations[var]

# -------------
# ITERATOR API
# -------------

Base.iterate(solution::SimulationSolution, state=1) =
  state > solution.nreals ? nothing : (solution[state], state + 1)
Base.length(solution::SimulationSolution) = solution.nreals

# --------------
# INDEXABLE API
# --------------

function Base.getindex(solution::SimulationSolution, ind::Int)
  sdomain = solution.domain
  sreals  = solution.realizations
  idata   = DataFrame([var => reals[ind] for (var, reals) in sreals])
  georef(idata, sdomain)
end
Base.getindex(solution::SimulationSolution, inds::AbstractVector{Int}) =
  [getindex(solution, ind) for ind in inds]
Base.firstindex(solution::SimulationSolution) = 1
Base.lastindex(solution::SimulationSolution) = length(solution)

# ------------
# IO methods
# ------------
function Base.show(io::IO, solution::SimulationSolution)
  N = ncoords(solution.domain)
  print(io, "$(N)D SimulationSolution")
end

function Base.show(io::IO, ::MIME"text/plain", solution::SimulationSolution)
  println(io, solution)
  println(io, "  domain: ", solution.domain)
  println(io, "  variables: ", join(keys(solution.realizations), ", ", " and "))
  print(  io, "  N° reals:  ", solution.nreals)
end


text:
"""
    mutable struct Structure_VMD  <: Str

This contains information for `converter_vmd()` to assign all potential parameters.

"""
mutable struct Structure_VMD <: Str
    atom_name
    num_atom_types::Int64
    num_bond_types::Int64
    num_angle_types::Int64
    num_dihedral_types::Int64
    num_improper_types::Int64
    para_mass
    para_pair
    para_bond
    para_angle
    para_dihedral
    para_improper
end


text:
using Test

@testset "kpconv" begin
 @test kpatomequal(KPAtom(:f, Dict{Symbol,Any}(:a=>:b)), kpconv(:(f(a=b))))
 @test kpatomequal(KPAtom(:f, Dict{Symbol,Any}(:a=>:b,:c=>:(f(x)))), kpconv(:(f(c=f(x),a=b))))

 @test kpatomequal(KPAtom(:f, KParam()), kpconv(:(f())))
end

@testset "kptype" begin
 @test kplitequal(KPLiteral(:+, KPAtom(:P, Dict{Symbol,Any}())), lit2kplit(Meta.parse("+P()")))
 @test kplitequal(KPLiteral(:+, KPAtom(:P, Dict{Symbol,Any}(:b=>:(f(x))))), lit2kplit(Meta.parse("+P(b=f(x))")))
 @test kplitequal(KPLiteral(:+, KPAtom(:P, Dict{Symbol,Any}(:a=>:x,:b=>:(f(x))))), lit2kplit(Meta.parse("+P(a=x,b=f(x))")))

end

@testset "kpapply" begin
 @test kpatomequal(kpconv(:(好き(ha=僕, ga=君))), kpapply([:x], kpconv(:(好き(ha=僕, ga=君))), [:誰]))
 @test kpatomequal(kpconv(:(好き(ha=僕, ga=誰))), kpapply([:x], kpconv(:(好き(ha=僕, ga=x ))), [:誰]))

 @test kpatomequal(kpconv(:(好き(ha=僕, ga=f(誰)))), kpapply([:x], kpconv(:(好き(ha=僕, ga=x ))), [:(f(誰))]))
 @test kpatomequal(kpconv(:(好き(ha=僕(h(花子)), ga=f(誰)))), kpapply([:x,:y], kpconv(:(好き(ha=僕(y), ga=x ))), [:(f(誰)),:(h(花子))]))

end

@testset "kp unify" begin
 @test [] ==  kpunify([], kpconv(:(f())), kpconv(:(f())))
 @test [] ==  kpunify([], kpconv(:(f(x))), kpconv(:(f(x))))
 @test [:x] ==  kpunify([:x], kpconv(:(f(x))), kpconv(:(f(x))))
 @test [] ==  kpunify([], kpconv(:(f(a=x))), kpconv(:(f(a=x))))

 @test [:x] ==  kpunify([:x], kpconv(:(f(a=x))), kpconv(:(f(a=x))))
 @test [:k] ==  kpunify([:x], kpconv(:(f(a=x))), kpconv(:(f(a=k))))
 @test [:k] ==  kpunify([:x], kpconv(:(f(a=k))), kpconv(:(f(a=x))))
 @test_throws ICMP kpunify([:x], kpconv(:(f(a=k))), kpconv(:(f(a=p))))

 @test [:k] ==  kpunify([:x], kpconv(:(f(a=(f(x))))), kpconv(:(f(a=(f(k))))))

 @test [:(g(k))] ==  kpunify([:y], kpconv(:(f(a=y))), kpconv(:(f(a=g(k)))))
 @test [:(g(k))] ==  kpunify([:y], kpconv(:(f(a=(f(y))))), kpconv(:(f(a=f(g(k))))))
 @test [:(g(k)),:k] ==  kpunify([:y,:w], kpconv(:(f(a=f(k,y)))), kpconv(:(f(a=f(w,g(w))))))

 @test [] == kpunify([], kpconv(:(f(a=x,b=y))), kpconv(:(f(a=x,b=y))))
 @test [:x,:y] == kpunify([:x,:y], kpconv(:(f(a=x,b=y))), kpconv(:(f(a=x,b=y))))
 @test [:p,:q] == kpunify([:x,:y], kpconv(:(f(a=x,b=y))), kpconv(:(f(a=p,b=q))))

 @test_throws ICMP kpunify([:x], kpconv(:(f(a=x,b=x))), kpconv(:(f(a=p,b=q))))

 @test [:p,:(f(p))] == kpunify([:x,:y], kpconv(:(f(a=x,b=f(x)))), kpconv(:(f(a=p,b=y))))
 @test_throws ICMP kpunify([:x,:y], kpconv(:(f(a=x,b=g(x)))), kpconv(:(f(a=p,b=g(q)))))
 @test_throws ICMP kpunify([:x,:y], kpconv(:(f(a=x,b=g(x)))), kpconv(:(f(a=p,b=h(y)))))
 @test [:p,:p] == kpunify([:x,:y], kpconv(:(f(a=g(x),b=h(x)))), kpconv(:(f(a=g(p),b=h(y)))))
 @test [:(g(p)),:(g(p))] == kpunify([:x,:y], kpconv(:(f(a=g(x),b=h(x)))), kpconv(:(f(a=g(g(p)),b=h(y)))))
 @test [:p,:p] == kpunify([:x,:y], kpconv(:(f(a=x,b=x))), kpconv(:(f(a=y,b=p))))
 @test [:p,:p] == kpunify([:x,:y], kpconv(:(f(a=y,b=y))), kpconv(:(f(b=p, a=x))))

 @test [:p,:p] == kpunify([:x,:y], kpconv(:(f(a=x,b=x))), kpconv(:(f(b=p, a=y))))

 @test [:k, :k, :(h(k,k)), :(p(k))] == kpunify([:x,:y,:z,:w], kpconv(:(f(a=x, b=g(h(x,y),y), c=p(y)))), kpconv(:(f(a=k, b=g(z,k),c=w))))
 @test [:p,:p] == kpunify([:x,:y], kpconv(:(f(a=x,b=x))), kpconv(:(f(a=p,b=y))))
 @test [:k, :k, :p] == kpunify([:x,:y,:z], kpconv(:(f(a=x, b=g(p,y)))), kpconv(:(f(a=k, b=g(z,k)))))
 @test [:k, :k, :(h(k,k))] == kpunify([:x,:y,:z], kpconv(:(f(a=x, b=g(h(x,y),y)))), kpconv(:(f(a=k, b=g(z,k)))))
 @test [:x, :k, :(h(x,k))] == kpunify([:x,:y,:z], kpconv(:(f(d=x, b=g(h(x,y),y)))), kpconv(:(f(a=k, b=g(z,k),c=:(g(x))))))

# common key args only unified
 @test [:k, :k, :(h(k,k))] == kpunify([:x,:y,:z], kpconv(:(f(a=x, b=g(h(x,y),y)))), kpconv(:(f(a=k, b=g(z,k),c=:(g(x))))))

println()
println("HERE")

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03md(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    NMAX = 20
    Q = Array{Float64,1}(undef, NMAX)
    E = Array{Float64,1}(undef, NMAX-1)
    Q2 = Array{Float64,1}(undef, NMAX)
    E2 = Array{Float64,1}(undef, NMAX-1)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    THETA = parse(Float64, replace(vs[2],'D'=>'E'))
    L = parse(BlasInt, vs[3])
    TOL = parse(Float64, replace(vs[4],'D'=>'E'))
    RELTOL = parse(Float64, replace(vs[5],'D'=>'E'))
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end
    if ( L<0 || L>N )
        @error "Illegal L=$L"
    end

    vs = String[]
    _isz = N
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    Q[1:_isz] .= parsex.(Float64, vs)

    vs = String[]
    _isz = N-1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    E[1:_isz] .= parsex.(Float64, vs)
    close(f)

    println(io, "J:")
    show(io, "text/plain", Bidiagonal(Q[1:N],E[1:N-1],:U))
    println(io)
    Q2[N] = Q[N]^2
    PIVMIN = Q2[N]
    for i in 1:N-1
        Q2[i] = Q[i]^2
        E2[i] = E[i]^2
        PIVMIN = max( PIVMIN, Q2[i], E2[i] )
    end
    SAFMIN = floatmin(1.0)
    PIVMIN = max( PIVMIN*SAFMIN, SAFMIN )
    TOL = max( TOL, ZERO )
    if RELTOL <= 0
        RELTOL = 2.0 * eps(1.0)
    end
# interp call 1

    L, THETA, INFO, IWARN = SLICOT.mb03md!(N, L, THETA, Q, E, Q2, E2, PIVMIN, TOL, RELTOL)
    @test INFO == 0
    println(io, "L = $L")
    println(io, "THETA = $THETA")
    println(io, "IWARN = $IWARN")

end # run_mb03md()


text:
# Unicode iteration and decoding for StringView, copied from the corresponding String functions in Base

@inline between(b::T, lo::T, hi::T) where {T<:Integer} = (lo ≤ b) & (b ≤ hi)

@inline function Base.iterate(s::StringView, i::Int=firstindex(s))
    (i % UInt) - 1 < ncodeunits(s) || return nothing
    b = @inbounds codeunit(s, i)
    u = UInt32(b) << 24
    between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return iterate_continued(s, i, u)
end

function iterate_continued(s::StringView, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

Base.@propagate_inbounds function Base.getindex(s::StringView, i::Int)
    b = codeunit(s, i)
    u = UInt32(b) << 24
    between(b, 0x80, 0xf7) || return reinterpret(Char, u)
    return getindex_continued(s, i, u)
end

function getindex_continued(s::StringView, i::Int, u::UInt32)
    if u < 0xc0000000
        # called from `getindex` which checks bounds
        @inbounds isvalid(s, i) && @goto ret
        Base.string_index_err(s, i)
    end
    n = ncodeunits(s)

    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 1
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16

    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 2
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8

    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 3
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b)
@label ret
    return reinterpret(Char, u)
end

Base.getindex(s::StringView, r::UnitRange{<:Integer}) = s[Int(first(r)):Int(last(r))]

@inline function Base.getindex(s::StringView, r::UnitRange{Int})
    isempty(r) && return ""
    i, j = first(r), last(r)
    @boundscheck begin
        checkbounds(s, r)
        @inbounds isvalid(s, i) || string_index_err(s, i)
        @inbounds isvalid(s, j) || string_index_err(s, j)
    end
    j = nextind(s, j) - 1
    return StringView(@view s.data[i:j])
end

Base.length(s::StringView) = length_continued(s, 1, ncodeunits(s), ncodeunits(s))

@inline function Base.length(s::StringView, i::Int, j::Int)
    @boundscheck begin
        0 < i ≤ ncodeunits(s)+1 || throw(BoundsError(s, i))
        0 ≤ j < ncodeunits(s)+1 || throw(BoundsError(s, j))
    end
    j < i && return 0
    @inbounds i, k = thisind(s, i), i
    c = j - i + (i == k)
    length_continued(s, i, j, c)
end

@inline function length_continued(s::StringView, i::Int, n::Int, c::Int)
    i < n || return c
    @inbounds b = codeunit(s, i)
    @inbounds while true
        while true
            (i += 1) ≤ n || return c
            0xc0 ≤ b ≤ 0xf7 && break
            b = codeunit(s, i)
        end
        l = b
        b = codeunit(s, i) # cont byte 1
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xe0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 2
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xf0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 3
        c -= (b & 0xc0 == 0x80)
    end
end


text:
using OodleAnalytics
using Test

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


text:
using LinearAlgebra
include("../Utils/MeshGeometryUtils.jl")

function getnodesgeometricalquantities(nodes, faces)
end

function getfacesgeometricalquantities(nodes, faces)
    for face in faces
        inodes = face.inodes

        # Collect the nodes of the face
        facepoints = []
        for inode in inodes
            push!(facepoints, nodes[inode].centroid)
        end

        # Compute the geometric center of the face
        facegc = computegeometriccenter(facepoints)

        # Start computing face centroid, surface normal, and area --------------
        facecentroid = [0, 0, 0]
        facesf = [0, 0, 0]
        facearea = 0

        # Loop over the entire sub-triangles (local faces)
        for i in 1:length(facepoints)
            r1 = facepoints[i]
            if i < length(facepoints)
                r2 = facepoints[i + 1]
            else
                r2 = facepoints[1]
            end

            # Collect points that make up the local face
            localfacepoints = [r1, r2, facegc]

            # Compute the local face centroid
            localfacecentroid = computegeometriccenter(localfacepoints)

            # Compute the local face surface normal and area
            localfacesf = computesurfacenormal(localfacepoints)
            localfacearea = computemagnitude(localfacesf)

            # Update the face properties
            facecentroid += localfacearea*localfacecentroid
            facesf += localfacesf
            facearea += localfacearea
        end

        # Store face properties
        face.centroid = facecentroid/facearea
        face.sf = facesf
        face.area = facearea

        # End computing total face centorid, surface normal, and area ----------
    end
end

function getcellsgeometricalquantities(nodes, faces, cells)
    for cell in cells
        inodescell = cell.inodes

        # Collect the nodes of the cell
        cellpoints = []
        for inode in inodescell
            push!(cellpoints, nodes[inode].centroid)
        end

        # Compute the geometric center of the cell
        cellgc = computegeometriccenter(cellpoints)

        # Start computing cell centroid and volume -----------------------------
        cellcentroid = [0, 0, 0]
        cellvolume = 0

        # Loop over the entire faces of the cell (local cells)
        for iface in cell.ifaces
            # Get the face geometric quantities that have been calculated
            facearea = faces[iface].area
            facecentroid = faces[iface].centroid
            faceinodes = faces[iface].inodes

            # Collect the nodes of the face
            facepoints = []
            for inode in faceinodes
                push!(facepoints, nodes[inode].centroid)
            end
            push!(facepoints, cellgc)

            # Compute the geometric center of the local cell (pyramid)
            localcellgc = computegeometriccenter(facepoints)

            # Compute the centroid of the local cell
            localcellcentroid = 0.75*facecentroid + 0.25*localcellgc

            # Compute the volume of the local cell
            localcellbase = facearea
            localcellheight = computemagnitude(cellgc - facecentroid)
            localcellvolume = (1/3)*localcellbase*localcellheight

            # Update the cell properties
            cellcentroid += localcellvolume*localcellcentroid
            cellvolume += localcellvolume
        end
        # Store the computed cell properties
        cell.centroid = cellcentroid / cellvolume
        cell.volume = cellvolume
    end
end

function getfacesremainingproperties(nodes, faces, cells)
    for face in faces
        iowner = face.iowner
        ineighbour = face.ineighbour
        sf = face.sf
        ef = sf/(computemagnitude(sf))
        dof = face.centroid - cells[iowner].centroid

        if ineighbour != -1
            # Compute geometric factor (gf)
            dfn = cells[ineighbour].centroid - face.centroid
            gf = dot(dof, ef) / (dot(dof, ef) + dot(dfn, ef))

            # Compute distance vector T joining owner and neighbour
            T = dfn + dof

            # Compute wall distance
            walldist = 0
        else
            # Compute geometric factor (gf)
            gf = 1

            # Compute distance vector T joining owner and neighbour
            T = dof

            # Compute wall distance
            walldist = dot(dof, ef)
        end

        # Store the face remaining properties
        face.gf = gf
        face.cn = dof
        face.T = T
        face.walldist = walldist

    end
end


text:
# TODO parameterise time type?

# TODO Using an abstractvector for values requires some thought. The problem is whether or
#   not it would permit us to have views of an array.
#   Possibly we need a more general type here which can package other representations?
#   See Stheno ColVecs / RowVecs - they have the same problem.
#
#   Also see: https://juliaarrays.github.io/ArraysOfArrays.jl/stable/

function _is_strictly_increasing(x::AbstractVector)
    length(x) < 2 && return true

    previous_value = @inbounds first(x)
    for value in @inbounds(x[2:end])
        if previous_value >= value
            return false
        end
        previous_value = value
    end
    return true
end

"""
Sentinel type for use in the constructor for `Block` below.
"""
struct UncheckedConstruction end
const unchecked = UncheckedConstruction()

"""
    Block{T}()
    Block(times::AbstractVector{DateTime}, values::AbstractVector{T})
    Block(unchecked, times, values)

Represent some data in timeseries.

Conceptually this is a list of `(time, value)` pairs, or "knots". Times must be strictly
increasing — i.e. no repeated timestamps are allowed.

The constructor `Block(times, values)` will verify that the input data satisfies this
constraint, however `Block(unchecked, times, values)` will skip the checks. This is
primarily intended for internal use, where the caller assumes responsibility for the
validity of `times` & `values`.

!!! danger
    `TimeDag` considers instances of `Block` to be completely immutable. Thus, when working
    with functions that accept blocks (e.g. [`TimeDag.run_node!`](@ref)), you must not
    modify `times` or `values` members.
"""
struct Block{T,VTimes<:AbstractVector{DateTime},VValues<:AbstractVector{T}}
    times::VTimes
    values::VValues

    function Block(
        ::UncheckedConstruction, times::VTimes, values::VValues
    ) where {T,VTimes<:AbstractVector{DateTime},VValues<:AbstractVector{T}}
        return new{T,VTimes,VValues}(times, values)
    end

    function Block(
        times::VTimes, values::VValues
    ) where {T,VTimes<:AbstractVector{DateTime},VValues<:AbstractVector{T}}
        if length(times) != length(values)
            throw(
                ArgumentError(
                    "Times have length $(length(times)), values $(length(values))"
                ),
            )
        end

        if !_is_strictly_increasing(times)
            throw(ArgumentError("Times are not strictly increasing."))
        end

        # TODO
        #   -> How make times & values immutable? Subclass of AbstractVector that doesn't
        #       implement setindex! ?
        return new{T,VTimes,VValues}(times, values)
    end
end

Block{T}() where {T} = Block(unchecked, DateTime[], T[])

function Block(
    knots::Union{AbstractVector{Tuple{DateTime,T}},AbstractVector{Pair{DateTime,T}}}
) where {T}
    n = length(knots)
    times = _allocate_times(n)
    values = _allocate_values(T, n)
    for (i, (time, value)) in enumerate(knots)
        @inbounds times[i] = time
        @inbounds values[i] = value
    end
    return Block(times, values)
end

value_type(::Block{T}) where {T} = T

# Blocks are considered immutable from the point of view of `duplicate`.
duplicate_internal(x::Block, ::IdDict) = x

# Equality & hash need to be defined, since we have defined an internal constructor.
Base.hash(a::Block, h::UInt) = hash(a.values, hash(a.times, hash(:Block, h)))
function Base.isequal(a::Block, b::Block)
    a === b && return true
    return isequal(a.times, b.times) && isequal(a.values, b.values)
end
function Base.:(==)(a::Block, b::Block)
    a === b && return true
    return a.times == b.times && a.values == b.values
end

# Implement approximate equality in terms of exact equality for timestamps, but approximate
# for values. Only blocks with the same value type should compare approximately equal.
function Base.isapprox(a::Block{T}, b::Block{T}; kwargs...) where {T}
    isempty(a) && isempty(b) && return true
    a.times == b.times || return false
    return isapprox(a.values, b.values; kwargs...)
end

# Make a block behave like a table with two columns, primarily for printing purposes.
Tables.istable(::Type{<:Block}) = true
Tables.columnaccess(::Type{<:Block}) = true
Tables.columns(block::Block) = block

const _BLOCK_COLUMNNAMES = (:time, :value)
Tables.schema(::Block{T}) where {T} = Tables.Schema(_BLOCK_COLUMNNAMES, (DateTime, T))
Tables.getcolumn(block::Block, i::Int) = Tables.getcolumn(block, _BLOCK_COLUMNNAMES[i])
function Tables.getcolumn(block::Block, nm::Symbol)
    if nm == :time
        return block.times
    elseif nm == :value
        return block.values
    else
        throw(ArgumentError("Unknown column $nm"))
    end
end
Tables.columnnames(::Block) = _BLOCK_COLUMNNAMES

function Base.show(io::IO, block::Block)
    return pretty_table(
        io,
        block;
        title="Block{$(value_type(block))}($(length(block)) knots)",
        tf=tf_markdown,
    )
end

"""Slice the block so we contain values in the interval [time_start, time_end)."""
function _slice(block::Block{T}, time_start::DateTime, time_end::DateTime) where {T}
    i_first = searchsortedfirst(block.times, time_start)
    i_last = searchsortedfirst(block.times, time_end) - 1

    if (i_first > length(block) || i_last < 1)
        return Block{T}()
    elseif (i_first == 1 && i_last == length(block))
        # Avoid constructing views if we want the whole block.
        return block
    else
        return Block(
            unchecked,
            @inbounds(@view(block.times[i_first:i_last])),
            @inbounds(@view(block.values[i_first:i_last])),
        )
    end
end

function Base.vcat(blocks::Block{T}...) where {T}
    return if length(blocks) == 1
        # Fast path for a single block.
        only(blocks)
    else
        # TODO Rather than using vcat here, we could use a function that more intelligently
        # merges views and ranges etc.
        # TODO we could use unchecked here, if we were to check that all the adjacent
        #   elements between blocks were correct.
        Block(
            vcat((block.times for block in blocks)...),
            vcat((block.values for block in blocks)...),
        )
    end
end

# Indexing, iteration, etc.
# NB: constraining length to return an Int here is very important for performance!
Base.length(block::Block)::Int = length(block.times)
Base.isempty(block::Block)::Bool = length(block) == 0
Base.firstindex(block::Block) = 1
Base.lastindex(block::Block) = length(block)
Base.getindex(block::Block, i::Int) = (block.times[i], block.values[i])
Base.getindex(block::Block, i::OrdinalRange{Int}) = Block(block.times[i], block.values[i])

Base.eltype(::Block{T}) where {T} = Tuple{DateTime,T}
function Base.iterate(block::Block, state::Int=1)
    return if (state > length(block))
        nothing
    else
        (block[state], state + 1)
    end
end

_allocate_times(n::Int) = Vector{DateTime}(undef, n)

"""
    _allocate_values(T, n::Int) -> AbstractVector{T}

Allocate some uninitialized memory that looks like a vector{T} of length `n`.
"""
function _allocate_values(T, n::Int)
    # TODO This is not necessarily optimal. See comment near definition of
    #   Block regarding what to do if we're viewing a dense array.
    #   The correct solution will probably involve dispatching based on the type of
    #   input.values.
    return Vector{T}(undef, n)
end

"""Resize the vector to length `n` and attempt to reclaim unused space."""
function _trim!(x::AbstractVector, n::Int)
    # TODO It sounds like this currently doesn't actually free any of the buffer, which
    #   could be a bit inefficient. Maybe sizehint! is required too?
    return resize!(x, n)
end

"""
    _equal_times(a::Block, b::Block) -> Bool

true => the times in blocks `a` and `b` are the same.
false => the times in blocks `a` and `b` may be different, or the same.

This function will try to return true for as many cases as possible, with the guarantee
that it will always run in constant time; i.e. it will never explicitly compare time values.
"""
function _equal_times(a::Block, b::Block)::Bool
    # FIXME This doesn't account for the case where e.g. a.times is a vector, and b.times is
    #   a view of the entirety of a.times.
    return a.times === b.times
end


text:
using LinearAlgebra, SparseArrays, Preconditioners, IterativeSolvers, Random

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

Random.seed!(1)

function test_matrix(A, F, atol)
    n = size(A, 1)
    if F === CholeskyPreconditioner
        C1 = EmptyCholeskyPreconditioner(A)
        UpdatePreconditioner!(C1, A, 2)

        C2 = CholeskyPreconditioner(A, 2)
        @test isapprox(C1.L, C2.L, atol=atol)

        C3 = CholeskyPreconditioner(Symmetric(A), 2)
        @test isapprox(C2.L, C3.L, atol=atol)

        C4 = CholeskyPreconditioner(A, n)
        @test isapprox(norm(C4 \ ones(n) - Symmetric(A) \ ones(n), Inf), 0.0, atol=0.001)
    end
    if F === RugeStuben || F === SmoothedAggregation
        p = AMGPreconditioner(F, A)
    else
        p = F(A)
    end
    @test isapprox(cg(A, A*ones(n), Pl=p), ones(n), atol=atol)
    if F === RugeStuben
        p = AMGPreconditioner(F, A)
        @test isapprox(cg(A, A*ones(n), Pl=p), ones(n), atol=atol)
    end
    A = sprand(n, n, 10/n)
    A = A + A' + 30I
    UpdatePreconditioner!(p, A)
    @test isapprox(cg(A, A*ones(n), Pl=p), ones(n), atol=atol)
end

@testset "$T preconditioner" for (T, F) in (("Diagonal", DiagonalPreconditioner), ("Incomplete Cholesky", CholeskyPreconditioner), ("AMG Ruge-Stuben", RugeStuben), ("AMG Smoothed Aggregation", SmoothedAggregation))
    n = 100
    atol = sqrt(eps(Float64))*n
    A = sprand(n, n, 10/n)
    A = A + A' + 50I
    test_matrix(A, F, atol)
    test_matrix(Symmetric(A), F, atol)
    test_matrix(Hermitian(A), F, atol)
end


text:
using Base.Test
using ArrayViewsAPL

A = reshape(float64(1:300), 5, 5, 3, 4)
S = subview(A, 1:2, 1:2, 1:2, 1:2)
@test isa(S, View{Float64, 4})
@test size(S) == (2,2,2,2)
@test strides(S) == (1,5,25,75)
S = sliceview(A, 1:2, 1:2, 1:2, 1:2)
@test isa(S, View{Float64, 4})
@test size(S) == (2,2,2,2)
@test strides(S) == (1,5,25,75)
S = subview(A, 1:2, 1:2, 1:2, 2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},UnitRange{Int},Int)})
@test size(S) == (2,2,2)
@test strides(S) == (1,5,25)
S = sliceview(A, 1:2, 1:2, 1:2, 2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},UnitRange{Int},Int)})
@test size(S) == (2,2,2)
@test strides(S) == (1,5,25)
S = subview(A, 1:2, 1:2, 2, 1:2)
@test isa(S, View{Float64, 4, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},UnitRange{Int},UnitRange{Int})})
@test size(S) == (2,2,1,2)
@test strides(S) == (1,5,25,75)
S = sliceview(A, 1:2, 1:2, 2, 1:2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,UnitRange{Int})})
@test size(S) == (2,2,2)
@test strides(S) == (1,5,75)
P = S
S = subview(P, 1:2, 1:2, 1:2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,UnitRange{Int})})
@test size(S) == (2,2,2)
@test strides(S) == (1,5,75)
S = sliceview(P, 1:2, 1:2, 1:2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,UnitRange{Int})})
@test size(S) == (2,2,2)
@test strides(S) == (1,5,75)
S = subview(P, 1:2, 1:2, 2)
@test isa(S, View{Float64, 2, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,Int)})
@test size(S) == (2,2)
@test strides(S) == (1,5)
S = sliceview(P, 1:2, 1:2, 2)
@test isa(S, View{Float64, 2, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,Int)})
@test size(S) == (2,2)
@test strides(S) == (1,5)
S = subview(P, 1:2, 2, 1:2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,UnitRange{Int})})
@test size(S) == (2,1,2)
@test strides(S) == (1,5,75)
S = sliceview(P, 1:2, 2, 1:2)
@test isa(S, View{Float64, 2, Array{Float64,4}, (UnitRange{Int},Int,Int,UnitRange{Int})})
@test size(S) == (2,2)
@test strides(S) == (1,75)
S = subview(P, 2, 1:2, 1:2)
@test isa(S, View{Float64, 3, Array{Float64,4}, (UnitRange{Int},UnitRange{Int},Int,UnitRange{Int})})
@test size(S) == (1,2,2)
@test strides(S) == (1,5,75)
S = sliceview(P, 2, 1:2, 1:2)
@test isa(S, View{Float64, 2, Array{Float64,4}, (Int,UnitRange{Int},Int,UnitRange{Int})})
@test size(S) == (2,2)
@test strides(S) == (5,75)

## 2d
# Whole-array views
A = reshape(1:24, 6, 4)
@test typeof(A[1:2,1:2]) == Array{Int, 2}  # so we won't be fooled when this changes
B = subview(A, 1:6, 1:4)
@test B[2] == A[2]
@test B[3,1] == A[3,1]
@test B[4,2,1] == A[4,2]
@test B == A
@test strides(B) == (1,6)
@test pointer(B) == pointer(A)
B = sliceview(A, 1:6, 1:4)
@test B == A
@test B[2] == A[2]
@test B[3,1] == A[3,1]
@test B[4,2,1] == A[4,2]
@test pointer(B) == pointer(A)
B = subview(A, :, :)
@test B == A
B = sliceview(A, :, :)
@test B == A
# Views with extra dimensions
B = subview(A, 1:6, 1:4, 1)
@test ndims(B) == 2
@test size(B) == (6,4)
@test strides(B) == (1,6)
@test B[7] == A[7]
@test B[6,4] == A[6,4]
@test B[3,3,1] == A[3,3]
@test B[4,2,1,1] == A[4,2]
@test B == A
B = sliceview(A, 1:6, 1:4, 1)
@test ndims(B) == 2
@test size(B) == (6,4)
@test strides(B) == (1,6)
@test B[7] == A[7]
@test B[6,4] == A[6,4]
@test B[3,3,1] == A[3,3]
@test B[4,2,1,1] == A[4,2]
@test B == A
B = subview(A, 1:6, 1:4, 2)
@test_throws BoundsError B[1,1]
@test_throws BoundsError B[1]
B = sliceview(A, 1:6, 1:4, 2)
@test_throws BoundsError B[1,1]
@test_throws BoundsError B[1]
# Partial views
B = subview(A, 2:2:6, 2)
@test ndims(B) == 1
@test strides(B) == (2,)
@test pointer(B) == pointer(A)+7*sizeof(Int)
@test B == [8:2:12]
@test B[1] == 8
@test B[1,1] == 8
@test B[3] == 12
@test B[3,1] == 12
B = sliceview(A, 2:2:6, 2)
@test ndims(B) == 1
@test pointer(B) == pointer(A)+7*sizeof(Int)
@test B == [8:2:12]
@test B[1] == 8
@test B[1,1] == 8
@test B[3] == 12
@test B[3,1] == 12
B = subview(A, 2, 2:4)
@test ndims(B) == 2
@test strides(B) == (1,6)
@test pointer(B) == pointer(A)+7*sizeof(Int)
@test B == A[2,2:4]
@test B[1] == 8
@test B[1,1] == 8
@test B[1,1,1] == 8
@test B[3] == 20
@test B[1,3] == 20
@test B[1,3,1] == 20
B = sliceview(A, 2, 2:4)
@test ndims(B) == 1
@test strides(B) == (6,)
@test pointer(B) == pointer(A)+7*sizeof(Int)
@test B == squeeze(A[2,2:4], 1)
@test B[1] == 8
@test B[1,1] == 8
@test B[1,1,1] == 8
@test B[3] == 20
@test B[3,1] == 20
# Views created with fewer dimensions (linear indexing)
B = subview(A, 5:13)
@test ndims(B) == 1
@test strides(B) == (1,)
@test pointer(B) == pointer(A) + 4*sizeof(Int)
@test B == [5:13]
@test B[1] == 5
@test B[9] == 13
@test B[2,1] == 6
@test B[8,1] == 12
B = sliceview(A, 5:13)
@test ndims(B) == 1
@test pointer(B) == pointer(A) + 4*sizeof(Int)
@test B == [5:13]
@test B[1] == 5
@test B[9] == 13
@test B[2,1] == 6
@test B[8,1] == 12
B = subview(A, 1:3:6, 2:2:4)
@test pointer(B) == pointer(A) + 6*sizeof(Int)
@test size(B) == (2,2)
@test strides(B) == (3,12)
@test B[1] == 7
@test B[3] == 19
@test B[2,1] == 10
@test B[2,2] == 22
@test B[1,2,1] == 19
# Flipped axes
B = subview(A, 5:-1:3, 1:4)
@test size(B) == (3,4)
@test strides(B) == (-1,6)
@test pointer(B) == pointer(A) + 4*sizeof(Int)
@test B[1,1] == A[5,1]
@test B[2] == A[4,1]
@test B[2,3,1] == A[4,3]
B = subview(A, 5:-1:3, 4:-2:1)
@test size(B) == (3,2)
@test strides(B) == (-1,-12)
@test pointer(B) == pointer(A) + 22*sizeof(Int)
@test B[1] == 23
@test B[5] == 10
@test B[3,1] == 21
@test B[2,2] == 10
@test B[3,2,1] == 9
# Views with Vector{Int}
B = subview(A, [1,3,4], 1:2:4)
@test B[2,2] == A[3,3]
@test_throws ErrorException strides(B)
@test_throws MethodError pointer(B)

I1 = 2:2:8
I2 = 5:3:12
dims = (13,)   # could supply a 2nd size, but it's not used so why bother
I12 = [54,56,58,60,93,95,97,99,132,134,136,138]
L = length(I12)
ind = ArrayViewsAPL.merge_indexes((I1, I2), dims, 1:L)
@test ind == I12
ind = ArrayViewsAPL.merge_indexes((I1, I2), dims, 3:7)
@test ind == I12[3:7]
ind = ArrayViewsAPL.merge_indexes_div((I1, I2), dims, 1:L)
@test ind == I12
ind = ArrayViewsAPL.merge_indexes_div((I1, I2), dims, 3:7)
@test ind == I12[3:7]

## 2d, views-of-views
# Whole-array views
B = subview(A, :, :)
C = subview(B, :, :)
@test C == A
C = sliceview(B, :, :)
@test C == A
# with extra indexes
C = subview(B, 1:6, 1:4, 1)
@test ndims(C) == 2
@test size(C) == (6,4)
@test strides(C) == (1,6)
@test pointer(C) == pointer(A)
@test C[2] == A[2]
@test C[3,2] == A[3,2]
@test C[4,3,1] == A[4,3]
C = sliceview(B, 1:6, 1:4, 1)
@test ndims(C) == 2
@test size(C) == (6,4)
@test strides(C) == (1,6)
@test pointer(C) == pointer(A)
@test C[2] == A[2]
@test C[3,2] == A[3,2]
@test C[4,3,1] == A[4,3]
C = subview(B, 1:6, 1:4, 2)
@test_throws BoundsError C[1]
C = subview(B, 1:6, 1:4, 1:1)
@test ndims(C) == 3
@test size(C) == (6,4,1)
@test strides(C) == (1,6,24)
@test pointer(C) == pointer(A)
@test C[2] == A[2]
@test C[3,2] == A[3,2]
@test C[4,3,1] == A[4,3]
C = sliceview(B, 1:6, 1:4, 1:1)
@test ndims(C) == 3
@test size(C) == (6,4,1)
@test strides(C) == (1,6,24)
@test pointer(C) == pointer(A)
@test C[2] == A[2]
@test C[3,2] == A[3,2]
@test C[4,3,1] == A[4,3]

# Views of partial views
# For this, we need a bigger A
A = reshape(1:13*14, 13, 14)
I1, I2 = 3:2:13, 4:3:11
B = subview(A, I1, I2)
As = A[I1, I2]
@test B == As
C = subview(B, 1:6, 1:3)
@test C == As
@test strides(C) == (2,13*3)
C = sliceview(B, 1:6, 1:3)
@test C == As
C = subview(B, 2:2:5, 2:3)
@test C == As[2:2:5, 2:3]
@test strides(C) == (4,13*3)
C = sliceview(B, 2:2:5, 2:3)
@test C == As[2:2:5, 2:3]
C = subview(B, 3:5, 1:2:3)
@test C == As[3:5, 1:2:3]
@test strides(C) == (2,13*6)
# With fewer indexes (linear indexing)
C = subview(B, 2:2:8)
@test ndims(C) == 1
@test size(C) == (4,)
@test_throws ErrorException strides(C)
@test C[1] == As[2]
@test C[2,1] == As[4]
C = sliceview(B, 2:2:8)
@test ndims(C) == 1
@test size(C) == (4,)
@test_throws ErrorException strides(C)
@test C[1] == As[2]
@test C[2,1] == As[4]

# with extra indexes
C = sliceview(B, 1:6, 3, 1:1)
@test ndims(C) == 2
@test size(C) == (6,1)
@test strides(C) == (2,13*14)
@test C[2] == As[2,3]
@test C[4,1] == As[4,3]
@test C[3,1,1] == As[3,3]
@test unsafe_load(pointer(C)) == As[1,3]


text:
export KProbe, USDT, Tracepoint, PerfEvent, XDP

using Libdl
using ObjectFile

abstract type AbstractProbe end
objects(p::AbstractProbe) = [p.obj]

"Merges maps from `other` into `obj`."
function merge_maps!(obj, other::AbstractProbe)
    for other_obj in objects(other)
        merge_maps!(obj, other_obj)
    end
end
function merge_maps!(obj, other_obj::API.Object)
    for other_map in API.maps(other_obj)
        for our_map in API.maps(obj)
            if API.name(our_map) == API.name(other_map)
                API.reuse_fd(our_map, API.fd(other_map))
            end
        end
    end
end

struct ProbeSet <: AbstractProbe
    probes::Vector{AbstractProbe}
end
ProbeSet() = ProbeSet(AbstractProbe[])
objects(p::ProbeSet) = vcat(map(objects, p.probes)...)
function API.findmap(p::ProbeSet, name::String)
    for obj in objects(p)
        try
            return API.findmap(obj, name)
        catch err
            err isa AssertionError || rethrow(err)
        end
    end
    throw(ArgumentError("Failed to find map $name"))
end

struct KProbe <: AbstractProbe
    obj::API.Object
    kfunc::String
    retprobe::Bool
end
function KProbe(f::Function, ::Type{T}, kfunc; merge_with=(), retprobe=false, kwargs...) where {T<:Tuple}
    obj = API.Object(bpffunction(f, T; kwargs...))
    foreach(prog->API.set_kprobe!(prog), API.programs(obj))
    for other in merge_with
        merge_maps!(obj, other)
    end
    KProbe(obj, kfunc, retprobe)
end
KProbe(f::Function, kfunc; kwargs...) =
    KProbe(f, Tuple{API.pointertype(API.pt_regs)}, kfunc; kwargs...)

const NOTE_ARG_TYPE = @NamedTuple{size::Int,
                                  signed::Bool,
                                  reg::Union{Symbol,Nothing},
                                  offset::Union{Int,Nothing}}
const NOTE_TYPE = @NamedTuple{owner::String,
                              location::UInt64,
                              base::UInt64,
                              semaphore::UInt64,
                              provider::String,
                              func::String,
                              args::Vector{NOTE_ARG_TYPE}}

struct UProbe#={F<:Function,T}=# <: AbstractProbe
    obj::API.Object
    pid::UInt32
    binpath::String
    addr::UInt64
    #=
    func::F
    sig::T
    =#
    retprobe::Bool
end
#=
function UProbe(f::Function, ::Type{T}, method, sig; merge_with=(), retprobe=false, kwargs...) where {T<:Tuple}
    obj = API.Object(bpffunction(f, T; kwargs...))
    #foreach(prog->API.set_uprobe!(prog), API.programs(obj))
    foreach(prog->API.set_kprobe!(prog), API.programs(obj))
    for other in merge_with
        merge_maps!(obj, other)
    end
    UProbe(obj, method, sig, retprobe)
end
UProbe(f::Function, method, sig; kwargs...) =
    UProbe(f, Tuple{API.pointertype(API.pt_regs)}, method, sig; kwargs...)
function UProbe(f::Function, binpath, addr; pid=0, retprobe=false, kwargs...)
    obj = API.Object(bpffunction(f, Tuple{API.pointertype(API.pt_regs)}; kwargs...))
    #foreach(prog->API.set_uprobe!(prog), API.programs(obj))
    foreach(prog->API.set_kprobe!(prog), API.programs(obj))
    UProbe(obj, pid, binpath, addr, retprobe)
end
=#
struct USDT <: AbstractProbe
    obj::API.Object
    pid::UInt32
    binpath::String
    func::String
    addr::UInt64
    args::Vector{NOTE_ARG_TYPE}
    retprobe::Bool
end
primitive type USDTContext{P,T1,T2,T3,T4,T5} Sys.WORD_SIZE end
@inline _context_reg(::Val{T}) where T = T
function infer_reg_type(size, signed)
    if size == 1
        signed ? Int8 : UInt8
    elseif size == 2
        signed ? Int16 : UInt16
    elseif size == 4
        signed ? Int32 : UInt32
    elseif size == 8
        signed ? Int64 : UInt64
    else
        NTuple{size,UInt8}
    end
end
@inline function API.get_param(ctx::USDTContext{P,T1,T2,T3,T4,T5}, ::Val{idx}) where {P,T1,T2,T3,T4,T5,idx}
    @assert 1<=idx<=5
    reg = idx == 1 ? T1 :
          idx == 2 ? T2 :
          idx == 3 ? T3 :
          idx == 4 ? T4 : T5
    reg = _context_reg(reg)
    ptr = reinterpret(P, ctx)
    reg, offset, size, signed = reg
    T = infer_reg_type(size, signed)
    if offset !== nothing
        ptr = unsafe_load(getproperty(ptr, reg))
        ptr = reinterpret(Ptr{T}, ptr) + offset
        dest = Ref{T}()
        GC.@preserve dest begin
            dest_ptr = Base.unsafe_convert(Ptr{T}, dest)
            if RT.probe_read_user(dest_ptr, sizeof(T), ptr) == 0
                return dest[]
            else
                return nothing
            end
        end
    else
        if reg === nothing
            # TODO: Broken argument, return garbage
            return zero(T)
        end
        return unsafe_load(reinterpret(Ptr{T}, getproperty(ptr, reg)))
    end
end
function USDT(f::Function, pid, binpath, note::NOTE_TYPE; merge_with=(), retprobe=false, kwargs...)
    func = note.func
    addr = note.location
    args = note.args
    Ts = [Val(length(args) >= i ? (args[i].reg, args[i].offset, args[i].size, args[i].signed) : nothing) for i in 1:5]
    P = API.pointertype(API.cpu_user_regs)
    C = USDTContext{P,Ts...}
    obj = API.Object(bpffunction(f, Tuple{C}; kwargs...))
    for other in merge_with
        merge_maps!(obj, other)
    end
    #foreach(prog->API.set_uprobe!(prog), API.programs(obj))
    foreach(prog->API.set_kprobe!(prog), API.programs(obj))
    USDT(obj, pid, binpath, func, addr, args, retprobe)
end

"Reads and returns the STAPSDT notes in the binary file `bin`."
function read_notes(bin)
    open(bin) do io
        elf = readmeta(io)
        sec = only(filter(x->section_name(x)==".note.stapsdt",
                          collect(Sections(elf))))
        sec_size = section_size(sec)

        notes = NOTE_TYPE[]
        seek(sec, 0)
        pos = 0
        off = position(io)
        while pos < sec_size
            # Extract fields
            type = read(io, UInt32)
            size = read(io, UInt32)
            _ = read(io, UInt32)
            owner = readuntil(io, '\0')
            location = read(io, UInt64)
            base = read(io, UInt64)
            semaphore = read(io, UInt64)
            provider = readuntil(io, '\0')
            func = readuntil(io, '\0')
            _args = readuntil(io, '\0')

            # Parse arguments
            _args = split(_args, ' ')
            args = NOTE_ARG_TYPE[]
            for arg in _args
                m = match(r"([\-0-9]*)@([$\-_0-9a-z]*)\(?%?([a-z]*[0-9]*)?\)?", arg)
                if m !== nothing
                    sz, offset, reg = m.captures
                    sz = parse(Int, sz)
                    signed = sz < 0
                    sz = abs(sz)
                    if startswith(offset, '$')
                        # TODO: Integer literal, what does this mean?
                        offset = nothing
                        reg = nothing
                    else
                        reg = Symbol(reg)
                        if offset == ""
                            offset = nothing
                        else
                            offset = try
                                parse(Int, offset)
                            catch
                                0 # TODO: Is this right? This is probably PC-relative?
                            end
                        end
                    end
                    push!(args, (;size=sz, signed, reg, offset))
                else
                    @assert isempty(arg) "Failed to parse argument: $arg"
                end
            end

            push!(notes, (;owner, location, base, semaphore, provider, func, args))

            # Seek to next note
            pos += size + sizeof(owner)+1 + (4*3)
            pos = cld(pos, 4) * 4
            seek(io, off+pos)
        end

        notes
    end
end
function USDT(f::Function, pid, binpath, provider::String, func::String; retprobe=false, multi=true, kwargs...)
    # FIXME: Find libs without bpftrace
    probes = String(read(`bpftrace -p $pid -l`))
    probe_rgx = Regex("^usdt:/proc/$pid/root(.*):$provider:$func\$")
    for probe in split(probes, '\n')
        startswith(probe, "usdt:") || continue
        if occursin(Regex("$provider:$func\$"), probe)
            m = match(probe_rgx, probe)
            @assert m !== nothing "Unexpected bpftrace probe format"
            probe_file = m.captures[1]
            notes = filter(x->x.func==func, read_notes(probe_file))
            if multi
                probes = ProbeSet()
                for note in filter(x->x.func==func, notes)
                    probe = USDT(f, pid, probe_file, note; retprobe, kwargs...)
                    push!(probes.probes, probe)
                end
                return probes
            else
                @assert length(notes) <= 1 "Multiple probe points found for $func"
                note = only(notes)
                return USDT(f, pid, probe_file, note; retprobe, kwargs...)
            end
        end
    end
    throw(ArgumentError("Failed to find $func in $binpath for process $pid"))
end
struct Tracepoint <: AbstractProbe
    obj::API.Object
    category::String
    name::String
end
function Tracepoint(f::Function, ::Type{T}, category, name; merge_with=(), kwargs...) where {T<:Tuple}
    obj = API.Object(bpffunction(f, T; kwargs...))
    foreach(prog->API.set_tracepoint!(prog), API.programs(obj))
    for other in merge_with
        merge_maps!(obj, other)
    end
    Tracepoint(obj, category, name)
end
Tracepoint(f::Function, category, name; kwargs...) =
    Tracepoint(f, Tuple{API.pointertype(API.pt_regs)}, category, name; kwargs...)

function Base.show(io::IO, p::ProbeSet)
    n = length(p.probes)
    println(io, "ProbeSet ($n probes):")
    for idx in 1:n
        print(io, "  ")
        show(io, p.probes[idx])
        idx < n && println(io)
    end
end
Base.show(io::IO, p::KProbe) =
    print(io, "KProbe ($(p.kfunc))")
Base.show(io::IO, p::UProbe) =
    print(io, "UProbe ($(p.addr) @ $(repr(p.binpath)))")
Base.show(io::IO, p::USDT) =
    print(io, "USDT ($(p.func) @ $(p.binpath) (pid $(p.pid)))")
Base.show(io::IO, p::Tracepoint) =
    print(io, "Tracepoint ($(p.category)/$(p.name))")

API.load(p::ProbeSet) = foreach(API.load, p.probes)
function API.load(p::KProbe)
    API.load(p.obj)
    foreach(prog->API.attach_kprobe!(prog, p.retprobe, p.kfunc),
            API.programs(p.obj))
end
#=
function API.load(p::UProbe)
    binpath = "$(Sys.BINDIR)/julia" # FIXME: Detect julia vs julia-debug
    ms = Base.method_instances(p.func, p.sig)
    @assert length(ms) == 1
    faddr = UInt64(ms[1].cache.specptr)
    iaddr = UInt64(dlsym(dlopen(binpath), :_init))
    addr = faddr - (iaddr - 0x1000) # FIXME: This doesn't work
    API.load(p.obj)
    foreach(prog->API.attach_uprobe!(prog, p.retprobe, UInt32(0), binpath, addr), API.programs(p.obj))
    API.load(p.obj)
    foreach(prog->API.attach_uprobe!(prog, p.retprobe, p.pid, p.binpath, p.addr), API.programs(p.obj))
end
=#
function API.load(p::USDT)
    API.load(p.obj)
    foreach(prog->API.attach_uprobe!(prog, p.retprobe, p.pid, p.binpath, p.addr),
            API.programs(p.obj))
end
function API.load(p::Tracepoint)
    API.load(p.obj)
    foreach(prog->API.attach_tracepoint!(prog, p.category, p.name),
            API.programs(p.obj))
end

function API.unload(p::AbstractProbe)
    API.unload(p.obj)
end
API.unload(p::ProbeSet) = foreach(API.unload, p.probes)


text:
using Distributed, StatsBase, Distributions, PyPlot
using TransitionalMCMC

lb  = -15
ub  = 15

fT(x) = pdf(Uniform(lb,ub), x[1,:])
sample_fT(Nsamples) = rand(Uniform(lb,ub), Nsamples, 1)

function log_fD_T(x)
    return log.(pdf(MvNormal([0],1), x) + pdf(MvNormal([5],0.2), x))
end

Nsamples = 2000
samps, acc =tmcmc_par(log_fD_T, fT, sample_fT, Nsamples)

plt.hist(samps,50)
plt.scatter(samps[:,1], samps[:,2])


text:
using SpWord2Vec
using ArgParse

function push_arg!(args, parsed_args, arg)
    if parsed_args[arg] != nothing
        push!(args, (Symbol(arg), parsed_args[arg]))
    end
end

function parse_commandline()
    s = ArgParseSettings(autofix_names=true)
    @add_arg_table s begin
        "--model", "-m"
            help = "Choose the model to train. Possible values are A, B, C, D and E"
            default = "A"
        "--corpus", "-c"
            help = "corpus file name"
            arg_type = String
            required = true
        "--lambda"
            help = "sparseness weight"
            arg_type = Float32
            default = 1.0f-2
        "--varepsilon"
            help = "L2 weight for matrix A in model [D]"
            arg_type = Float32
            default = 1.0f-2
        "--lr", "-a"
            help = "learning rate"
            arg_type = Float32
            default = 1.0f-3
        "--min-lr", "-b"
            help = "Minimum learning rate"
            arg_type = Float32
            default = 1.0f-4
        "--negative"
            help = "number negative samples per positive sample"
            arg_type = Int
            default = 15
        "--pretrain", "-p"
            help = "The ratio of data used to pre-train the sparse representations"
            arg_type = Float32
            default = 0.0f0
            dest_name = "preratio"
        "--every", "-e"
            help = "Dictionary parameters update frequencey"
            arg_type = Int
            default = 1
        "--normalize-sp"
            help = "normalize sparse parameter gradients or not"
            action = :store_true
            default = false
        "--normalize-dict"
            help = "normalize dictionary parameter gradients or not"
            action = :store_true
            default = false
        "--save-name"
            help = "save name"
            dest_name = "save_basename"
        "--mini-batch"
            help = "mini-batch size"
            arg_type = Int
            default = 1024
        "--saveratio"
            help = "save ratio"
            arg_type = Float32
            default = 0.05f0
        "--part"
            help = "training data ratio"
            arg_type = Float32
            default = 1.0f0
        "--vecdim", "-d"
            help = "words dimensions"
            arg_type = Int
            default = 200
        "--num-atom"
            help = "number of atoms"
            arg_type = Int
            default = 1024
        "--num-vocab"
            help = "number of vocabulary to keep"
            arg_type = Int
            default = 20_000
    end
    return parse_args(ARGS, s)
end

parsed_args = parse_commandline()
@eval SpSkipGram = SpWord2Vec.$(Symbol(:SpSkipGram, parsed_args["model"]))

train = SpSkipGram.train
init = SpSkipGram.init
close_device = SpSkipGram.close_device
@eval SparseSkipGram = SpSkipGram.$(Symbol(:SparseSkipGram, parsed_args["model"]))

corpus = parsed_args["corpus"]
d = parsed_args["vecdim"]
nbase = parsed_args["num_atom"]
nvocab = parsed_args["num_vocab"]
mini_batch = parsed_args["mini_batch"]
negative = parsed_args["negative"]
lr = parsed_args["lr"]
lambda = parsed_args["lambda"]
varepsilon = parsed_args["varepsilon"]

init(corpus)
args = Any[]

if parsed_args["model"] == "D"
    model = SparseSkipGram{Array}(
        nbase,              # number of atom words
        nvocab,             # number of vocabulary
        mini_batch,         # mini-batch size
        negative,           # number of negative samples per positive sample
        lr,                 # learning rate
        lambda,             # sparseness weight
        varepsilon,
        SpSkipGram.init_model)  # dumped init_model.
    for arg in String["min_lr", "normalize_sp", "saveratio", "part",
                      "save_basename"]
        push_arg!(args, parsed_args, arg)
    end
else
    model = SparseSkipGram{Array}(
        d,                  # vector dimension
        nbase,              # number of atom words
        nvocab,             # number of vocabulary
        mini_batch,         # mini-batch size
        negative,           # number of negative samples per positive sample
        lr,                 # learning rate
        lambda,             # sparseness weight
        SpSkipGram.init_model)  # dumped init_model.

    for arg in String["min_lr", "normalize_sp", "normalize_dict", "every",
                      "preratio", "saveratio", "part", "save_basename"]
        push_arg!(args, parsed_args, arg)
    end
end

train(model; args...)

close_device()


text:
using DataStructures
using GeometricExamples

const out_dir = "../../docs/src/diagnostics/poincare_1st/guiding_center_4d/"

# create list of figures to put in list
figures = (
    "_poincare_1st_q",
    "_poincare_1st_p",
    "_poincare_1st_l"
)

# suffix of figure files
fig_suff = ".png"

# collect simulation parameters
function get_params(Δt, ntime, problem)
    params = OrderedDict()
    params["Δt"] = Δt
    params["nt"] = ntime
    params["nloop"] = nloop
    params["nsave"] = nsave
    params["nplot"] = nplot
    params["GeometricIntegrators.jl"] = package_version("GeometricIntegrators")
    params["ChargedParticleDynamics.jl"] = package_version("ChargedParticleDynamics")
    params["MagneticEquilibria.jl"] = package_version("MagneticEquilibria")
    params
end

function run_tableau_list(tableaus, pinv, prefix, out_dir, fig_dir, nplot)
    # create output directory
    if !isdir(out_dir * "/" * fig_dir)
        mkpath(out_dir * "/" * fig_dir)
    end

    for tab in tableaus
        integrator = tab[1]
        tableau    = tab[2]
        runid      = tab[3]
        filename   = out_dir * "/" * fig_dir * "/" * prefix*runid*"_nl"*string(pinv.nloop)*"_nt"*string(pinv.ntime)

        sim = Simulation(pinv.equ, integrator, tableau, pinv.Δt, runid, filename * "_solution.h5", pinv.ntime, pinv.nsave)
        sol = run!(sim)

        info("Computing loop integrals...")

        I, J, L = evaluate_poincare_invariant(pinv, sol)

        write_to_hdf5(pinv, sol, filename * "_poincare_1st.h5")

        info("Plotting results...")

        plot_integral_error(sol.t.t, I, filename * "_poincare_1st_q" * fig_suff;
                            plot_title=L"$\Delta \int_\gamma \theta_i (q) \, dq^i$")

        if isdefined(sol, :p)
            plot_integral_error(sol.t.t, J, filename * "_poincare_1st_p" * fig_suff;
                                plot_title=L"$\Delta \int_\gamma p_i \, dq^i$")
        end

        if isdefined(sol, :λ)
            plot_integral_error(sol.t.t, L, filename * "_poincare_1st_l" * fig_suff;
                                plot_title=L"$\Delta \int_\gamma ( \theta_i (q) - \Delta t \, \lambda^{k} \theta_{k,i} (q) ) ( dq^i - \Delta t \, d\lambda^i )$")
        end

        plot_poincare_loop(sol, nplot, filename * "_loop" * fig_suff; xmin=XMIN, xmax=XMAX, ymin=YMIN, ymax=YMAX)
        plot_poincare_trajectories(sol, nplot, filename * "_trajectories" * fig_suff; xmin=XMIN, xmax=XMAX, ymin=YMIN, ymax=YMAX)

    end
end

run_tableau_list(tableau_list, pinv, prefix, out_dir, fig_dir, nplot)

# generate HTML file


text:
"""
# Available Functions

- `initialize`
- `terminate`
- `get_data`
- `lpf_on`
- `lpf_off`
- `get_lpf_state`
- `set_impedance_1Mohm`
- `set_impedance_50ohm`
- `get_impedance`

# Example
```
> scope_h = initialize(AgilentDSOX4034A)
```

Grab data from channel 1
```
> data = get_data(scope_h, 1)
```

Grab data from channel 2 and 4
```
> data_array = get_data(scope_h, [2,4])
```

Grab data from all available channels
Plot waves from all available channels
```
> using Plots
> data_array = get_data(scope_h)
> plot(data_struct)
```

Low Pass Filter
Turn on Low Pass Filter 25
```
> lpf_on(scope)
```

Check if low pass filter is on
```
> get_lpf_state(scope) == "1"
```

Turn on Low Pass Filter 25MHz
```
> lpf_off(scope)
> get_lpf_state(scope) == "0"
```

Impedance

```
> set_impedance_1Mohm(scope_h)
> get_impedance(scope_h) == ONEM
> set_impedance_50ohm(scope_h)
> get_impedance(scope_h) == FIFT
```

Terminate TCP connection
```
> terminate(scope_h)
```
"""
struct AgilentDSOX4034A <: Oscilloscope end


text:
export FunctionalScore

"""
    struct FunctionalScore{I} <: Score{I}

A score defined by a function.
"""
struct FunctionalScore{I} <: Score{I}
    fn :: Function # Function I => Double
end

@impl begin
    struct FunctionalScoreGetScore end
    function get_score(sf::FunctionalScore{I}, i::I)::AbstractFloat where {I}
        return sf.fn(i)
    end
end


text:
using Crazyflie
using MeshCat, GeometryTypes, CoordinateTransformations

function plotcflie(vis, uri="radio://0/80/2M")
    play(uri) do cf
        lconfig = logger.LogConfig(name="Configuration", period_in_ms=100)
        lconfig.add_variable("stateEstimateZ.x", "int16_t")
        lconfig.add_variable("stateEstimateZ.y", "int16_t")
        lconfig.add_variable("stateEstimateZ.z", "int16_t")
        lconfig.add_variable("stateEstimate.qx", "float")
        lconfig.add_variable("stateEstimate.qy", "float")
        lconfig.add_variable("stateEstimate.qz", "float")
        lconfig.add_variable("stateEstimate.qw", "float")
        record(cf, lconfig) do logs
            for entry in logs
                data = entry[2]
                pos = Translation(data["stateEstimateZ.x"]/1.0e3,
                                  data["stateEstimateZ.y"]/1.0e3,
                                  data["stateEstimateZ.z"]/1.0e3)
                quat = LinearMap(Quat(data["stateEstimate.qw"],
                                      data["stateEstimate.qx"],
                                      data["stateEstimate.qy"],
                                      data["stateEstimate.qz"]))
                print("qw = ", data["stateEstimate.qz"], "\n")
                settransform!(vis, pos ∘ quat)
            end
        end
    end
end

function drawcflie()
    vis = Visualizer()
    open(vis)
    setobject!(vis, HyperRectangle(Vec(-0.1, -0.1, 0), Vec(0.2, 0.2, 0.01)))
    return vis
end


text:
import HierarchicalUtils: NodeType, LeafNode, InnerNode, noderepr, children

NodeType(::Type{<:Union{Missing, ArrayNode, ArrayModel, LazyNode}}) = LeafNode()
NodeType(::Type{<:AbstractNode}) = InnerNode()
NodeType(::Type{<:AbstractMillModel}) = InnerNode()
NodeType(::Type{<:LazyModel}) = InnerNode()

noderepr(::T) where T <: Union{AbstractNode, AbstractMillModel} = "$(Base.nameof(T))"
noderepr(::Missing) = "∅"
noderepr(n::ArrayNode) = "ArrayNode$(size(n.data))"
noderepr(n::ArrayModel) = "ArrayModel($(n.m))"
noderepr(n::BagNode) = if ismissing(n.data) || nobs(n.data) == 0
    "BagNode with $(length(n.bags)) empty bag(s)"
else
    "BagNode with $(length(n.bags)) bag(s)"
end
noderepr(n::BagModel) = "BagModel ↦ $(repr("text/plain", n.a)) ↦ $(repr("text/plain", n.bm))"
noderepr(n::WeightedBagNode) = "WeightedNode with $(length(n.bags)) bag(s) and weights Σw = $(sum(n.weights))"
noderepr(n::AbstractProductNode) = "ProductNode"
noderepr(n::ProductModel) = "ProductModel ↦ $(noderepr(n.m))"
noderepr(n::MissingNode) = "Missing"
noderepr(n::MissingModel) = "Missing"
noderepr(n::LazyNode{N,D}) where {N,D} = (data_len = length(n.data); "$N $data_len item$(data_len > 1 ? "s" : "")")
noderepr(n::LazyNode{N,D}) where {N,D<:Nothing} = "$(N) ∅"
noderepr(n::LazyModel{Name}) where {Name} = "Lazy$(Name)Model"

children(n::AbstractBagNode) = (n.data,)
children(n::BagModel) = (n.im,)
children(n::ProductNode) = n.data
children(n::ProductModel) = n.ms
children(n::MissingNode) = (n.data,)
children(n::MissingModel) = (n.data,)
children(n::LazyModel) = (n.m,)