text:
# open-loop-bank-model.jl
# SimLynx Open Loop Processing Example

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const N_TELLERS = 2
const MEAN_INTERARRIVAL_TIME = 4.0
const MIN_SERVICE_TIME = 2.0
const MAX_SERVICE_TIME = 10.0

"Process to generate n customers arriving into the system."
@process generator(n::Integer) begin
    dist = Exponential(MEAN_INTERARRIVAL_TIME)
    for i = 1:n
        @schedule now customer(i)
        # We don't want to wait after the last customer, which could extend the
        # simulation time past the last customer leaving.
        if i < n
            work(rand(dist))
        end
    end
end

"The ith customer into the system."
@process customer(i::Integer) begin
    dist = Uniform(MIN_SERVICE_TIME, MAX_SERVICE_TIME)
    @with_resource tellers begin
        work(rand(dist))
    end
end

"Run the simulation for n customers."
function run_simulation(n1::Integer, n2::Integer)
    @assert n1 > 0 && n2 > 0
    println("SimLynx Open Loop Processing Example")
    println("Number of runs = $n1, number of customers = $n2")
    println("Single Queue, $N_TELLERS Teller Bank Model")
    @simulation begin
        avg_wait = Variable{Float64}(data=:tally, history=true)
        for i = 1:n1
            @simulation begin
                global tellers = Resource(N_TELLERS, "tellers")
                @schedule at 0.0 generator(n2)
                start_simulation()
                avg_wait.value = tellers.wait.stats.mean
            end
        end
        print_stats(avg_wait, title="Average Wait Statistics")
        plot_history(avg_wait, title="Average Wait History")
    end
end

run_simulation(10_000, 100)


text:
module Aircrafts
    include("aircrafts/C310.jl")
    include("aircrafts/F16.jl")

    include("trimmer.jl")
end

ext: jl


text:
using DiffEqBayes, BenchmarkTools

using OrdinaryDiffEq, RecursiveArrayTools, Distributions, ParameterizedFunctions, Mamba, BenchmarkTools
using Plots

gr(fmt=:png)

fitz = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

prob_ode_fitzhughnagumo = ODEProblem(fitz,[1.0,1.0],(0.0,10.0),[0.7,0.8,1/12.5,0.5])
sol = solve(prob_ode_fitzhughnagumo, Tsit5())

t = collect(range(1,stop=10,length=10))
sig = 0.20
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))

scatter(t, data[1,:])
scatter!(t, data[2,:])
plot!(sol)

priors = [Truncated(Normal(1.0,0.5),0,1.5),Truncated(Normal(1.0,0.5),0,1.5),Truncated(Normal(0.0,0.5),-0.5,0.5),Truncated(Normal(0.5,0.5),0,1)]

@time bayesian_result_stan = stan_inference(prob_ode_fitzhughnagumo,t,data,priors;reltol=1e-5,abstol=1e-5,vars =(StanODEData(),InverseGamma(3,2)))

plot_chain(bayesian_result_stan)

@time bayesian_result_turing = turing_inference(prob_ode_fitzhughnagumo,Tsit5(),t,data,priors)

plot_chain(bayesian_result_turing)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Pangomm"
version = v"2.49.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://download.gnome.org/sources/pangomm/$(version.major).$(version.minor)/pangomm-$(version).tar.xz",
                  "a2272883152618fddea016a62f50eb23b9b056ab3c08f3b64422591e6a507bd5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pangomm*/
if [[ "${target}" == "${MACHTYPE}" ]]; then
    # Delete host libexpat.so to avoid confusion
    rm /usr/lib/libexpat*
fi
mkdir output && cd output
meson --cross-file=${MESON_TARGET_TOOLCHAIN} ..
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))
# We don't have armv6l Pango at the moment
filter!(p -> arch(p) != "armv6l", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libpangomm-$(version.major)", "libpangomm-$(version.major).48"], :libpangomm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Cairomm_jll", uuid="af74c99f-f0eb-54aa-aecc-a10e8fc65c17"); compat="~1.16.1")
    Dependency(PackageSpec(name="Glibmm_jll", uuid="5d85a9da-21f7-5855-afec-cdc5039c46e8"); compat="~2.68.1")
    Dependency(PackageSpec(name="Pango_jll", uuid="36c8627f-9965-5494-a995-c6b170f724f3"); compat="1.47.0")
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
The **pcl_surface** library deals with reconstructing the original surfaces
from 3D scans.

http://docs.pointclouds.org/trunk/group__surface.html

## Exports

$(EXPORTS)
"""
module PCLSurface

using DocStringExtensions
using LibPCL
using PCLCommon
using Cxx

const libpcl_surface = LibPCL.find_library_e("libpcl_surface")
try
    Libdl.dlopen(libpcl_surface, Libdl.RTLD_GLOBAL)
catch e
    warn("You might need to set DYLD_LIBRARY_PATH to load dependencies proeprty.")
    rethrow(e)
end

end # module


text:
function r2cf(n1::Integer, n2::Integer)
    ret = Int[]
    while n2 != 0
        n1, (t1, n2) = n2, divrem(n1, n2)
        push!(ret, t1)
    end
    ret
end
r2cf(r::Rational) = r2cf(numerator(r), denominator(r))

function r2cf(n1, n2, maxiter=20)
    ret = Int[]
    while n2 != 0 && maxiter > 0
        n1, (t1, n2) = n2, divrem(n1, n2)
        push!(ret, t1)
        maxiter -= 1
    end
    ret
end

mutable struct NG
    a1::Int
    a::Int
    b1::Int
    b::Int
end

function ingress(ng, n)
    ng.a, ng.a1= ng.a1, ng.a + ng.a1 * n
    ng.b, ng.b1 = ng.b1, ng.b + ng.b1 * n
end

needterm(ng) = ng.b == 0 || ng.b1 == 0 || !(ng.a // ng.b == ng.a1 // ng.b1)

function egress(ng)
    n = ng.a // ng.b
    ng.a, ng.b = ng.b, ng.a - ng.b * n
    ng.a1, ng.b1 = ng.b1, ng.a1 - ng.b1 * n
    r2cf(n)
end

egress_done(ng) = (if needterm(ng) ng.a, ng.b = ng.a1, ng.b1 end; egress(ng))

done(ng) = ng.b == 0 && ng.b1 == 0

function testng()
    data = [["[1;5,2] + 1/2",      [2,1,0,2], [13,11]],
        ["[3;7] + 1/2",        [2,1,0,2], [22, 7]],
        ["[3;7] divided by 4", [1,0,0,4], [22, 7]],
        ["[1;1] divided by sqrt(2)", [0,1,1,0], [1,sqrt(2)]]]

    for d in data
        str, ng, r = d[1], NG(d[2]...), d[3]
        print(rpad(str, 25), "->")
        for n in r2cf(r...)
            if !needterm(ng)
                print(" $(egress(ng))")
            end
            ingress(ng, n)
        end
        while true
            print(" $(egress_done(ng))")
            if done(ng)
                println()
                break
            end
        end
    end
end

testng()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PyPlot
using SpecialFunctions: erfcx

"""
    ψ1(α, ϵ, N=100)

Approximates `(1+ϵ)^α - 1` by its Taylor expansion about `ϵ=0` using up to
`N` terms.
"""
function ψ1(α::T, ϵ::Complex{T}, N=100) where T <: AbstractFloat
    s = Complex(α, zero(T))
    term = s
    for n = 1:N
        term *= -( (n-α)/(n+1) ) * ϵ
        s += term
        if abs(term) < eps(T)
            break
        end
    end
    if abs(term) > 10*eps(T)
        error("ψ1 expansion did not converge after $N terms")
    end
    return s
end

"""
    ψ2(α, ϵ, N=100)

Approximates `(1+ϵ)^α - (1+αϵ)` by its Taylor expansion about `ϵ=0` using up
to `N` terms.
"""
function ψ2(α::T, ϵ::Complex{T}, N=100) where T <: AbstractFloat
    s = Complex(-α * (1-α) / 2, zero(T))
    term = s
    for n = 1:N
        term *= -( (n+1-α)/(n+2) ) * ϵ
        s += term
        if abs(term) < eps(T)
            break
        end
    end
    if abs(term) > 10*eps(T)
        error("ψ1 expansion did not converge after $N terms")
    end
    return s
end

function H(α::T, β::T, w::Complex{T}, x::T, sep::T) where T <: AbstractFloat
    xa = x^(1/α)
    ϵ = ( w - xa ) / xa
    if abs(ϵ) > sep
        return w^(α-β) / (w^α-x) - 1 / (α*ϵ*x^(β/α))
    else
        return ( ψ1(α-β,ϵ) - (1+ϵ)^(α-β) * ψ2(α,ϵ)/ψ1(α,ϵ) ) / (α*x^(β/α))
    end
end

a(ϕ) = acosh(2ϕ/((4ϕ-π)*sin(ϕ)))
b(ϕ) = ( (4π*ϕ-π^2) / a(ϕ) )

function Qpos(α::T, β::T, x::T, N::Integer, sep::T) where T <: AbstractFloat
    ϕ = parse(T, "1.17210")
    h = a(ϕ) / N
    μ = b(ϕ) * N
    w0r = μ * ( 1 - sin(ϕ) )
    w0 = Complex(w0r, zero(T))
    s = ( exp(w0r) * cos(ϕ) / 2 ) * real(H(α, β, w0, x, sep))
    for n = 1:N
        un = n * h
        iunmϕ = Complex(-ϕ, un)
        wn = μ * ( 1 + sin(iunmϕ) )
        s += real( exp(wn) * H(α, β, wn, x, sep) * cos(iunmϕ) )
    end
    return ( x^((1-β)/α) * exp(x^(1/α)) / α ) + ( a(ϕ) * b(ϕ) / π ) * s
end

function Hplus(α::T, β::T, w::Complex{T}, x::T, sep::T) where T <: AbstractFloat
    γp = x^(1/α) * exp(complex(zero(T), π/α))
    ϵp = ( w - γp ) / γp
    if abs(ϵp) > sep
        return w^(α-β) / ( w^α + x ) - γp^(1-β) / ( α * ϵp * γp^β )
    else
        return ( ψ1(α-β,ϵp) - (1+ϵp)^(α-β)*ψ2(α,ϵp) / ψ1(α,ϵp) ) / (α*γp^β)
    end
end

function Hminus(α::T, β::T, w::Complex{T}, x::T, sep::T
               ) where T <: AbstractFloat
    γm= x^(1/α) * exp(complex(zero(T), -π/α))
    ϵm = ( w - γm) / γm
    if abs(ϵm) > sep
        return w^(α-β) / ( w^α + x ) - γm^(1-β) / ( α * ϵm * γm^β )
    else
        return ( ψ1(α-β,ϵm) - (1+ϵm)^(α-β)*ψ2(α,ϵm) / ψ1(α,ϵm) ) / (α*γm^β)
    end
end

function Qneg(α::T, β::T, x::T, N::Integer, sep::T) where T <: AbstractFloat
    @assert 1 < α < 2
    ϕ = parse(T, "1.17210")
    h = a(ϕ) / N
    μ = b(ϕ) * N
    w0r = μ * ( 1 - sin(ϕ) )
    w0 = Complex(w0r, zero(T))
    s = ( exp(w0r) * cos(ϕ) / 2 ) * real(Hplus(α, β, w0, x, sep))
    for n = 1:N
        un = n * h
        iunmϕ = Complex(-ϕ, un)
        wn = μ * ( 1 + sin(iunmϕ) )
        s += real( exp(wn) * ( Hplus(α, β, wn, x, sep)
                             + Hminus(α, β, wn, x, sep) ) * cos(iunmϕ) ) / 2
    end
    sum_residues = ( ( x^((1-β)/α) / α ) * exp(x^(1/α)*cos(π/α))
                     *cos( π*(1-β)/α + x^(1/α)*sin(π/α) ) )
    return sum_residues + ( a(ϕ)*b(ϕ)/π ) * s
end

α = 0.5
β = 1.0
sep = 0.2
N = 10

figure(1)
x = 2.0
t = range(-2.0, 2.0, length=201)
θ = π/4
w = x^(1/α) .+ t * exp(Complex(0.0, θ))
Hvals = H.(α, β, w, x, sep)
plot(t, real.(Hvals), t, imag.(Hvals))
grid(true)

figure(2)
x = range(0.0, 1.0, length=201)
plot(x, Qpos.(α, β, x, N, sep) - erfcx.(-x))
grid(true)

figure(3)
α = 3/2
x = 2.0
γplus = x^(1/α) * exp(Complex(0.0,π/α))
t = range(-2.0, 2.0, length=201)
θ = π/4
w = x^(1/α) .+ t * exp(Complex(0.0, θ))
Hvals = Hplus.(α, β, w, x, sep)
plot(t, real.(Hvals), t, imag.(Hvals))
grid(true)

figure(4)
x = range(0.0, 5.0, length=201)
plot(x, Qneg.(α, β, x, N, sep))
grid(true)


text:
"""
Additional activity selector

returns a category of additional activity choosen by agent.

**Arguments**
* `feature_classes` : dictionary containing all possible types of features in the simulation

** Assumptions
 - there is maximum one waypoint
 -probability of driving directly to work is equal to 0.5
 -otherwise waypoint type is chosen randomly from the features categories
"""
function additional_activity(agent_profile::DataFrames.DataFrame, before::Bool,
	                         sim_data::OpenStreetMapXSim.SimData)
	if rand() < 0.5
		return nothing
	else
		return rand(keys(sim_data.feature_classes))
	end
end


text:
function POMDPs.transition(mdp::GenerativeMergingMDP, s::AugScene, a::Int64)
    # Pr(sp |s, a)
    rng = MersenneTwister(1)
    spred = gen(DDNOut(:sp), mdp, s, a, rng)
    s_vec = extract_features(mdp, spred)
    sig_p = diagm(0 => [1.0 for i=1:length(s_vec)])
    return  d = MultivariateNormal(s_vec, sig_p)
end

function POMDPs.observation(mdp::GenerativeMergingMDP, a::Int64, s)
    sig_o = diagm(0 => [0.5 for i=1:length(s)])
    return MultivariateNormal(s, sig_o)
end

"""
    MergingBelief
A type to represent belief state. It consists of the current observation `o` and the estimated driver types represented
by a dictionnary mapping ID to cooperation levels.
# fields
- `o::AugScene`
- `driver_types::Dict{Int64, Float64}`
"""
struct MergingBelief
    o::AugScene
    driver_types::Dict{Int64, Float64}
end

function POMDPs.convert_s(t::Type{V}, b::MergingBelief, mdp::GenerativeMergingMDP) where V<:AbstractArray
    ovec = convert_s(t, b.o, mdp)
    fore, merge, fore_main, rear_main = get_neighbors(mdp.env, b.o.scene, EGO_ID)
    if fore.ind != nothing
        fore_id = b.o.scene[fore.ind].id
        ovec[6] = b.driver_types[fore_id] > 0.5
    end
    if merge.ind != nothing
        merge_id = b.o.scene[merge.ind].id
        ovec[9] = b.driver_types[merge_id] > 0.5
    end
    if fore_main.ind != nothing
        fore_main_id = b.o.scene[fore_main.ind].id
        ovec[12] = b.driver_types[fore_main_id] > 0.5
    end
    if rear_main.ind != nothing
        rear_main_id = b.o.scene[rear_main.ind].id
        ovec[15] = b.driver_types[rear_main_id] > 0.5
    end
    return ovec
end

function belief_weight(t::Type{V}, b::MergingBelief, mdp::GenerativeMergingMDP, state) where V <: AbstractArray
    ovec = convert_s(t, b.o, mdp)
    fore, merge, fore_main, rear_main = get_neighbors(mdp.env, b.o.scene, EGO_ID)
    weight = 1.0
    if fore.ind != nothing
        fore_id = b.o.scene[fore.ind].id
        ovec[6] = state[1] #b.driver_types[fore_id] > 0.5
        weight *= b.driver_types[fore_id]*state[1] + (1 - b.driver_types[fore_id])*(1 - state[1])
    end
    if merge.ind != nothing
        merge_id = b.o.scene[merge.ind].id
        ovec[9] = state[2] #b.driver_types[merge_id] > 0.5
        weight *= b.driver_types[merge_id]*state[2] + (1 - b.driver_types[merge_id])*(1 - state[2])
    end
    if fore_main.ind != nothing
        fore_main_id = b.o.scene[fore_main.ind].id
        ovec[12] = state[3] #b.driver_types[fore_main_id] > 0.5
        weight *= b.driver_types[fore_main_id]*state[3] + (1 - b.driver_types[fore_main_id])*(1 - state[3])
    end
    if rear_main.ind != nothing
        rear_main_id = b.o.scene[rear_main.ind].id
        ovec[15] = state[4] #b.driver_types[rear_main_id] > 0.5
        weight *= b.driver_types[rear_main_id]*state[4] + (1 - b.driver_types[rear_main_id])*(1 - state[4])
    end
    return ovec, weight
end

const driver_type_states =  collect(Iterators.product([[0,1] for i=1:4]...))

"""
    MergingUpdater <: Updater
A belief updater for `MergingBelief`. It sets `o` to the current observation and updates
the belief on the drivers cooperation level using Bayes' rule
"""
struct MergingUpdater <: Updater
    # must set observe_cooperation to false
    # must fill in all the driver models and make a deepcopy
    mdp::GenerativeMergingMDP

    function MergingUpdater(mdp::GenerativeMergingMDP)
        mdp_ = deepcopy(mdp)
        mdp_.observe_cooperation = false
        for i=1:mdp.max_cars
            mdp_.driver_models[i+1] = CooperativeIDM()
        end
        return new(mdp_)
    end
end

function POMDPs.update(up::MergingUpdater, b_old::MergingBelief, a::Int64, o::AugScene)
    # b_neigh = update_neighbors(up.mdp, b_old, o)
    driver_types = deepcopy(b_old.driver_types)
    for i=2:up.mdp.max_cars+1
        update_proba!(up.mdp, b_old, driver_types, a, o, i)
    end
    return MergingBelief(o, driver_types)
end

function update_neighbors(mdp::GenerativeMergingMDP, b::MergingBelief, o::AugScene)
    fore, merge, fore_main, rear_main = get_neighbors(mdp.env, o.scene, EGO_ID)
    bfore = b.fore
    bmerge = b.merge
    bforemain = b.fore_main
    brearmain = b.rear_main
    current_neighbors = [bfore.id, bmerge.id, bforemain.id, brearmain.id]
    if fore.ind == nothing
        bfore = (id=nothing, prob=0.5)
    elseif o.scene[fore.ind].id != b.fore.id
        bfore = (id=o.scene[fore.ind].id, prob=0.5)
    end
    if merge.ind == nothing
        bmerge = (id=nothing, prob=0.5)
    elseif o.scene[merge.ind].id != b.merge.id
        bmerge = (id=o.scene[merge.ind].id, prob=0.5)
    end
    if fore_main.ind == nothing
        bforemain = (id=nothing, prob=0.5)
    elseif o.scene[fore_main.ind].id != b.fore_main.id
        bforemain = (id=o.scene[fore_main.ind].id, prob=0.5)
    end
    if rear_main.ind == nothing
        bforerear = (id=nothing, prob=0.5)
    elseif o.scene[rear_main.ind].id != b.rear_main.id
        brearmain = (id=o.scene[rear_main.ind].id, prob=0.5)
    end
    return MergingBelief(b.o, bfore, bmerge, bforemain, brearmain)
end

function update_proba!(mdp::GenerativeMergingMDP, b::MergingBelief, driver_types::Dict, a::Int64, o::AugScene, id::Int64)
    sp_vec = extract_features(mdp, o)
    probs = zeros(2)
    sp_vec = extract_features(mdp, o)
    old_prob = driver_types[id]
    if maximum(old_prob) ≈ 1.0
        maximum(old_prob)
        return driver_types
    end
    probs = zeros(2)
    for c in [0, 1]
        mdp.driver_models[id].c = c
        v_des = 5.0
        set_desired_speed!(mdp.driver_models[id], v_des)
        d = transition(mdp, b.o, a)
        probs[Int(c + 1)] += pdf(d, sp_vec)*(c*old_prob + (1 - c)*(1 - old_prob))
    end
    # @printf("probs = %s id=%d \n", probs, id)
    if sum(probs) ≈ 0.0
        probs = [0.5, 0.5]
    end
    normalize!(probs, 1)
    driver_types[id] = probs[2]
    return driver_types
end

function BeliefUpdaters.initialize_belief(up::MergingUpdater, s0::AugScene)
    driver_types = Dict{Int64, Float64}()
    for i=2:up.mdp.max_cars+1
        driver_types[i] = 0.5
    end
    b0 = MergingBelief(s0, driver_types)
    return b0
end


text:
# ---
# title: 480. Sliding Window Median
# id: problem480
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Sliding Window
# link: <https://leetcode.com/problems/sliding-window-median/description/>
# hidden: true
# ---
#
# Median is the middle value in an ordered integer list. If the size of the list
# is even, there is no middle value. So the median is the mean of the two middle
# value.
#
# Examples:
#
# `[2,3,4]` , the median is `3`
#
# `[2,3]`, the median is `(2 + 3) / 2 = 2.5`
#
# Given an array _nums_ , there is a sliding window of size _k_ which is moving
# from the very left of the array to the very right. You can only see the _k_
# numbers in the window. Each time the sliding window moves right by one
# position. Your job is to output the median array for each window in the
# original array.
#
# For example,
# Given _nums_ = `[1,3,-1,-3,5,3,6,7]`, and _k_ = 3.
#
#
#
#     Window position                Median
#     ---------------               -----
#     [1  3  -1] -3  5  3  6  7       1
#      1 [3  -1  -3] 5  3  6  7       -1
#      1  3 [-1  -3  5] 3  6  7       -1
#      1  3  -1 [-3  5  3] 6  7       3
#      1  3  -1  -3 [5  3  6] 7       5
#      1  3  -1  -3  5 [3  6  7]      6
#
#
# Therefore, return the median sliding window as `[1,-1,-1,3,5,6]`.
#
# **Note:**
# You may assume `k` is always valid, ie: `k` is always smaller than input
# array's size for non-empty array.
# Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

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


text:
# module for testing each common modules

module TestCommon
  using Test

  # target modules
  include(joinpath(split(@__FILE__, "test")[1], "src/common/covariance_ellipse/covariance_ellipse.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/common/error_calculation/error_calculation.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/common/state_transition/state_transition.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/common/observation_function/observation_function.jl"))

  function main()
    @testset "Common" begin
      @testset "CovarianceEllipse" begin
        @test_nowarn draw_covariance_ellipse!([0.0, 0.0, 0.0],
                                              [1.0 0.0 0.0;
                                               0.0 1.0 0.0;
                                               0.0 0.0 1.0], 3)
      end
      @testset "ErrorCalculation" begin
        error = calc_lon_lat_error([5.0, 5.0, 5.0], [0.0, 0.0, 0.0])
        @test error[1] == 5.00
        @test error[2] == 5.00
        @test error[3] == 5.00
      end
      @testset "StateTransition" begin
        @test state_transition(0.1, 0.0, 1.0, [0, 0, 0]) == [0.1, 0.0, 0.0]
        @test state_transition(0.1, 10.0/180*pi, 9.0, [0, 0, 0]) == [0.5729577951308232, 0.5729577951308231, 1.5707963267948966]
        @test state_transition(0.1, 10.0/180*pi, 18.0, [0, 0, 0]) == [7.016709298534876e-17, 1.1459155902616465, 3.141592653589793]
      end
      @testset "ObservationFunction" begin
        @test observation_function([-2, -1, pi/5*6], [2.0, -2.0]) == [4.123105625617661, 2.2682954597449703]
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using BayesianDataFusion
using Distributed
using SparseArrays
using Test

X = sprand(50, 100, 0.1)
Y1 = rand(50, 3)
Y2 = rand(100, 2)

Z1 = X' * Y1
Z2 = X  * Y2
addprocs(2)
@everywhere using BayesianDataFusion

fp  = psparse(X, workers())
Z1p = transpose(fp) * Y1
Z2p = fp * Y2

@test Z1 ≈ Z1p
@test Z2 ≈ Z2p

W  = sprand(50, 10, 0.1)
rd = RelationData(W, class_cut = 0.5, feat1 = fp)
assignToTest!(rd.relations[1], 2)

result = macau(rd, burnin = 2, psamples = 2, num_latent=5, verbose=false)

@test size(rd.entities[1].model.beta) == (size(fp,2), 5)

########### parallel sparse with CSR ###########
fp2 = psparse(SparseMatrixCSC(X'), workers())

Z1r = transpose(fp2) * Y1
Z2r = fp2 * Y2

@test Z1 ≈ Z1r
@test Z2 ≈ Z2r

rd2 = RelationData(W, class_cut = 0.5, feat1 = fp2)
assignToTest!(rd2.relations[1], 2)

result = macau(rd2, burnin = 2, psamples = 2, num_latent=5, verbose=false)

rmprocs( workers() )


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    democurve()

Shows two simple test curves.

# Examples:
```julia
julia> democurve()
```
"""
function democurve()
crv = nrbtestcrv();

# plot the control points
plot(crv.coefs[1,:],crv.coefs[2,:],
     title = "Arbitrary Test 2D Curves.",
     m = (7, :transparent, stroke(1, :red)),
     lw =1.0,
     c = :red,
     linestyle =:dash,
     framestyle=:box,
     legend = false);

# plot the nurbs curve
nrbplot!(crv,48);

# modify the curve
crv.knots[4] = 0.1;

# plot the nurbs curve
nrbplot!(crv,48);
end # democurve


text:
using Test
using CLIMAParameters
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Planet

@testset "Planet (Earth)" begin
  struct EarthParameterSet <: AbstractEarthParameterSet end

  ps = EarthParameterSet()

  # Planet
  @test cp_d(ps)           ≈ R_d(ps) / kappa_d(ps)
  @test cv_d(ps)           ≈ cp_d(ps) - R_d(ps)
  @test molmass_ratio(ps)  ≈ molmass_dryair(ps) / molmass_water(ps)
  @test cv_v(ps)           ≈ cp_v(ps) - R_v(ps)
  @test cv_l(ps)           ≈ cp_l(ps)
  @test cv_i(ps)           ≈ cp_i(ps)
  @test T_0(ps)            ≈ T_triple(ps)
  @test LH_f0(ps)          ≈ LH_s0(ps) - LH_v0(ps)
  @test e_int_v0(ps)       ≈ LH_v0(ps) - R_v(ps) * T_0(ps)
  @test e_int_i0(ps)       ≈ LH_f0(ps)
  @test R_d(ps)            ≈ gas_constant() / molmass_dryair(ps)
  @test R_v(ps)            ≈ gas_constant() / molmass_water(ps)
  @test year_anom(ps)      ≈ 365.26 * day(ps)
  @test orbit_semimaj(ps)  ≈ 1 * astro_unit()

  @test !isnan(T_freeze(ps))
  @test !isnan(T_min(ps))
  @test !isnan(T_max(ps))
  @test !isnan(T_icenuc(ps))
  @test !isnan(pow_icenuc(ps))
  @test !isnan(press_triple(ps))
  @test !isnan(surface_tension_coeff(ps))

end


text:
const AGGR2FUN = Dict(:add => sum, :sub => sum, :mul => prod, :div => prod,
                      :max => maximum, :min => minimum, :mean => mean)

abstract type GraphNet end

@inline update_edge(gn::T, e, vi, vj, u) where {T<:GraphNet} = e
@inline update_vertex(gn::T, ē, vi, u) where {T<:GraphNet} = vi
@inline update_global(gn::T, ē, v̄, u) where {T<:GraphNet} = u

@inline function update_batch_edge(gn::T, adj, E, V, u) where {T<:GraphNet}
    edge_idx = edge_index_table(adj)
    E_ = Vector[]
    for (i, js) = enumerate(adj)
        for j = js
            k = edge_idx[(i,j)]
            e = update_edge(gn, get_feature(E, k), get_feature(V, i), get_feature(V, j), u)
            push!(E_, e)
        end
    end
    hcat(E_...)
end

@inline function update_batch_vertex(gn::T, Ē, V, u) where {T<:GraphNet}
    V_ = Vector[]
    for i = 1:size(V,2)
        v = update_vertex(gn, get_feature(Ē, i), get_feature(V, i), u)
        push!(V_, v)
    end
    hcat(V_...)
end

@inline function aggregate_neighbors(gn::T, aggr::Symbol, E, accu_edge, num_V, num_E) where {T<:GraphNet}
    @assert !iszero(accu_edge) "accumulated edge must not be zero."
    cluster = generate_cluster(E, accu_edge, num_V, num_E)
    pool(aggr, cluster, E)
end

@inline function aggregate_neighbors(gn::T, aggr::Nothing, E, accu_edge, num_V, num_E) where {T<:GraphNet}
    @nospecialize E accu_edge num_V num_E
end

@inline function aggregate_edges(gn::T, aggr::Symbol, E) where {T<:GraphNet}
    u = vec(AGGR2FUN[aggr](E, dims=2))
    aggr == :sub && (u = -u)
    aggr == :div && (u = 1 ./ u)
    u
end

@inline function aggregate_edges(gn::T, aggr::Nothing, E) where {T<:GraphNet}
    @nospecialize E
end

@inline function aggregate_vertices(gn::T, aggr::Symbol, V) where {T<:GraphNet}
    u = vec(AGGR2FUN[aggr](V, dims=2))
    aggr == :sub && (u = -u)
    aggr == :div && (u = 1 ./ u)
    u
end

@inline function aggregate_vertices(gn::T, aggr::Nothing, V) where {T<:GraphNet}
    @nospecialize V
end

function propagate(gn::T, fg::FeaturedGraph, naggr=nothing, eaggr=nothing, vaggr=nothing) where {T<:GraphNet}
    adj = adjacency_list(fg)
    num_V = nv(fg)
    accu_edge = accumulated_edges(adj)
    num_E = accu_edge[end]
    E = edge_feature(fg)
    V = node_feature(fg)
    u = global_feature(fg)

    E = update_batch_edge(gn, adj, E, V, u)

    Ē = aggregate_neighbors(gn, naggr, E, accu_edge, num_V, num_E)

    V = update_batch_vertex(gn, Ē, V, u)

    ē = aggregate_edges(gn, eaggr, E)

    v̄ = aggregate_vertices(gn, vaggr, V)

    u = update_global(gn, ē, v̄, u)

    FeaturedGraph(graph(fg), V, E, u)
end


text:
# Squared Exponential Covariance Function

include("se_iso.jl")
include("se_ard.jl")

"""
    SE(ll::Union{Float64,Vector{Float64}}, lσ::Float64)

Create squared exponential kernel with length scale `exp.(ll)` and signal standard deviation
`exp(lσ)`.

See also [`SEIso`](@ref) and [`SEArd`](@ref).
"""
SE(ll::Float64, lσ::Float64) = SEIso(ll, lσ)
SE(ll::Vector{Float64}, lσ::Float64) = SEArd(ll, lσ)


text:
# module Bukdu.Plug

struct Parsers <: AbstractPlug
end

"""
    plug(::Type{Parsers}, parsers::Vector{Symbol} = ContentParsers.default_content_parsers; decoders...)
"""
function plug(::Type{Parsers}, parsers::Vector{Symbol} = ContentParsers.default_content_parsers; decoders...)
    ContentParsers.env[:decoders] = merge(ContentParsers.default_content_decoders, Dict{Symbol,Type{<:ContentParsers.AbstractDecoder}}(decoders))
    ContentParsers.env[:parsers] = union(parsers, first.(collect(decoders)))
end

# module Bukdu.Plug


text:
"""
    box_counting_dim(xg, yg, basin; kwargs...)
This function compute the box-counting dimension of the boundary. It is related to the uncertainty dimension.

[C. Grebogi, S. W. McDonald, E. Ott, J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983]

## Arguments
* `basin` : the matrix containing the information of the basin.
* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test.

## Keyword arguments
* `kwargs...` these arguments are passed to `generalized_dim` see `ChaosTools.jl` for the docs.

"""
function box_counting_dim(xg, yg, basin; kwargs...)
    # get points coordinates in the boudary
    bnd = get_boundary_filt(basin)
    I1 = findall(bnd .== 1);
    v = Dataset{2,eltype(xg)}()
    for  ind in I1; push!(v, [xg[ind[1]] , yg[ind[2]]]); end;
    return generalized_dim(v; q=0, kwargs...)
end

function box_counting_dim(xg, yg, bsn::BasinInfo; kwargs...)
    ind  = findall(iseven.(bsn.basin) .== true)
    # Not sure it is necessary to deepcopy
    basin_test = deepcopy(bsn.basin)
    # Remove attractors from the picture
    for k in ind; basin_test[k] = basin_test[k]+1; end
    return box_counting_dim(xg, yg, basin_test; kwargs...)
end

"""
    uncertainty_exponent(bsn::BasinInfo, integ; precision=1e-3) -> α,e,ε,f_ε
This function estimates the uncertainty exponent of the boundary. It is related to the uncertainty dimension.

[C. Grebogi, S. W. McDonald, E. Ott, J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983]

## Arguments
* `bsn` : structure that holds the information of the basin.
* `integ` : handle to the iterator of the dynamical system.

## Keyword arguments
* `precision` variance of the estimator of the uncertainty function.

"""
function uncertainty_exponent(bsn::BasinInfo, integ; precision=1e-3)
    xg = bsn.xg; yg = bsn.yg;
    nx=length(xg)
    ny=length(yg)
    xi = xg[1]; xf = xg[end];
    grid_res_x=xg[2]-xg[1]

    num_step=10
    N_u = zeros(1,num_step) # number of uncertain box
    N = zeros(1,num_step) # number of boxes
    ε = zeros(1,num_step) # resolution

    # First we compute a crude estimate to have an idea at which resolution
    D,e,f = static_estimate(xg,yg,bsn);
    println("First estimate: α=",D)

    # resolution in log scale
    # estimate the minimum resolution to estimate a f_ε ∼ 10^-3
    min_ε = -3/D[1]
    @show min_ε = max(min_ε,-7) # limit resolution to 10^-7 to avoid numerical inconsistencies.
    max_ε = log10(grid_res_x*5) # max radius of the ball is roughly 5 pixels.

    r_ε = 10 .^ range(min_ε,max_ε,length=num_step)

    for (k,eps) in enumerate(r_ε)
        Nb=0; Nu=0; μ=0; σ²=0; M₂=0;
        completed = 0;
        # Find uncertain boxes
        while completed == 0
            k1 = rand(1:nx)
            k2 = rand(1:ny)

            c1 = bsn.basin[k1,k2]
            c2 = get_color_point!(bsn, integ, [xg[k1]+eps,yg[k2]])
            c3 = get_color_point!(bsn, integ, [xg[k1]-eps,yg[k2]])

            if length(unique([c1,Int(c2),Int(c3)]))>1
                Nu = Nu + 1
            end
            Nb += 1

            # Welford's online average estimation and variance of the estimator
            M₂ = wel_var(M₂, μ, Nu/Nb, Nb)
            μ = wel_mean(μ, Nu/Nb, Nb)
            σ² = M₂/Nb
            # If the process is binomial, the variance of the estimator is Nu/Nb·(1-Nu/Nb)/Nb (p·q/N)
            # simulations matches
            #push!(sig,Nu/Nb*(1-Nu/Nb)/Nb)

            # Stopping criterion: variance of the estimator of the mean bellow 1e-3
            if Nu > 50 && σ² < precision
                completed = 1
                @show Nu,Nb,σ²
            end

            if Nu < 10 && Nb>10000
                # skip this value, we don't find the boundary
                # corresponds to roughly f_ε ∼ 10^-4
                @warn "Resolution may be too small for this basin, skip value"
                completed = 1
                Nu = 0
            end

        end
        N_u[k]=Nu
        N[k]=Nb
        ε[k]=eps
    end

    # uncertain function
    f_ε = N_u./N

    # remove zeros in case there are:
    ind = f_ε .> 0
    f_ε =  f_ε[ind]
    ε = ε[ind]

    # @show i1,D=linear_region(vec(log10.(ε)), vec(log10.(f_ε)))
    # Dϵ=0.
    # @show i1
    # get exponent from a linear region
    @. model(x, p) = p[1]*x+p[2]
    fit = curve_fit(model, vec(log10.(ε)), vec(log10.(f_ε)), [2., 2.])
    D = coef(fit)
    Dϵ = estimate_errors(fit)
    return D[1], Dϵ[1], vec(log10.(ε)),vec(log10.(f_ε))
end

function wel_var(M₂, μ, xₙ, n)

    μ₂ = μ + (xₙ - μ)/n
    M₂ = M₂ + (xₙ - μ)*(xₙ - μ₂)
    return M₂
end

function wel_mean(μ, xₙ, n)
    return μ + (xₙ - μ)/n
end

function static_estimate(xg,yg,bsn; precision=1e-4)

    xg = bsn.xg; yg = bsn.yg;
    y_grid_res=yg[2]-yg[1]
    nx=length(xg)
    ny=length(yg)

    num_step=10
    N_u = zeros(1,num_step) # number of uncertain box
    N = zeros(1,num_step) # number of boxes
    ε = zeros(1,num_step) # resolution

    # resolution in pixels
    min_ε = 1;
    max_ε = 10;

    r_ε = min_ε:max_ε
    #@show r_ε = 10 .^ range(log10(min_ε),log10(max_ε),length=num_step)
    for (k,eps) in enumerate(r_ε)
        Nb=0; Nu=0; μ=0; σ²=0; M₂=0;
        completed = 0;
        # Find uncertain boxes
        while completed == 0
            kx = rand(1:nx)
            ky = rand(ceil(Int64,eps+1):floor(Int64,ny-eps))

            indy = range(ky-eps,ky+eps,step=1)
            c = [bsn.basin[kx,ky] for ky in indy]

            if length(unique(c))>1
                Nu = Nu + 1
            end
            Nb += 1

            # Welford's online average estimation and variance of the estimator
            M₂ = wel_var(M₂, μ, Nu/Nb, Nb)
            μ = wel_mean(μ, Nu/Nb, Nb)
            σ² = M₂/Nb

            # Stopping criterion: variance of the estimator of the mean bellow 1e-3
            if Nu > 50 && σ² < precision
                completed = 1
                #@show Nu,Nb,σ²
            end

        end
        N_u[k]=Nu
        N[k]=Nb
        ε[k]=eps*y_grid_res
    end

    # uncertain function
    f_ε = N_u./N

    # remove zeros in case there are:
    ind = f_ε .> 0.
    f_ε =  f_ε[ind]
    ε = ε[ind]
    # get exponent
    @. model(x, p) = p[1]*x+p[2]
    fit = curve_fit(model, vec(log10.(ε)), vec(log10.(f_ε)), [2., 2.])
    D = coef(fit)
    @show estimate_errors(fit)
    #D = linear_region(vec(log10.(ε)), vec(log10.(f_ε)))
    return D[1],vec(log10.(ε)), vec(log10.(f_ε))

end

function ball_estimate(xg,yg,bsn; precision = 1e-4)

    xg = bsn.xg; yg = bsn.yg;
    xf=xg[end];xi=xg[1];yf=yg[end];yi=yg[1];
    y_grid_res=yg[2]-yg[1]
    nx=length(xg)
    ny=length(yg)

    num_step=10
    N_u = zeros(1,num_step) # number of uncertain box
    N = zeros(1,num_step) # number of boxes
    ε = zeros(1,num_step) # resolution

    # resolution in pixels
    min_ε = 1;
    max_ε = 10;
    @show r_ε = 10 .^ range(log10(min_ε),log10(max_ε),length=num_step)
    #r_ε = (min_ε:max_ε)*y_grid_res

    for (k,eps) in enumerate(r_ε)
        Nb=0; Nu=0; μ=0; σ²=0; M₂=0;
        completed = 0;
        c_ind = get_ind_in_circle(eps)
        # Find uncertain boxes
        while completed == 0

            kx = rand()*((nx-eps)-(1+eps))+1+eps
            ky = rand()*((ny-eps)-(1+eps))+1+eps

            # I = get_indices_in_range(kx, ky, eps, c_ind)
            # c = [ bsn.basin[n[1],n[2]] for n in eachrow(I)]

            Ix,Iy = get_indices_in_range(kx, ky, eps, c_ind)
            c = [ bsn.basin[n,m] for n in Ix, m in Iy]

            if length(unique(c))>1
                Nu = Nu + 1
            end
            Nb += 1

            # Welford's online average estimation and variance of the estimator
            M₂ = wel_var(M₂, μ, Nu/Nb, Nb)
            μ = wel_mean(μ, Nu/Nb, Nb)
            σ² = M₂/Nb

            # Stopping criterion: variance of the estimator of the mean bellow 1e-3
            if Nu > 50 && σ² < precision
                completed = 1
                @show Nu,Nb,σ²
            end

        end
        N_u[k]=Nu
        N[k]=Nb
        ε[k]=eps*y_grid_res
    end

    # uncertain function
    f_ε = N_u./N

    # remove zeros in case there are:
    ind = f_ε .> 0.
    f_ε =  f_ε[ind]
    ε = ε[ind]
    # get exponent
    @. model(x, p) = p[1]*x+p[2]
    fit = curve_fit(model, vec(log10.(ε)), vec(log10.(f_ε)), [2., 2.])
    D = coef(fit)
    @show estimate_errors(fit)
    return D[1],vec(log10.(ε)), vec(log10.(f_ε))

end

function get_indices_in_range(kx, ky, eps,c_ind)

    # center and scale

    #indx = ceil.(Int64,(kx-eps):1.:(kx+eps))
    #indy = ceil.(Int64,(ky-eps):1.:(ky+eps))
    indx = range(ceil.(Int64,kx-eps),ceil(Int64,kx+eps)-1,step=1)
    indy = range(ceil.(Int64,ky-eps),ceil(Int64,ky+eps)-1,step=1)

    #I = deepcopy(c_ind)

    #Find indices in a circle of radius eps
    #I[:,1] .= I[:,1] .+ ceil(Int64,kx)
    #I[:,2] .= I[:,2] .+ ceil(Int64,ky)
    #@show [[n[1],n[2]] for n in I]

    return indx,indy
    #return I
end

function get_ind_in_circle(r)
    #Find indices in a circle of radius eps
    I = Array{Int64,1}[]
    r_l = ceil(Int16,r)
    for n in -r_l:r_l, m in -r_l:r_l
        if n^2+m^2 ≤ r^2
            push!(I, [n; m])
        end
    end
    v= hcat(I...)'
    return v
end


text:
let

    roadway = gen_straight_roadway(3, 1000.0, lane_width=1.0)
    rec = SceneRecord(1, 0.1, 5)
    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2( 0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
        ]))

    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 1)) == 0.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 2)) == 0.0

    @test isapprox(convert(Float64, get(MARKERDIST_LEFT, rec, roadway, 1)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT, rec, roadway, 2)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, 1)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, 2)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT_LEFT, rec, roadway, 1)), 1.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT_LEFT, rec, roadway, 2)), 1.5)
    @test isnan(convert(Float64, get(MARKERDIST_RIGHT_RIGHT, rec, roadway, 1)))
    @test isnan(convert(Float64, get(MARKERDIST_RIGHT_RIGHT, rec, roadway, 2)))

    @test isapprox(convert(Float64, get(DIST_FRONT, rec, roadway, 1)), 10.0 - 5.0)
    @test convert(Float64, get(DIST_FRONT, rec, roadway, 2, censor_hi=100.0)) == 100.0

    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2( 1.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            Vehicle(VehicleState(VecSE2(12.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 3),
            Vehicle(VehicleState(VecSE2( 0.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 4),
        ]))

    @test isapprox(convert(Float64, get(MARKERDIST_LEFT, rec, roadway, 3)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT, rec, roadway, 3)), 0.5)
    @test isapprox(convert(Float64, get(MARKERDIST_LEFT_LEFT, rec, roadway, 3)), 1.5)
    @test isapprox(convert(Float64, get(MARKERDIST_RIGHT_RIGHT, rec, roadway, 3)), 1.5)

    @test isapprox(convert(Float64, get(DIST_FRONT, rec, roadway, 1)), 9.0 - 5.0)
    @test isapprox(convert(Float64, get(DIST_FRONT_LEFT, rec, roadway, 1)), 11.0 - 5.0)
    @test convert(Float64, get(DIST_FRONT_RIGHT, rec, roadway, 1)) == 100.0
    @test isapprox(convert(Float64, get(DIST_FRONT, rec, roadway, 4)), 12.0 - 5.0)
    @test convert(Float64, get(DIST_FRONT_LEFT, rec, roadway, 4)) == 100.0
    @test isapprox(convert(Float64, get(DIST_FRONT_RIGHT, rec, roadway, 4)), 1.0 - 5.0)

    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 1)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 2)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 3)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 4)) == 1.0

    update!(rec, Scene([
            Vehicle(VehicleState(VecSE2(  0.0,0.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 1),
            Vehicle(VehicleState(VecSE2(-10.0,1.0,0.0), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 2),
            Vehicle(VehicleState(VecSE2(  0.0,1.0,0.5), roadway, 10.0), VehicleDef(AgentClass.CAR, 5.0, 2.0), 3),
        ]))

    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 1)) == 1.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 2)) == 0.0
    @test convert(Float64, get(IS_COLLIDING, rec, roadway, 3)) == 1.0

    # super-simple test
    for f in allfeatures()
        get(f, rec, roadway, 3)
    end
end


text:
using MatrixAlgebra
using BandedMatrices
using MiscUtilities
using DoubleFloats
using Plots

n = 10000

function sband(n::Integer, a::Vector{T}) where T
    k = length(a)
    SH = T<:Complex ? Hermitian : Symmetric
    SH(BandedMatrix([-i => fill(a[i+1], n-i) for i = 0:min(k,n)-1]...), :L)
end

#A, B = sband(n, [6.0, -4, 1, 0.1]), sband(n, [2.0, 1, 0, 0])
#ws = MatrixAlgebra.make_workspace(A, B)
cws(T, A=A, B=B) = MatrixAlgebra.make_workspace(copy_elementtype(T, A), copy_elementtype(T, B))

eig(k) = MatrixAlgebra.eigval!(ws, k, lb[k], ub[k], 1)

function pl(a::T, b::T, v::AbstractVector = T[]) where T
    pf(x) = clamp(x - f(x), a, b)
    ph1(x) = clamp(x - h(x, 1), a, b)
    phh0(x) = clamp(x - hh(x, 1, (1.0, 0.0, 0.0, 1.0)), a, b)
    phh1(x) = clamp(x - hh(x, 1, (-1.0, 1.0, 1.0, 1.0)), a, b)
    ph2(x) = clamp(x - h(x, 2), a, b)
    ph3(x) = clamp(x - h(x, 3), a, b)
    ph10(x) = clamp(x - h(x, 10), a, b)
    ph100(x) = clamp(x - h(x, 100), a, b)
    ph1000(x) = clamp(x - h(x, 1000), a, b)
    pg1(x) = clamp(x - g(x, v), a, b)
    p = plot(legend=nothing)
    plot!(p, [a; b], [a; b])
    plot!(p, range(a,stop=b,length=2001), [pf; phh0; phh1])
    #display(p)
end

function pl(a::Integer,b::Integer = a + 1)
    d = (eve[b] - eve[a] ) * 0.5
    pl(eve[a] - d, eve[b] + d)
end

function plfh(a, b)
    p = plot(legend=nothing)
    plot!(p, [a; b], [0; 0])
    plot!(p, range(a,stop=b,length=1001), [f; h])
    #display(p)
end

function h(x::AbstractFloat, r=1)
    κ, η, ζ = detderivates(A, B, x)
    n = size(A, 1)
    MatrixAlgebra.laguerre1(η, ζ, n, r)
end

function g(x::AbstractFloat, v::AbstractVector)
    n = size(A, 1)
    p = length(v)
    y = sum(inv.(x .- v))
    z = sum(inv.(x .- v) .^ 2)
    κ, η, ζ = detderivates(A, B, x)
    r1 = float(1)
    r2 = clamp(η^2 / ζ, 1, n)
    w = 0 # η^2 / ( ζ + η^2 )
    r = r1 + ( r2 - r1 ) * w^n
    r = clamp(rest(x), 1, n)
    MatrixAlgebra.laguerre1(η - y, ζ - z, n - p, r)
end

function hh(x::AbstractFloat, r = 1, (a,b,c,d) = (1.0, 0.0, 0.0, 1.0), A=A, B=B)
    AA = A * a + B * c
    BB = B * d + A * b
    ws = cws(typeof(x), AA, BB)
    t = (a * x + c) / (d + b * x)
    κ, η, ζ = MatrixAlgebra.detderivates!(ws, t)
    n = size(A, 1)
    α = ζ / η^2
    dt = n / η / ( 1 + sqrt((n-r)/r*max(n * α - 1, 0.0)))
    tt = t - dt
    xx = (d * tt - c) / (a - b * tt)
    dx = x - xx
    dx
end

function f(x)
    κ, η, ζ = detderivates(A, B, x)
    1 / η
end

function kappa(x, A=A, B=B)
    κ, η, ζ = detderivates(A, B, x)
    κ
end

function rest(x, A=A, B=B)
    κ, η, ζ = detderivates(A, B, x)
    η ^ 2 / ζ
end

nothing


text:
@testset "393.utf-8-validation.jl" begin
    @test valid_utf8([197, 130, 1]) == true
    @test valid_utf8([235, 140, 4]) == false
    @test valid_utf8([255]) == false
    @test valid_utf8([237]) == false
    @test valid_utf8([0]) == true
end


text:
"""
    drawmandelbrotR2(f [, x0, y0; hasescaped, maxiterations])

Return the drawing of the Mandelbrot set of a family of functions
    \$f_{a,b}:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$ with \$a,b\\in\\mathbb{R}\$,
    in a rectangular region in \$\\mathbb{R}^2\$,
    using the escape time of iterations algorithm,
    over a given "critical" point \$(x_0,y_0)\$.

The Mandelbrot set of \$f_{a,b}\$ over \$(x_0,y_0)\$ is defined as
\$\\mathcal{M}(f_{a,b},(x_0,y_0))=\\{(a,b)\\in\\mathbb{R}^2\\,|\\,|f_{a,b}^n(x_0,y_0)|\\nrightarrow\\infty\\,n\\rightarrow\\infty\\}\$

#### Arguments
- `f::Function`: A function \$f:\\mathbb{R}^2\\times\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$.
- `x0::Real`: X coordinate of the "critical" point.
- `y0::Real`: Y coordinate of the "critical" point.
- `hasescaped::Function`: A boolean function to check if the iterations has escaped.
- `maxiterations::Integer`: Maximum number of iterations to check.
"""
function drawmandelbrotR2(f::Function, x0::Real=0, y0::Real=0;
    hasescaped::Function=(x::Real,y::Real) -> x*x+y*y > 4, maxiterations::Int=100)

    # Veryfying if graphics backend supports functions
    SDDGraphics.supported(:drawpixel)

    # Verifying functions
    @assert typeof(f(0., 0., 1., 1.)) <: Tuple{Real,Real}
    @assert typeof(hasescaped(1., 1.)) <: Bool

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(maxiterations)

    @sweeprectregion SDDGraphics.xlims() SDDGraphics.ylims() SDDGraphics.canvassize() begin
        xn,yn = x0,y0
        escapetime = maxiterations
        for n in 0:maxiterations
            if hasescaped(xn,yn)
                escapetime = n
                break
            end # if hasescaped
            xn,yn = f(x,y,xn,yn)
        end # for n maxiterations
        SDDGraphics.color(escapetime)
        SDDGraphics.drawpixel(i,j)
    end # Implemented algorithm

    SDDGraphics.drawing()

end # function drawtrappedpointsR2

"""
    drawmandelbrotC(f [, z0; hasescaped, maxiterations])

Return the drawing of the Mandelbrot set of a family of functions
    \$f_{\\lambda}:\\mathbb{C}\\rightarrow\\mathbb{C}\$,
    in a rectangular region in \$\\mathbb{C}\$,
    using the escape time of iterations algorithm,
    over a given "critical" point \$z_0\$.

The Mandelbrot set of \$f_{\\lambda}\$ over \$z_0\$ is defined as
\$\\mathcal{M}(f_{\\lambda},z_0)=\\{\\lambda\\in\\mathbb{C}\\,|\\,|f_{\\lambda}^n(z_0)|\\nrightarrow\\infty\\,n\\rightarrow\\infty\\}\$

#### Arguments
- `f::Function`: A function \$f:\\mathbb{C}\\times\\mathbb{C}\\rightarrow\\mathbb{C}\$.
- `z0::Number`: A "critical" point.
- `hasescaped::Function`: A boolean function to check if the iterations has escaped.
- `maxiterations::Integer`: Maximum number of iterations to check.
"""
function drawmandelbrotC(f::Function, z0::Number=0;
    hasescaped::Function=z::Number -> abs2(z) > 4, maxiterations::Int=100)

    # Veryfying if graphics backend supports functions
    SDDGraphics.supported(:drawpixel)

    # Verifying functions
    @assert typeof(f(0, 1.0im)) <: Number
    @assert typeof(hasescaped(1.0im)) <: Bool

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(maxiterations)

    @sweeprectregion SDDGraphics.xlims() SDDGraphics.ylims() SDDGraphics.canvassize() begin
        λ = complex(x,y)
        zn = z0
        escapetime = maxiterations
        for n in 0:maxiterations
            if hasescaped(zn)
                escapetime = n
                break
            end # if hasescaped
            zn = f(λ, zn)
        end # for n maxiterations
        SDDGraphics.color(escapetime)
        SDDGraphics.drawpixel(i,j)
    end # Implemented algorithm

    SDDGraphics.drawing()

end # function drawtrappedpointsC


text:
#handle scope
N = 3
f1(a, b) = a+b
g1(a, b) = (a*b, a-b)
h1(a, b) = a^b
k(a, b, c) = (a + b) / c

@testset "NNTopo" begin
    using Transformers.Stacks: print_topo
    f(x) = x+1
    g(x) = x+2
    h(x) = x+3

    topo1 = @nntopo x => a => b => y
    @test topo1((f,g,h), 10) == h(g(f(10)))

    topo2 = @nntopo x => 4 => y
    @test topo2((f,f,f,f, g), 10) == g(f(f(f(f(10)))))

    (x1, x2, x3, x4) = [2, 3, 7, 5]
    t = f1(x1, x2)
    z1, z2 = g1(t, x3)
    w = h1(x4, z1)
    y = k(x2, z2, w)

    topo3 = @nntopo (x1, x2, x3, x4):(x1, x2) => t:(t, x3) => (z1, z2):(x4, z1) => w:(x2, z2, w) => y
    @test topo3((f1, g1, h1, k), x1, x2, x3, x4) ≈ y

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()

        print_topo(outWrite, topo3; models=(f1, g1, h1, k))
        close(outWrite)

        topo_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test topo_string == """topo_func(model, x1, x2, x3, x4)
	t = f1(x1, x2)
	(z1, z2) = g1(t, x3)
	w = h1(x4, z1)
	y = k(x2, z2, w)
	y
end
"""
    end

    topo = @nntopo((e, m, mask):e → pe:(e, pe) → t → (t:(t, m, mask) → t:(t, m, mask)) → $N:t → c)

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()

        print_topo(outWrite, topo)
        close(outWrite)

        topo_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test topo_string == """topo_func(model, e, m, mask)
	pe = model[1](e)
	t = model[2](e, pe)
	t = model[3](t)
	t = model[4](t, m, mask)
	t = model[5](t, m, mask)
	t = model[6](t, m, mask)
	c = model[7](t)
	c
end
"""
    end

    localn = 3
    topo = @nntopo_str "(e, m, mask):e → pe:(e, pe) → t → (t:(t, m, mask) → t:(t, m, mask)) → $localn:t → c"

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()

        print_topo(outWrite, topo)
        close(outWrite)

        topo_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test topo_string == """topo_func(model, e, m, mask)
	pe = model[1](e)
	t = model[2](e, pe)
	t = model[3](t)
	t = model[4](t, m, mask)
	t = model[5](t, m, mask)
	t = model[6](t, m, mask)
	c = model[7](t)
	c
end
"""
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()

        topo = @nntopo x => ((y => z => t) => 3 => w) => 2
        print_topo(outWrite, topo)
        close(outWrite)

        topo_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)

        @test topo_string == """topo_func(model, x)
	y = model[1](x)
	z = model[2](y)
	t = model[3](z)
	z = model[4](t)
	t = model[5](z)
	z = model[6](t)
	t = model[7](z)
	w = model[8](t)
	z = model[9](w)
	t = model[10](z)
	z = model[11](t)
	t = model[12](z)
	z = model[13](t)
	t = model[14](z)
	w = model[15](t)
	w
end
"""
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()

        topo = @nntopo x => y' => 3 => z
        print_topo(outWrite, topo)
        close(outWrite)

        topo_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)
        if v"1.0.0" <= VERSION  <= v"1.2.0"
            @test topo_string == """topo_func(model, x)
	y = model[1](x)
	%1 = y
	y = model[2](y)
	%2 = y
	y = model[3](y)
	%3 = y
	y = model[4](y)
	%4 = y
	z = model[5](y)
	(z, (%1, %2, %3, %4))
end
"""
        else
            @test topo_string == """topo_func(model, x)
	y = model[1](x)
	%1 = y
	y = model[2](y)
	%2 = y
	y = model[3](y)
	%3 = y
	y = model[4](y)
	%4 = y
	z = model[5](y)
	(z, (var"%1", var"%2", var"%3", var"%4"))
end
"""
        end
    end

    let STDOUT = stdout
        (outRead, outWrite) = redirect_stdout()

        topo = @nntopo (x,y) => (a,b,c,d') => (w',r',y) => (m,n)' => z
        print_topo(outWrite, topo)
        close(outWrite)

        topo_string = String(readavailable(outRead))
        close(outRead)

        redirect_stdout(STDOUT)
        if v"1.0.0" <= VERSION  <= v"1.2.0"
            @test topo_string == """topo_func(model, x, y)
	(a, b, c, d) = model[1](x, y)
	%1 = d
	(w, r, y) = model[2](a, b, c, d)
	%2 = (w, r)
	(m, n) = model[3](w, r, y)
	%3 = (m, n)
	z = model[4](m, n)
	(z, (%1, %2, %3))
end
"""
        else
            @test topo_string == """topo_func(model, x, y)
	(a, b, c, d) = model[1](x, y)
	%1 = d
	(w, r, y) = model[2](a, b, c, d)
	%2 = (w, r)
	(m, n) = model[3](w, r, y)
	%3 = (m, n)
	z = model[4](m, n)
	(z, (var"%1", var"%2", var"%3"))
end
"""
        end
    end
end


text:
module BayesOpt

export bayesopt

if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    include("../deps/deps.jl")
else
    error("BayesOpt was not properly installed. Please run Pkg.build(\"BayesOpt\")")
end

include("params.jl")
include("opt.jl")

# This doesn't work on v0.4 since closures aren't c-callable functions.
function bayesopt(
        f::Function,
        x0::Vector{Float64},
        lb::Vector{Float64},
        ub::Vector{Float64};
        n_iterations::Integer=190,
        verbose_level::Integer=1)
    # wrap the function f (assumed to take an argument of Vector{Float64}) with a closure
    function f′(n, x, grad, fdata)
        return f(unsafe_wrap(Array, x, n))
    end
    params = initialize_parameters_to_default()
    params.n_iterations = n_iterations
    params.verbose_level = verbose_level
    minf = Ref{Float64}(0)
    x = copy(x0)
    r = bayes_optimization(length(x), f′, C_NULL, lb, ub, x, minf, params)
    if r != 0
        error("optimization failed")
    end
    return minf[], x
end

end # module


text:
code = ARGS[1]
code = replace(code, "⏎" => "\n")

highlighter = nothing
highlighter_str = get(ENV, "_INTERACTIVECODESEARCH_JL_HIGHLIGHTER", nothing)
if highlighter_str !== nothing
    highlighter = @eval @cmd $highlighter_str
end

file = line = nothing
if (m = match(r"(.*) in \w+ at (.*):([0-9]+)$", code)) !== nothing
    # code = String(m[1])
    file = m[2]
    line = parse(Int, m[3])
    if isabspath(file)
        print(basename(file), " (")
        printstyled(file; color = :light_black)
        println(")")
    else
        print(file, " ")
    end
    println("at line ", line)
    if isfile(file)
        open(pipeline(highlighter, stdin = IOBuffer(read(file)), stderr = stderr)) do io
            width, height = displaysize(stdout)
            for (i, str) in enumerate(eachline(io))
                if i > line + height * 2
                    close(io)
                    break
                elseif i > line - 2
                    if i == line
                        printstyled(lpad(i, 5), " "; color = :magenta, bold = true)
                        printstyled(">"; color = :red)
                    else
                        print(lpad(i, 5), " ")
                        printstyled(":"; color = :light_black)
                    end
                    print(str)
                    printstyled("\u200b")  # zero-width space
                    println()
                end
            end
        end
        exit()
    end
end

if highlighter === nothing
    print(code)
else
    run(pipeline(highlighter, stdin = IOBuffer(code), stdout = stdout, stderr = stderr))
end


text:
mutable struct Config
    train::Bool
end
const CONFIGS = [Config(true) for i=1:nthreads()]

getconfig() = CONFIGS[threadid()]
setconfig(config::Config) = CONFIGS[threadid()] = config

istrain() = getconfig().train
settrain(b::Bool) = getconfig().train = b


text:
using YaoCompiler
using YaoCompiler.Intrinsics
target = YaoTargetQASM.OpenQASMTarget()
interp = YaoInterpreter(;target)
code_typed(Intrinsics.apply, Tuple{AnyReg, typeof(Adder3.adder3())}; interp)
ast = code_qasm(typeof(Adder3.adder3()))

Adder3.t()

qelib = joinpath(pkgdir(YaoTargetQASM), "test", "qelib1.inc")
@macroexpand qasm"""
OPENQASM 2.0;
include "$qelib";
gate majority a,b,c
{
  cx c,b;
  cx c,a;
  ccx a,b,c;
}
"""

using OpenQASM
using YaoTargetQASM
using YaoCompiler
using YaoCompiler.Intrinsics
using CompilerPluginTools

qasm"""
OPENQASM 2.0;
gate cx c,t { CX c,t; }
gate u1(lambda) q { U(0,0,lambda) q; }
gate t a { u1(pi/4) a; }
gate h a { u2(0,pi) a; }
gate tdg a { u1(-pi/4) a; }

gate ccx a,b,c
{
  h c;
  cx b,c; tdg c;
  cx a,c; t c;
  cx b,c; tdg c;
  cx a,c; t b; t c; h c;
  cx a,b; t a; tdg b;
  cx a,b;
}

gate majority a,b,c
{
  cx c,b;
  cx c,a;
  ccx a,b,c;
}
"""

ast = code_qasm(typeof(majority()))

@device function test_inline()
    1:3 => ccx()
end

ast = code_qasm(typeof(test_inline()))
code_qasm(typeof(ccx()))

code_typed(Intrinsics.apply, (AnyReg, typeof(majority())); interp=YaoInterpreter())

interp = YaoInterpreter()
mi = method_instance(Intrinsics.apply, (AnyReg, typeof(test_inline())))
result = Core.Compiler.InferenceResult(mi)
frame = Core.Compiler.InferenceState(result, false, interp)
Core.Compiler.typeinf(interp, frame)

opt_params = Core.Compiler.OptimizationParams(interp)
opt = Core.Compiler.OptimizationState(frame, opt_params, interp)
preserve_coverage = Core.Compiler.coverage_enabled(opt.mod)
ci = opt.src; nargs = opt.nargs - 1;
ir = Core.Compiler.convert_to_ircode(ci, Core.Compiler.copy_exprargs(ci.code), preserve_coverage, nargs, opt)
ir = Core.Compiler.slot2reg(ir, ci, nargs, opt)
ir = compact!(ir)
sv = opt
todo = Core.Compiler.assemble_inline_todo!(ir, sv.inlining)
state = sv.inlining
todo = Core.Compiler.Pair{Int, Any}[]
et = state.et
method_table = state.method_table
ir
r = Core.Compiler.process_simple!(ir, todo, 6, state)
ir = Core.Compiler.ssa_inlining_pass!(ir, ir.linetable, sv.inlining, sv.src.propagate_inbounds)

# ir = compact!(ir)
# ir = Core.Compiler.getfield_elim_pass!(ir)
# ir = Core.Compiler.adce_pass!(ir)
# ir = Core.Compiler.type_lift_pass!(ir)
# ir = compact!(ir)
# ir = inline_const!(ir)
# ir = YaoCompiler.elim_location_mapping!(ir)

# code_typed(Intrinsics.apply, (AnyReg, typeof(ccx()), typeof(Locations(1:3))); interp)
# op = tdg()
# todo = assemble_inline_todo!(ir, state)
# ir


text:
# julia1_2.6.jl String Operations and Pattern matching
#
# String character addressing 1 to end
# The first byte of a UTF-8 multibyte character
# will yield the character, the other bytes cause exceptions
s="Abcπ∑z\u2200\n" # or s=string("Abcπ∑z∀")
# ∀ for each operator (eg. ∀x∊S for each x belonging to set S)

println(sizeof(s)) # total number of bytes in string
println(length(s)) # length in characters
# best way to deal with each character
for c in s
    print(c)
end

# and this can become a one liner using ";"
for c in s; print(c); end

println(s[5]) # this will cause an exception since continuation byte

name="Mike"
n=76
s=string("Hi ",name,n,30,"  ",π,",\n") # Converts and concatenates

# julia uses * for concatenation
#  rather than + (+ implies communtative)
s="Hi "*name*",/n"

s="hi $name,\n"   # so using String is not necessary
# The $ operator can also be used to evaluate expressions
n=63
s="Hi $(sqrt(n+1)),\n"

# triple quote used for embedded quotes and indentation
lem="""
         The lead lemming said, "Jump!"
         and they all did.
    """
println(lem)

for i=0x2200:0x2240; print(Char(i)); end # interesting symbols

###

t='a'*'z'
u='z'-'a'
v='a'+1
w='a'+'b' # meaningless so error

# character comparisons
'A'<='X'<='Z'   # true comparison
'A'<='a'<='Z'   #false comparison

# lexicographic ordering using < > <= >= == !=
"abacus"<"zoo"   # true
"Zoo"<"abacus"   # true because 'Z'<'a'Char'
"Zoo"<"Abacus"   # false
"α"<"𝛂"          # true because first alpha from lower greek alphabet.

s="πAbc∑z\u2200\n"
t=chomp(s)
chop(t)
lowercase(t)
uppercase(t)
t[end]
first(t)  # all strings start at 1.
Last(t)
split(t,'∑')
strip("   "*t*"   ")
textwidth(s)
i=1   # all strings start at 1.
t[i]
i=nextind(t,i)
t[i]
textwidth(s)

###

# regular expressions (type Regex)- Perl compatible
# see www.pcre.org
# Regex common escape sequences - table
# \a   An alarm bell
# \c?  The control character you specify for ?
# \d   A digit between 0 and 9
# \D   A non digit character
# \e   The character generated by pressing escape
# \E   The end of an \L or \U sequence
# \f   A form feed
# \l   The next lower case character
# \L   The next lower case characters until \E
# \n   A new line of data
# \Q   The next metacharacter until \E is found
# \r   a carriage return
# \s   A whitespace character
# \S   A non-white space character
# \t   A tab
# \u   The next Uppercase character
# \U   The next uppercase characters until \E
# \w   A word character (alphanumeric characters or underscores)
# \W   A non-word characters
# \O?  The octal character you specify for ?
# \x?  The hexadecimal character you specify for ?

Other metacharacters
# ^    Match at beginning of String
# $    Match at end of String
# +    1 or more of preceding type
# *    0 or more of preceding type
# [?-?]  match range, eg. [A-G], [1-5] and [A-Za-z]

# Test to make sure password not all alpha
i=0
while(i==0)
    global i
    passwd=readline()
    if match(r"^[A-Za-zα-ω]*$",passwd) === nothing
        println("password good")
        i=1
    else
        println("all alpha, try again")
    end
end

# pick time out of text
linewithtime="a good holiday weekend, its 15:45 and sunny."
m=match(r"(?<hour>\d+):(?<minute>\d+)",linewithtime)
println(m[:hour]*":"*m[:minute])

# For more read Redex Cheat Sheet at:
#   www.rexegg.com/regex-quickstart.html


text:
export MaximalIndependentSet

struct MaximalIndependentSet end

"""
    independent_set(g, MaximalIndependentSet(); seed=-1)

Find a random set of vertices that are independent (no two vertices are adjacent to each other) and
it is not possible to insert a vertex into the set without sacrificing the independence property.

### Implementation Notes
Performs [Approximate Maximum Independent Set](https://en.wikipedia.org/wiki/Maximal_independent_set#Sequential_algorithm) once.
Returns a vector of vertices representing the vertices in the independent set.

### Performance
Runtime: O(|V|+|E|)
Memory: O(|V|)
Approximation Factor: maximum(degree(g))+1

### Optional Arguments
- If `seed >= 0`, a random generator is seeded with this value.
"""
function independent_set(
    g::AbstractGraph{T},
    alg::MaximalIndependentSet;
    seed::Int=-1
    ) where T <: Integer

    nvg = nv(g)
    ind_set = Vector{T}()
    sizehint!(ind_set, nvg)
    deleted = falses(nvg)

    for v in randperm(getRNG(seed), nvg)
        (deleted[v] || has_edge(g, v, v)) && continue

        deleted[v] = true
        push!(ind_set, v)
        deleted[neighbors(g, v)] .= true
    end

    return ind_set
end


text:
## Graph layouts ##

"""
    layout_random(g::Graph)

Place nodes at random positions. They are uniformly distributed along each axis.
"""
layout_random(g::Graph) = (rand(nodecount(g)), rand(nodecount(g)))

"""
    layout_line(g::Graph)

Place nodes on a horisontal line.
Combining this layout with `edge_curve=1` when plotting yields an
[arc diagram](https://en.wikipedia.org/wiki/Arc_diagram).
"""
layout_line(g::Graph) = (
        [Float64(i) for i = nodes(g)], ones(Float64, nodecount(g)))

"""
    layout_circle(g::Graph[, clockwise=true])

Place nodes on a circle, a.k.a. [circular layout](https://en.wikipedia.org/wiki/Circular_layout).
Set `clockwise` to change node ordering.
"""
function layout_circle(g::Graph, clockwise=true)
    Δα = 2 * pi / nodecount(g) * (clockwise ? 1 : -1)
    α0 = -pi / 2 - Δα
    return (
        Float64[cos(α0 + Δα * i) for i = nodes(g)],
        Float64[sin(α0 + Δα * i) for i = nodes(g)])
end

"""
    layout_fruchterman_reingold(g::Graph[; maxiter, scale...)

Layout network with Fruchterman-Reingold force-directed algorithm.

# References

Fruchterman, Thomas MJ, and Edward M Reingold. 1991.
“Graph Drawing by Force-Directed Placement.” Software: Practice and Experience 21 (11):1129–64.
"""
function layout_fruchterman_reingold(g::Graph; maxiter=250, scale=sqrt(nodecount(g)), init_temp=sqrt(nodecount(g)))
    n = nodecount(g)
    x = scale / 2 .* (rand(n) .- 0.5)
    y = scale / 2 .* (rand(n) .- 0.5)
    @inbounds for iter = 1:maxiter
        force_x, force_y = zeros(n), zeros(n)
        for i = 1:n
            @inbounds for j = 1:i-1
                d_x, d_y = x[j] - x[i], y[j] - y[i]
                f = - 1 / (d_x^2 + d_y^2)           # repulsive force
                force_x[i] += d_x * f
                force_y[i] += d_y * f
                force_x[j] -= d_x * f
                force_y[j] -= d_y * f
            end
            @inbounds for j = map(x->x.node, g.nodes[i].forward)  # TODO: create generic way to handle this
                d_x, d_y = x[j] - x[i], y[j] - y[i]
                f = n * sqrt(d_x ^ 2 + d_y ^ 2)     # attractive force
                force_x[i] += d_x * f
                force_y[i] += d_y * f
                force_x[j] -= d_x * f
                force_y[j] -= d_y * f
            end
        end
        t = init_temp / iter
        @inbounds for i = 1:n
            force_mag = sqrt(force_x[i] ^ 2 + force_y[i] ^ 2)  # apply forces
            if force_mag > t
                coef = t / force_mag
                x[i] += force_x[i] * coef
                y[i] += force_y[i] * coef
            end
            mag = sqrt(x[i] ^ 2 + y[i] ^ 2)         # don't let points run away
            if mag > scale
                x[i] *= scale / mag
                y[i] *= scale / mag
            end
        end
    end
    return x, y
end

"""
    rescale(v, newmax, margin)

Rescale elements of `v` to fill range `margin:newmax - margin`.
"""
function rescale(v, newmax::Real, margin::Real) :: Vector{Int}
    if length(v) == 0
        return Int[]
    end
    oldmin, oldmax = extrema(v)
    if oldmax > oldmin
        k = (newmax - 2margin) / (oldmax - oldmin)
        return Int[round(Int, margin + (x - oldmin) * k) for x = v]
    else
        return Int[round(Int, newmax / 2) for x = v]
    end
end


text:
"""
    read_scalar(f::JLDFile, rr, header_offset::RelOffset)

Read raw data representing a scalar with read representation `rr` from the current position
of JLDFile `f`. `header_offset` is the RelOffset of the object header, used to resolve
cycles.
"""
function read_scalar end

"""
    read_array!(v::Array, f::JLDFile, rr)

Fill the array `v` with the contents of JLDFile `f` at the current position, assuming a
ReadRepresentation `rr`.
"""
function read_array! end

"""
    read_compressed_array!(v::Array, f::JLDFile, rr, data_length::Int)

Fill the array `v` with the compressed contents of JLDFile `f` at the current position,
assuming a ReadRepresentation `rr` and that the compressed data has length `data_length`.
"""
function read_compressed_array! end

#
# MmapIO
#

# Cutoff for using ordinary IO instead of copying into mmapped region
const MMAP_CUTOFF = 1048576

@inline function read_scalar(f::JLDFile{MmapIO}, rr, header_offset::RelOffset)
    io = f.io
    inptr = io.curptr
    obj = jlconvert(rr, f, inptr, header_offset)
    io.curptr = inptr + odr_sizeof(rr)
    obj
end

@inline function read_array!(v::Array{T}, f::JLDFile{MmapIO},
                             rr::ReadRepresentation{T,T}) where T
    io = f.io
    inptr = io.curptr
    n = length(v)
    nb = odr_sizeof(T)*n
    if nb > MMAP_CUTOFF && (!Sys.iswindows() || !f.written)
        # It turns out that regular IO is faster here (at least on OS X), but on Windows,
        # we shouldn't use ordinary IO to read, since coherency with the memory map is not
        # guaranteed
        mmapio = f.io
        regulario = mmapio.f
        seek(regulario, inptr - io.startptr)
        unsafe_read(regulario, pointer(v), nb)
    else
        unsafe_copyto!(pointer(v), convert(Ptr{T}, inptr), n)
    end
    io.curptr = inptr + odr_sizeof(T) * n
    v
end

@inline function read_array!(v::Array{T}, f::JLDFile{MmapIO},
                             rr::ReadRepresentation{T,RR}) where {T,RR}
    io = f.io
    inptr = io.curptr
    n = length(v)
    @simd for i = 1:n
        if !jlconvert_canbeuninitialized(rr) || jlconvert_isinitialized(rr, inptr)
            @inbounds v[i] = jlconvert(rr, f, inptr, NULL_REFERENCE)
        end
        inptr += odr_sizeof(RR)
    end
    io.curptr = inptr + odr_sizeof(RR) * n
    v
end

@inline function read_compressed_array!(v::Array{T}, f::JLDFile{MmapIO},
                                        rr::ReadRepresentation{T,RR},
                                        data_length::Int) where {T,RR}
    io = f.io
    inptr = io.curptr
    data = transcode(ZlibDecompressor, unsafe_wrap(Array, Ptr{UInt8}(inptr), data_length))
    @simd for i = 1:length(v)
        dataptr = Ptr{Cvoid}(pointer(data, odr_sizeof(RR)*(i-1)+1))
        if !jlconvert_canbeuninitialized(rr) || jlconvert_isinitialized(rr, dataptr)
            @inbounds v[i] = jlconvert(rr, f, dataptr, NULL_REFERENCE)
        end
    end
    io.curptr = inptr + data_length
    v
end

function write_data(io::MmapIO, f::JLDFile, data, odr::S, ::ReferenceFree,
                    wsession::JLDWriteSession) where S
    io = f.io
    ensureroom(io, odr_sizeof(odr))
    cp = io.curptr
    h5convert!(cp, odr, f, data, wsession)
    io.curptr == cp || throw(InternalError())
    io.curptr = cp + odr_sizeof(odr)
    nothing
end

function write_data(io::MmapIO, f::JLDFile, data, odr::S, ::HasReferences,
                    wsession::JLDWriteSession) where S
    io = f.io
    ensureroom(io, odr_sizeof(odr))
    p = position(io)
    cp = IndirectPointer(io, p)
    h5convert!(cp, odr, f, data, wsession)
    seek(io, p + odr_sizeof(odr))
    nothing
end

@static if Sys.isunix()
    function raw_write(io::MmapIO, ptr::Ptr{UInt8}, nb::Int)
        if nb > MMAP_CUTOFF
            pos = position(io)

            # Ensure that the current page has been flushed to disk
            msync(io, pos, min(io.endptr - io.curptr, nb))

            # Write to the underlying IOStream
            regulario = io.f
            seek(regulario, pos)
            unsafe_write(regulario, ptr, nb)

            # Invalidate cache of any pages that were just written to
            msync(io, pos, min(io.n - pos, nb), true)

            # Make sure the mapping is encompasses the written data
            ensureroom(io, nb + 1)

            # Seek to the place we just wrote
            seek(io, pos + nb)
        else
            unsafe_write(io, ptr, nb)
        end
        nothing
    end
else
    # Don't use ordinary IO to write files on Windows, since coherency with memory map is
    # not guaranteed
    function raw_write(io::MmapIO, ptr::Ptr{UInt8}, nb::Int)
        unsafe_write(io, ptr, nb)
        nothing
    end
end

write_data(io::MmapIO, f::JLDFile, data::Array{T}, odr::Type{T}, ::ReferenceFree,
           wsession::JLDWriteSession) where {T} =
    raw_write(io, Ptr{UInt8}(pointer(data)), odr_sizeof(odr) * length(data))

function write_data(io::MmapIO, f::JLDFile, data::Array{T}, odr::S, ::ReferenceFree,
                    wsession::JLDWriteSession) where {T,S}
    io = f.io
    ensureroom(io, odr_sizeof(odr) * length(data))
    cp = cporig = io.curptr
    @simd for i = 1:length(data)
        @inbounds h5convert!(cp, odr, f, data[i], wsession)
        cp += odr_sizeof(odr)
    end
    io.curptr == cporig || throw(InternalError())
    io.curptr = cp
    nothing
end

function write_data(io::MmapIO, f::JLDFile, data::Array{T}, odr::S, ::HasReferences,
                    wsession::JLDWriteSession) where {T,S}
    io = f.io
    ensureroom(io, odr_sizeof(odr) * length(data))
    p = position(io)
    cp = IndirectPointer(io, p)

    for i = 1:length(data)
        if isassigned(data, i)
            @inbounds h5convert!(cp, odr, f, data[i], wsession)
        else
            @inbounds h5convert_uninitialized!(cp, odr)
        end
        cp += odr_sizeof(odr)
    end

    seek(io, cp.offset)
    nothing
end

#
# IOStream/BufferedWriter
#

@inline function read_scalar(f::JLDFile{IOStream}, rr, header_offset::RelOffset)
    r = Vector{UInt8}(undef, odr_sizeof(rr))
    @GC.preserve r begin
        unsafe_read(f.io, pointer(r), odr_sizeof(rr))
        jlconvert(rr, f, pointer(r), header_offset)
    end
end

@inline function read_compressed_array!(v::Array{T}, f::JLDFile{IOStream},
                                        rr::ReadRepresentation{T,RR},
                                        data_length::Int) where {T,RR}
    io = f.io
    data_offset = position(io)
    n = length(v)
    data = read!(ZlibDecompressorStream(io), Vector{UInt8}(undef, odr_sizeof(RR)*n))
    @simd for i = 1:n
        dataptr = Ptr{Cvoid}(pointer(data, odr_sizeof(RR)*(i-1)+1))
        if !jlconvert_canbeuninitialized(rr) || jlconvert_isinitialized(rr, dataptr)
            @inbounds v[i] = jlconvert(rr, f, dataptr, NULL_REFERENCE)
        end
    end
    seek(io, data_offset + data_length)
    v
end

@inline function read_array!(v::Array{T}, f::JLDFile{IOStream},
                             rr::ReadRepresentation{T,T}) where T
    unsafe_read(f.io, pointer(v), odr_sizeof(T)*length(v))
    v
end

@inline function read_array!(v::Array{T}, f::JLDFile{IOStream},
                             rr::ReadRepresentation{T,RR}) where {T,RR}
    n = length(v)
    nb = odr_sizeof(RR)*n
    io = f.io
    data = read!(io, Vector{UInt8}(undef, nb))
    @simd for i = 1:n
        dataptr = Ptr{Cvoid}(pointer(data, odr_sizeof(RR)*(i-1)+1))
        if !jlconvert_canbeuninitialized(rr) || jlconvert_isinitialized(rr, dataptr)
            @inbounds v[i] = jlconvert(rr, f, dataptr, NULL_REFERENCE)
        end
    end
    v
end

function write_data(io::BufferedWriter, f::JLDFile, data, odr::S, ::DataMode,
                    wsession::JLDWriteSession) where S
    position = io.position[]
    h5convert!(Ptr{Cvoid}(pointer(io.buffer, position+1)), odr, f, data, wsession)
    io.position[] = position + odr_sizeof(odr)
    nothing
end

function write_data(io::BufferedWriter, f::JLDFile, data::Array{T}, odr::Type{T}, ::ReferenceFree,
                    wsession::JLDWriteSession) where T
    unsafe_write(io, Ptr{UInt8}(pointer(data)), odr_sizeof(odr) * length(data))
    nothing
end

function write_data(io::IOStream, f::JLDFile, data::Array{T}, odr::Type{T}, ::ReferenceFree,
                    wsession::JLDWriteSession) where T
    unsafe_write(io, Ptr{UInt8}(pointer(data)), odr_sizeof(odr) * length(data))
    nothing
end

function write_data(io::BufferedWriter, f::JLDFile, data::Array{T}, odr::S,
                    ::DataMode, wsession::JLDWriteSession) where {T,S}
    position = io.position[]
    cp = Ptr{Cvoid}(pointer(io.buffer, position+1))
    @simd for i = 1:length(data)
        if isassigned(data, i)
            @inbounds h5convert!(cp, odr, f, data[i], wsession)
        else
            @inbounds h5convert_uninitialized!(cp, odr)
        end
        cp += odr_sizeof(odr)
    end
    io.position[] = position + odr_sizeof(odr) * length(data)
    nothing
end

function write_data(io::IOStream, f::JLDFile, data::Array{T}, odr::S, wm::DataMode,
                    wsession::JLDWriteSession) where {T,S}
    nb = odr_sizeof(odr) * length(data)
    buf = Vector{UInt8}(undef, nb)
    pos = position(io)
    cp = Ptr{Cvoid}(pointer(buf))
    @simd for i = 1:length(data)
        if isassigned(data, i)
            @inbounds h5convert!(cp, odr, f, data[i], wsession)
        else
            @inbounds h5convert_uninitialized!(cp, odr)
        end
        cp += odr_sizeof(odr)
    end
    # We might seek around in the file as a consequence of writing stuff, so seek back. We
    # don't need to worry about this for a BufferedWriter, since it will seek back before
    # writing.
    !isa(wm, ReferenceFree) && seek(io, pos)
    write(io, buf)
    nothing
end


text:
"""
## module ASE

### Summary

Provides Julia wrappers for some of ASE's functionality. Currently:

* `ase.Atoms` becomes `type ASEAtoms`
* `ase.calculators.neighborlist.NeighborList`
  should become `ASENeighborList`, but it is ignored to use
  the `matscipy` neighbourlist instead.
  TODO: implement ASENeighbourList as backup in case `matscipy` is
  not available
* `Atoms.get_array` becomes `get_data`
* `Atoms.set_array` becomes `set_data!`
"""
module ASE

# the functions to be implemented
import JuLIP:
      positions, set_positions!,  unsafe_positions,  # ✓
      cell, set_cell!,             # ✓
      pbc, set_pbc!,               # ✓
      set_data!, get_data,         # ✓
      calculator, set_calculator!, # ✓
      constraint, set_constraint!, # ✓
      neighbourlist,                # ✓
      energy, forces

import Base.length, Base.deleteat!         # ✓

# from arrayconversions:
using JuLIP: mat, vecs, JVecF, JVecs, JVecsF, pyarrayref,
      AbstractAtoms, AbstractConstraint, NullConstraint,
      AbstractCalculator, NullCalculator

# extra ASE functionality:
import Base.repeat         # ✓
export ASEAtoms,      # ✓
      repeat, rnn, chemical_symbols, ASECalculator, extend!
using PyCall

@pyimport ase.io as ase_io
@pyimport ase.atoms as ase_atoms
@pyimport ase.build as ase_build

#################################################################
###  Wrapper for ASE Atoms object and its basic functionality
################################################################

"""
`type ASEAtoms <: AbstractAtoms`

Julia wrapper for the ASE `Atoms` class.

## Constructors

The default constructor is
```
at = ASEAtoms("Al"; repeat=(2,3,4), cubic=true, pbc = (true, false, true))
```

For internal usage there is also a constructor `ASEAtoms(po::PyObject)`
"""
type ASEAtoms <: AbstractAtoms
   po::PyObject       # ase.Atoms instance
   calc::AbstractCalculator
   cons::AbstractConstraint
end

ASEAtoms(po::PyObject) = ASEAtoms(po, NullCalculator(), NullConstraint())

function set_calculator!(at::ASEAtoms, calc::AbstractCalculator)
   at.calc = calc
   return at
end
calculator(at::ASEAtoms) = at.calc
function set_constraint!(at::ASEAtoms, cons::AbstractConstraint)
   at.cons = cons
   return at
end
constraint(at::ASEAtoms) = at.cons

function ASEAtoms( s::AbstractString;
                   repeatcell=nothing, cubic=false, pbc=(true,true,true) )
   at = bulk(s, cubic=cubic)
   repeatcell != nothing ? at = repeat(at, repeatcell) : nothing
   set_pbc!(at, pbc)
   return at
end

ASEAtoms(s::AbstractString, X::JVecsF) = ASEAtoms( ase_atoms.Atoms(s, mat(X)') )

"Return the PyObject associated with `a`"
pyobject(a::ASEAtoms) = a.po

get_array(a::ASEAtoms, name) = a.po[:get_array(name)]

set_array!(a::ASEAtoms, name, value) = a.po[:set_array(name, value)]

#
# TODO: write an explanation about storage layout here
#
positions(at::ASEAtoms) = copy( pyarrayref(at.po["positions"]) ) |> vecs
unsafe_positions(at::ASEAtoms) = pyarrayref(at.po["positions"]) |> vecs

function set_positions!(a::ASEAtoms, p::JVecsF)
   p_py = PyReverseDims(mat(p))
   a.po[:set_positions](p_py)
   return a
end

length(at::ASEAtoms) = length( unsafe_positions(at::ASEAtoms) )

set_pbc!(at::ASEAtoms, val::Bool) = set_pbc!(at, (val,val,val))

function set_pbc!(a::ASEAtoms, val::NTuple{3,Bool})
   a.po[:pbc] = val
   return a
end

pbc(a::ASEAtoms) = a.po[:pbc]

cell(at::ASEAtoms) = at.po[:get_cell]()

set_cell!(a::ASEAtoms, p::Matrix) = (a.po[:set_cell](p); a)

function deleteat!(at::ASEAtoms, n::Integer)
   at.po[:__delitem__](n-1) # delete in the actual array
   return at
end

"""
`repeat(a::ASEAtoms, n::(Int64, Int64, Int64)) -> ASEAtoms`

Takes an `ASEAtoms` configuration / cell and repeats is n_j times
into the j-th dimension. For example,
```
    atm = repeat( bulk("C"), (3,3,3) )
```
creates 3 x 3 x 3 unit cells of carbon.
"""
repeat(a::ASEAtoms, n::NTuple{3, Int64}) = ASEAtoms(a.po[:repeat](n))

import Base.*
*(at::ASEAtoms, n::NTuple{3, Int64}) = repeat(at, n)
*(n::NTuple{3, Int64}, at::ASEAtoms) = repeat(at, n)
*(at::ASEAtoms, n::Integer) = repeat(at, (n,n,n))
*(n::Integer, at::ASEAtoms) = repeat(at, (n,n,n))

export bulk, graphene_nanoribbon, nanotube, molecule

@doc ase_build.bulk[:__doc__] ->
bulk(args...; kwargs...) = ASEAtoms(ase_build.bulk(args...; kwargs...))

@doc ase_build.graphene_nanoribbon[:__doc__] ->
graphene_nanoribbon(args...; kwargs...) =
   ASEAtoms(ase_build.graphene_nanoribbon(args...; kwargs...))

"nanotube(n, m, length=1, bond=1.42, symbol=\"C\", verbose=False)"
nanotube(args...; kwargs...) =
      ASEAtoms(ase_build.nanotube(args...; kwargs...))

@doc ase_build.molecule[:__doc__] ->
      molecule(args...; kwargs...) =
         ASEAtoms(ase_build.molecule(args...; kwargs...))

############################################################
# matscipy neighbourlist functionality
############################################################

include("MatSciPy.jl")

neighbourlist(at::ASEAtoms, cutoff::Float64) = MatSciPy.NeighbourList(at, cutoff)

######################################################
#    Attaching an ASE-style calculator
#    ASE-style aliases
######################################################

type ASECalculator <: AbstractCalculator
   po::PyObject
end

function set_calculator!(at::ASEAtoms, calc::ASECalculator)
   at.po[:set_calculator](calc.po)
   at.calc = calc
   return at
end

set_calculator!(at::ASEAtoms, po::PyObject) =
      set_calculator!(at, ASECalculator(po))

forces(calc::ASECalculator, at::ASEAtoms) = at.po[:get_forces]()' |> vecs
energy(calc::ASECalculator, at::ASEAtoms) = at.po[:get_potential_energy]()

"""
Creates an `ASECalculator` that uses `ase.calculators.emt` to compute
energy and forces. This is very slow and is only included for
demonstration purposes.
"""
function EMTCalculator()
   @pyimport ase.calculators.emt as emt
   return ASECalculator(emt.EMT())
end

################### extra ASE functionality  ###################
# TODO: we probably want more of this
#       and a more structured way to translate
#

"""
`extend!(at::ASEAtoms, atadd::ASEAtoms)`

add `atadd` atoms to `at` and returns `at`; only `at` is modified

A short variant is
```julia
extend!(at, (s, x))
```
where `s` is a string, `x::JVecF` a position
"""
function extend!(at::ASEAtoms, atadd::ASEAtoms)
   at.po[:extend](atadd.po)
   return at
end

extend!{S <: AbstractString}(at::ASEAtoms, atnew::Tuple{S,JVecF}) =
   extend!(at, ASEAtoms(atnew[1], atnew[2]))

"return vector of chemical symbols as strings"
chemical_symbols(at::ASEAtoms) = pyobject(at)[:get_chemical_symbols]()

write(filename::AbstractString, at::ASEAtoms) = ase_io.write(filename, at.po)

# TODO: trajectory; see
#       https://wiki.fysik.dtu.dk/ase/ase/io/trajectory.html

# TODO: rnn should be generalised to compute a reasonable rnn estimate for
#       an arbitrary set of positions
"""
`rnn(species)` : returns the nearest-neighbour distance for a given species
"""
function rnn(species::AbstractString)
   X = unsafe_positions(bulk(species) * 2)
   return minimum( norm(X[n]-X[m]) for n = 1:length(X) for m = n+1:length(X) )
end

end # module